ℹ️ Skipped - page is already crawled
| Filter | Status | Condition | Details |
|---|---|---|---|
| HTTP status | PASS | download_http_code = 200 | HTTP 200 |
| Age cutoff | PASS | download_stamp > now() - 6 MONTH | 1.4 months ago |
| History drop | PASS | isNull(history_drop_reason) | No drop reason |
| Spam/ban | PASS | fh_dont_index != 1 AND ml_spam_score = 0 | ml_spam_score=0 |
| Canonical | PASS | meta_canonical IS NULL OR = '' OR = src_unparsed | Not set |
| Property | Value |
|---|---|
| URL | https://www.mdpi.com/1996-1944/19/5/964 |
| Last Crawled | 2026-03-04 05:00:49 (1 month ago) |
| First Indexed | not set |
| HTTP Status Code | 200 |
| Meta Title | Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading |
| Meta Description | To meet the requirements of effectiveness and strength in actual engineering, based on the dynamic fracture characteristics, the dynamic propagation of orthogonal anisotropic interface cracks in piezoelectric bimaterials was analyzed. By performing Laplace transformation and Fourier transformation on the governing equations, the problem was transformed into a singular integral equation. Using the Chebyshev point method and Laplace inversion, the stress and electric displacement intensity factors at the crack tip of the orthogonal anisotropic interface were obtained. The results show that the crack length affects the dimensionless function. The longer the crack, the larger the dimensionless function. Under certain conditions, the smaller the elastic parameters, the smaller the dimensionless dynamic stress intensity factor. At the same time, the impact time also affects the dynamic crack propagation. With the passage of time, the dimensionless function first increases, then reaches a peak, and finally oscillates and converges to the static value. On this basis, the response surface method was used for analysis and prediction. The R2 value of the random forest model is 0.9886, which indicates that the model has high predictive accuracy. When the optimal values of A (d1/a), B (cpt/a) and C (c44(2)/c44(1)) are 0.4045, 1.6797 and 1.9035 respectively, the stress intensity reaches its maximum value of 1.3375. |
| Meta Canonical | null |
| Boilerpipe Text | 1. Introduction
Piezoelectric materials have been widely used in advanced equipment, instruments and materials due to their electromechanical coupling properties. However, during manufacturing or usage, defects such as cracks are inevitable, which may lead to structural damage or device failure. Therefore, the crack problem of piezoelectric materials has attracted considerable attention [
1
,
2
,
3
,
4
,
5
]. In recent years, piezoelectric bimaterials have been widely applied in devices. The sudden change in materials on both sides of the interface may cause stress concentration, which is one of the important reasons for material failure [
6
,
7
,
8
]. Li et al. [
9
] analyzed the SH-wave scattering by interfacial cracks in piezoelectric ceramic-polymer composites. Their core objective was to establish a dynamic response analysis framework for cracks under a steady-state wave field, and they preliminarily investigated the effects of frequency, incident angle, and material properties on crack propagation. This work provides a fundamental wave-based analysis model for subsequent dynamic fracture research; however, it does not address the transient response under impact loading or the coupling behavior in complex anisotropic material systems.
Research on the behavior of interface crack propagation has achieved fruitful results. Loboda et al. conducted studies on the conductive interface cracks between one-dimensional hexagonal piezoelectric quasicrystals and ordinary piezoelectric materials, analyzing the coupling effect between the cracks and the far-field electrodes, and clarifying the energy-evolution law of piezoelectric hetero-interface cracks, providing theoretical references for the analysis of related issues [
10
,
11
,
12
]. Sun et al. [
13
] proposed the arbitrary-order generalized finite difference method (GFDM), providing a new numerical solution path for homogeneous and dual-material piezoelectric crack problems. The crack propagation behavior of multiple interfaces caused by circular holes under dynamic loading was studied, and the influence law of functional gradient characteristics on the dynamic response of cracks was clarified, confirming that the gradient distribution of material parameters would significantly change the stress-electric field concentration effect at the crack tip [
14
]. Based on the improved mechanical energy release rate criterion, the dynamic expansion characteristics of cracks in piezoelectric plates were revealed [
15
]. In terms of innovative analysis methods, a new interaction integral method was proposed, achieving precise analysis of dynamic impermeable interface cracks in heterogeneous piezoelectric materials, effectively improving the accuracy of field quantity calculation at the crack tip [
16
]. Liu et al. [
17
,
18
,
19
] based on the non-local piezoelectric theory, systematically studied the dynamic response of I-type penetrating cracks in piezoelectric materials, derived the dynamic analytical solutions for finite-penetrating and complete-penetrating cracks, and improved the theoretical system of dynamic fracture in piezoelectric materials. Afshar et al. [
20
] studied the dynamic fracture behavior of multiple intruded cracks in functional gradient piezoelectric strips [
20
]. Wünsche et al. [
21
] used the symmetrical Galerkin boundary element method to establish an analysis model for dynamic cracks in piezoelectric bodies under harmonic loading, providing an efficient solution method for piezoelectric fracture problems under harmonic loading. Yu et al. [
22
] based on the extended finite element method (XFEM), established a numerical solution framework for the dynamic crack of heterogeneous piezoelectric materials under force–electric coupling loading, providing an efficient calculation tool for interface fracture problems in complex configurations.
Through in-depth analysis of the existing research, we can find that the electromechanical coupling constitutive relationship of orthotropic anisotropic materials is more complex, and the direction dependence of their physical parameters leads to the inability to directly transplant the existing theoretical models for isotropic or functionally graded materials; at the same time, although some studies involve dynamic loads, they are mainly concentrated on harmonic loads or do not fully capture the transient response characteristics under impact loads; and the existing research mostly focuses on specific types of conductive cracks, non-conductive cracks, etc., and mainly focuses on the analysis of the influence of a single factor, without deeply exploring the coupling effect and quantitative influence laws of multiple parameters such as elastic parameters, crack length, and impact time, which is difficult to meet the systematic research requirements of this article on the dynamic expansion mechanism of interface cracks in orthotropic anisotropic piezoelectric dual materials.
In summary, although the existing research has laid a foundation for the crack problem of dual-layer materials, it has obvious limitations in adapting to the specific research object, load conditions, and analysis goals of this article.
The above research has constructed a research framework for the crack problem of piezoelectric materials from aspects such as material system, crack form, analysis method, and load type. However, the theory and methods still have obvious limitations in applying to the orthogonal anisotropic piezoelectric dual-material interface crack problem under dynamic loads. On one hand, the existing research on dual-material interface cracks mostly focuses on static/sinusoidal loads, isotropic/functional gradient piezoelectric media, or do not consider the orthogonal anisotropy characteristics of the materials. The electromechanical coupling constitutive relationship of orthogonal anisotropic piezoelectric dual-material is more complex, and the transient response laws under dynamic impact loads are significantly different from those of conventional media. The existing static theories and analysis methods under sinusoidal loads cannot be directly transplanted. On the other hand, most dynamic-crack research focuses on internal cracks of a single piezoelectric material or does not deeply explore the coupling effects of multiple parameters and quantitative influence rules. There is no complete analytical–numerical coupled solution system for the dynamic-fracture problem of orthogonal anisotropic piezoelectric dual-material interfaces, and the research on the quantitative influence of key parameters such as elastic parameters, crack length, and impact time and the interaction between parameters is still scarce. In addition, most existing studies adopt a single analytical or numerical method for crack analysis. The research of integrating the response-surface method into dynamic fracture analysis to achieve rapid prediction of multiple parameters and identification of optimal working conditions is still in the exploratory stage.
During manufacturing and usage, cracks may occur at the interface due to manufacturing level and other factors. Meanwhile, under dynamic loads, the material thickness, elastic constants, etc., of the piezoelectric bimaterials will affect the dynamic propagation of interface cracks, which may lead to device failure and losses. However, there are few articles on this topic. Moreover, the dynamic propagation of interface cracks in piezoelectric bimaterials has limitations in software simulation and experiments. Therefore, precise quantitative analysis of the dynamic fracture problem of interface cracks in piezoelectric bimaterials is of great significance for structural design and safety.
The motivation for this study stems from the urgent need to address the dynamic fracture behavior of interface cracks in piezoelectric bimaterials. This is a topic that has received limited attention in existing literature, especially concerning the complex scenario involving orthogonally anisotropic materials under dynamic loads. Piezoelectric materials play a critical role in advanced equipment, instruments, and sensors due to their unique electromechanical coupling properties. However, their brittle nature and inevitable defects, such as cracks, during manufacturing and use can lead to structural damage or even device failure. While research on static fracture problems has a certain foundation, many piezoelectric components in practical applications are often subjected to dynamic loads, making the study of dynamic fracture behavior particularly important and challenging. Currently, precise quantitative analysis of dynamic interface cracks in orthogonally anisotropic piezoelectric bimaterials, whether in theoretical analysis, numerical simulation, or experimental verification, faces significant deficiencies.
Our primary contribution lies in the innovative application and integration of advanced analytical and computational methods to deeply explore this complex problem. Specifically, we employed the powerful mathematical tools of Laplace transformation and Fourier transformation, successfully converting the governing equations that describe the dynamic expansion of orthogonally anisotropic interface cracks into a complex singular integral equation. Subsequently, we applied a high-precision Chebyshev point method for numerical discretization of this integral equation, combined with efficient Laplace inversion techniques, to accurately solve for the stress-intensity factors and electric-displacement-intensity factors at the crack tip. This combined analytical and numerical approach provides a solid foundation for understanding dynamic-fracture mechanisms. Furthermore, to further enhance analysis efficiency and predictive capability, we proactively integrated the Response Surface Methodology (RSM) and specifically developed and applied a high-precision Random Forest (Random Forest) machine learning model (with a determination coefficient R
2
as high as 0.9886, demonstrating the model’s excellent goodness of fit and predictive accuracy) for effective analysis and prediction of complex multi-parameter influences.
This multi-level comprehensive method offers significant added value. It not only precisely reveals the quantitative influence laws of key factors such as crack length, elastic parameters (e.g., modulus), and impact time on dynamic fracture behavior—for example, revealing the positive correlation between crack length and the dimensionless function, as well as the inverse relationship between elastic parameters and the dimensionless dynamic stress intensity factor under certain conditions—but also captures the dynamic evolution process induced by impact time, where the dynamic crack response first increases, reaches a peak, and then oscillates and converges to the static value. Moreover, it ultimately enables the identification of optimal conditions affecting stress intensity. This in-depth and quantitative understanding greatly enhances the prediction accuracy of dynamic crack propagation mechanisms, providing a key theoretical basis for evaluating and optimizing the dynamic reliability of piezoelectric devices, thereby significantly improving the operational reliability and structural safety of piezoelectric equipment in practical engineering applications. By directly addressing the limitations of current software simulations in handling dynamic fracture problems of complex anisotropic materials, as well as the challenges of experimental research in terms of cost, scale, and dynamic loading control, this study aims to fill a critical gap in the field and provide important scientific support and technical guidance for the future design of more robust, safer, and more reliable piezoelectric smart structures.
2. Establishment of the Model
As shown in
Figure 1
, the thicknesses of the two materials are
d
1
and
d
2
, respectively. The crack is from
a
1
to
b
1
. With the interface direction as the X-axis and the thickness direction as the Y-axis. Let us assume that the polarization direction of the piezoelectric dual materials is the Z-axis, the structure is orthotropic, and the symmetry axes of the materials are perfectly aligned at the interface.
The governing equation for the crack in the piezoelectric bimaterial orthotropic interface is:
c
55
(
k
)
∂
2
w
∂
X
2
+
c
44
(
k
)
∂
2
w
∂
Y
2
+
e
15
(
k
)
∂
2
ϕ
∂
X
2
+
e
24
(
k
)
∂
2
ϕ
∂
Y
2
=
ρ
∂
2
w
∂
t
2
(1)
e
15
(
k
)
∂
2
w
∂
X
2
+
e
24
(
k
)
∂
2
w
∂
Y
2
−
ε
11
(
k
)
∂
2
ϕ
∂
X
2
−
ε
22
(
k
)
∂
2
ϕ
∂
Y
2
=
0
(2)
In the above equations,
w
=
w
(
x
,
y
,
t
)
denotes the mechanical displacements,
φ
(
x
,
y
,
t
)
represents the electric potentials,
ρ
is the mass density,
c
55
(
k
)
and
c
44
(
k
)
are the elastic stiffness constants,
e
15
(
k
)
and
e
24
(
k
)
stand for the piezoelectric constants, and
ε
11
(
k
)
and
ε
22
(
k
)
are the dielectric constants. Here,
k
=
1
,
2
correspond to piezoelectric material 1 and piezoelectric material 2, respectively.
Introduce the transformation
x
=
X
−
v
t
,
y
=
Y
, then Equations (1) and (2) becomes:
c
55
(
k
)
c
44
(
k
)
−
v
2
c
k
¯
2
1
e
15
(
k
)
c
44
(
k
)
e
24
(
k
)
c
44
(
k
)
e
15
(
k
)
c
44
(
k
)
e
24
(
k
)
c
44
(
k
)
−
ε
11
(
k
)
c
44
(
k
)
−
ε
22
(
k
)
c
44
(
k
)
∂
2
w
∂
x
2
∂
2
w
∂
y
2
∂
2
ϕ
∂
x
2
∂
2
ϕ
∂
y
2
=
Ψ
(
k
)
0
(3)
Perform the Laplace transform on time
t
and the Fourier transform on
x
. The solution of Equation (3) can be expressed as
[
w
*
x
,
y
,
p
ϕ
*
x
,
y
,
p
]
=
1
2
π
∫
−
∞
+
∞
[
χ
11
(
k
)
χ
12
(
k
)
χ
21
(
k
)
χ
22
(
k
)
]
[
A
11
(
k
)
e
(
3
−
2
k
)
s
γ
1
(
k
)
y
A
12
(
k
)
e
(
3
−
2
k
)
s
γ
2
(
k
)
y
]
e
−
i
s
x
d
s
(4)
where
χ
1
η
(
k
)
=
e
24
(
k
)
c
44
(
k
)
γ
η
(
k
)
2
−
s
2
e
15
(
k
)
c
44
(
k
)
,
χ
2
η
(
k
)
=
γ
η
(
k
)
2
−
s
2
[
(
c
55
(
k
)
c
44
(
k
)
−
v
2
c
k
¯
2
)
+
2
v
p
i
c
k
¯
2
s
+
p
2
c
k
¯
2
s
2
]
,
η
=
1
,
2
.
The value of
γ
η
(
k
)
can be found in
Appendix A
.
By analyzing the model, it can be concluded that the boundary conditions that the cracks must satisfy are
σ
y
z
(
x
,
d
1
,
t
)
=
σ
y
z
(
x
,
−
d
2
,
t
)
=
0
(5)
w
(
x
,
0
+
,
t
)
=
w
(
x
,
0
−
,
t
)
,
D
y
(
x
,
0
+
,
t
)
=
D
y
(
x
,
0
−
,
t
)
(6)
σ
y
z
(
x
,
0
+
,
t
)
=
σ
y
z
(
x
,
0
−
,
t
)
=
−
τ
0
H
(
t
)
D
y
(
x
,
0
+
,
t
)
=
D
y
(
x
,
0
−
,
t
)
=
−
D
0
H
(
t
)
(7)
where
H
(
t
)
=
0
,
t
<
0
1
,
t
>
0
.
By analyzing the model, it can be concluded that the boundary conditions that the cracks must satisfy are
σ
y
z
*
(
x
,
d
1
,
p
)
=
σ
y
z
*
(
x
,
−
d
2
,
p
)
=
0
(8)
w
*
(
x
,
0
+
,
p
)
=
w
*
(
x
,
0
−
,
p
)
,
D
y
*
(
x
,
0
+
,
p
)
=
D
y
*
(
x
,
0
−
,
p
)
(9)
σ
y
z
*
(
x
,
0
+
,
p
)
=
σ
y
z
*
(
x
,
0
−
,
p
)
=
−
τ
0
p
,
D
y
*
(
x
,
0
+
,
p
)
=
D
y
*
(
x
,
0
−
,
p
)
=
−
D
0
p
(10)
The electric displacement corresponding to the stress in Equation (4) is
σ
y
z
*
D
y
*
=
1
2
π
∫
−
∞
+
∞
(
c
44
(
k
)
e
24
(
k
)
e
15
(
k
)
−
ε
22
(
k
)
χ
11
(
k
)
χ
12
(
k
)
χ
21
(
k
)
χ
22
(
k
)
γ
1
(
k
)
e
(
3
−
2
k
)
s
γ
1
(
k
)
y
A
11
(
k
)
γ
2
(
k
)
e
(
3
−
2
k
)
s
γ
2
(
k
)
y
A
12
(
k
)
(
3
−
2
k
)
s
)
e
−
i
s
x
d
s
(11)
3. Solution of the Model
In the process of crack propagation, the expansion of the crack is affected by a variety of complex factors. To more accurately describe the crack propagation behavior, introducing a density function is an effective method. Therefore, the following density function is introduced:
h
w
(
x
)
=
d
d
x
[
w
(
x
,
0
+
,
p
)
−
w
(
x
,
0
−
,
p
)
]
(12)
h
ϕ
(
x
)
=
d
d
x
[
ϕ
(
x
,
0
+
,
p
)
−
ϕ
(
x
,
0
−
,
p
)
]
(13)
Since displacement and potential are continuous, we can obtain:
∫
a
1
b
1
h
w
(
x
)
d
x
=
0
,
x
∉
(
a
1
,
b
1
)
(14)
∫
a
1
b
1
h
ϕ
(
x
)
d
x
=
0
,
x
∉
(
a
1
,
b
1
)
(15)
Based on the boundary conditions (8–10) and Equations (12) and (13), the following conclusion can be obtained:
χ
11
(
1
)
χ
12
(
1
)
χ
21
(
1
)
χ
22
(
1
)
A
11
B
11
−
χ
11
(
2
)
χ
12
(
2
)
χ
21
(
2
)
χ
22
(
2
)
A
12
B
12
=
i
s
∫
a
1
b
1
h
w
(
s
)
h
ϕ
(
s
)
e
i
s
f
d
f
=
H
1
H
2
(16)
A
11
(
1
)
A
12
(
1
)
A
11
(
2
)
A
12
(
2
)
=
1
A
*
A
1
B
1
A
2
A
2
A
1
*
B
1
*
A
2
*
B
2
*
H
1
H
2
(17)
The values of
A
1
,
B
1
,
A
2
,
B
2
,
A
1
*
,
B
1
*
,
A
2
*
,
B
2
*
,
A
*
can be found in
Appendix A
.
Substituting Equation (17) into Equation (11) yields
σ
y
z
*
(
x
,
0
+
,
p
)
D
y
*
(
x
,
0
+
,
p
)
=
lim
y
→
0
+
1
2
π
∫
−
∞
+
∞
s
e
s
γ
y
A
*
C
1
C
2
C
3
C
4
H
1
H
2
e
−
i
s
x
d
s
(18)
where
C
1
=
∑
i
=
1
2
(
c
44
(
1
)
χ
1
i
(
1
)
+
e
24
(
1
)
χ
2
i
(
1
)
)
A
i
γ
i
(
1
)
,
C
2
=
∑
i
=
1
2
(
c
44
(
1
)
χ
1
i
(
1
)
+
e
24
(
1
)
χ
2
i
(
1
)
)
B
i
γ
i
(
1
)
,
C
3
=
∑
i
=
1
2
(
e
15
(
1
)
χ
1
i
(
1
)
−
ε
22
(
1
)
χ
2
i
(
1
)
)
A
i
*
γ
i
(
1
)
,
C
4
=
∑
i
=
1
2
(
e
15
(
1
)
χ
1
i
(
1
)
−
ε
22
(
1
)
χ
2
i
(
1
)
)
B
i
*
γ
i
(
1
)
,
e
s
γ
y
=
∑
i
=
1
2
e
s
γ
i
(
1
)
y
On the crack surface, there is
∫
a
1
b
1
(
C
1
*
C
2
*
C
3
*
C
4
*
1
B
*
(
f
−
x
)
+
k
1
(
x
,
f
)
k
2
(
x
,
f
)
k
3
(
x
,
f
)
k
4
(
x
,
f
)
)
h
w
(
f
)
h
ϕ
(
f
)
d
f
=
−
π
p
σ
0
D
0
(19)
where
k
1
(
x
,
f
)
=
∫
0
∞
(
C
1
A
*
−
C
1
*
B
*
)
sin
(
s
(
f
−
x
)
)
d
s
,
k
2
(
x
,
f
)
=
∫
0
∞
(
C
2
A
*
−
C
2
*
B
*
)
sin
(
s
(
f
−
x
)
)
d
s
,
k
3
(
x
,
f
)
=
∫
0
∞
(
C
3
A
*
−
C
3
*
B
*
)
sin
(
s
(
f
−
x
)
)
d
s
,
k
4
(
x
,
f
)
=
∫
0
∞
(
C
4
A
*
−
C
4
*
B
*
)
sin
(
s
(
f
−
x
)
)
d
s
.
The value of
C
1
*
,
C
2
*
,
C
3
*
,
C
4
*
,
B
*
is referred to in
Appendix A
.
To simplify the calculation, we will transform the normalized quantity into the standard form and introduce the dislocation density function
h
w
(
x
)
=
q
1
(
x
)
,
h
ϕ
(
x
)
=
q
2
(
x
)
,
k
l
(
x
,
f
)
=
P
l
(
r
,
c
)
,
(
l
=
1
,
2
,
3
,
4
)
,
m
o
=
b
1
−
a
1
2
,
n
o
=
b
1
+
a
1
2
,
f
=
x
−
n
0
m
0
,
c
=
f
−
n
0
m
0
,
(20)
q
1
(
c
)
=
h
1
(
c
)
1
−
c
2
π
σ
0
p
,
q
2
(
c
)
=
h
2
(
c
)
1
−
c
2
π
D
0
p
.
According to Equation (20) and the Chebyshev placement method, Equation (19) can be transformed into
1
N
∑
Q
=
0
N
Φ
Q
(
C
1
*
C
2
*
C
3
*
C
4
*
Ξ
+
m
0
P
1
(
r
u
,
c
Q
)
P
2
(
r
u
,
c
Q
)
P
3
(
r
u
,
c
Q
)
P
4
(
r
u
,
c
Q
)
)
h
1
(
c
Q
)
h
2
(
c
Q
)
=
−
1
−
1
(21)
where
Ξ
=
m
0
B
*
(
m
0
c
Q
+
n
0
−
m
0
r
u
−
n
0
)
,
∑
Q
=
0
N
Φ
Q
h
1
(
c
Q
)
=
0
,
∑
Q
=
0
N
Φ
Q
h
2
(
c
Q
)
=
0
,
u
=
1
,
2
,
…
N
,
u
=
1
,
2
,
…
N
,
Φ
0
=
Φ
N
=
1
2
,
Φ
1
=
…
=
Φ
N
−
1
=
1
,
r
u
=
cos
(
2
u
−
1
)
π
2
N
,
c
Q
=
cos
Q
π
N
(N represents a node).
4. Intensity Factor
According to Equation (21), it can be obtained that
lim
x
→
a
1
−
σ
y
z
*
(
x
,
0
,
p
)
=
σ
0
C
1
*
p
B
*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
+
σ
0
C
2
*
p
B
*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(22)
lim
x
→
b
1
+
σ
y
z
*
(
x
,
0
,
p
)
=
σ
0
C
1
*
p
B
*
lim
c
→
1
+
−
h
1
(
1
)
c
2
−
1
+
σ
0
C
2
*
p
B
*
lim
c
→
1
+
−
h
2
(
1
)
c
2
−
1
(23)
lim
x
→
a
1
−
D
y
*
(
x
,
0
,
p
)
=
D
0
C
3
*
p
B
*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
+
D
0
C
4
*
p
B
*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(24)
lim
x
→
b
1
+
D
y
*
(
x
,
0
,
p
)
=
D
0
C
3
*
p
B
*
lim
c
→
1
+
−
h
1
(
1
)
c
2
−
1
+
D
0
C
4
*
p
B
*
lim
c
→
1
+
−
h
2
(
1
)
c
2
−
1
(25)
Using Miller and Guy, assuming that the Laplace transform of
f
(
t
)
is
f
*
(
p
)
,
p
can be given by the following discrete points
p
=
δ
(
β
+
1
+
k
)
(26)
where
δ
>
0
,
β
>
−
1
,
k
=
0
,
1
,
2
,
⋯
, Then it can be concluded that
f
(
t
)
=
∑
T
=
0
N
B
T
P
T
(
0
,
β
)
2
e
−
δ
t
−
1
,
δ
>
0
,
β
>
−
1
(27)
Here,
P
T
(
0
,
β
)
is a Jacobi polynomial, and the coefficients
B
T
are determined by the following system of equations
δ
f
*
(
δ
(
β
+
1
+
k
)
)
=
∑
m
=
0
k
k
(
k
−
1
)
⋯
k
−
(
m
−
1
)
(
k
+
β
+
1
)
(
k
+
β
+
2
)
⋯
(
k
+
β
+
1
+
m
)
B
T
(28)
Based on the definitions of stress and electric displacement, the stress intensity factors
K
I
I
I
and the electric displacement intensity factors
K
I
V
at the crack tips
x
=
a
1
and
x
=
b
1
can be obtained as
K
I
I
I
(
a
1
)
K
I
V
(
a
1
)
=
2
π
(
a
1
−
x
)
x
→
a
1
−
0
σ
y
z
(
x
,
0
,
t
)
D
y
(
x
,
0
,
t
)
(29)
K
I
I
I
(
b
1
)
K
I
V
(
b
1
)
=
2
π
(
x
−
b
1
)
x
→
b
1
+
0
σ
y
z
(
x
,
0
,
t
)
D
y
(
x
,
0
,
t
)
(30)
Here, a non-dimensional function is defined
K
=
K
I
I
I
(
t
,
v
)
/
K
I
I
I
(
0
,
v
)
(31)
Here,
K
I
I
I
(
0
,
v
)
represents the stress intensity factor under static load conditions.
5. Far-Field Forces and Electric Impact Loading
For the crack problem in piezoelectric media, as a cutting-edge branch of multi-physics field-coupled fracture mechanics, its core difficulty lies in accurately describing the complex boundary conditions at the crack tip and its surface. These conditions involve not only mechanical stress and strain but are also closely related to the electric field and electric displacement, forming a highly nonlinear boundary-value problem.
The classical strict electro-mechanical boundary conditions require that the normal displacement components at the surface of the piezoelectric medium (as the defect interface) and within the material must be continuous in the vicinity of the crack tip and at the crack interface, and the potential or electric field may also need to satisfy specific continuity or jump conditions. This strict condition stems from the principle of physical continuity and is currently the most accurate theoretical description. However, as stated in the original text, this greatly limits the scope of obtaining analytical solutions or closed-form solutions, making it extremely difficult to deeply understand and predict the fracture behavior of piezoelectric materials under complex conditions, severely restricting the application of related theories in engineering practice.
To overcome this challenge, the engineering and academic communities have developed a practical and effective approximation method. The core idea is to simplify the electric field based on the “transparency” of the crack surface. Currently, the most widely accepted and applied approach is to classify crack-boundary conditions into two categories: impermeable cracks and permeable cracks. The impermeable crack model, also known as the D-P boundary condition, is the most commonly used simplified model. It is based on a key assumption: the crack surface is completely insulated and does not allow any free charges, so the free charge density on the crack surface is zero. According to Gauss’s law, this means that the normal electric displacement component on the crack surface must be zero. Additionally, this model typically assumes that the potential on the crack surface remains constant. This boundary condition significantly simplifies the mathematical processing as it directly eliminates the influence of the crack surface as an electric field boundary, making the crack surface a “free surface” or an “equipotential surface”.
Although from a physical perspective this might be overly idealized (completely preventing the electric field from penetrating), in many engineering applications, especially when the crack propagation direction is not consistent with the electric field direction or when the crack size is much smaller than the characteristic wavelength, the D-P condition can provide results with reasonable accuracy, becoming a common benchmark for fracture mechanics analysis (such as calculating stress-intensity factors, electric-displacement-intensity factors) and structural-safety assessment. The permeable-crack model takes the opposite extreme, allowing the electric field to pass through the crack surface. In the permeable model, it is assumed that the potentials on the upper and lower surfaces of the crack are equal (i.e., the potential is continuous), or more generally, it allows for a certain jump in the potential on both sides of the crack surface, but the key is that the normal electric displacement components on both sides of the crack surface remain continuous. This means that the crack surface is no longer completely insulated but allows current to pass through. The permeable model typically better reflects the actual distribution of the electric field at the crack, especially in strong electric-field effects or specific geometric structures. However, its mathematical processing is usually more complex, and analytical solutions are more difficult to obtain.
By introducing these approximate boundary conditions, although the accuracy of the model is sacrificed to some extent, it significantly reduces the difficulty of solving the piezoelectric crack problem, enabling the use of analytical methods, semi-analytical methods to analyze fracture behavior, calculate stress intensity factors to be possible, and assess the safety of structures under complex load coupling. These approximate models provide important theoretical tools and engineering methods for understanding and predicting the failure mechanism of piezoelectric intelligent structures (such as piezoelectric sensors, detectors, and actuators).
In this study, we investigate the crack-growth problem subjected to stress and electric displacement impact loading at infinity, denoted by τ∞H(t) and D∞H(t), respectively. The boundary conditions described in Equations (5) and (7) are revised and transformed into the following forms
σ
y
z
(
x
,
d
1
,
t
)
=
σ
y
z
(
x
,
−
d
2
,
t
)
=
0
(32)
w
(
x
,
0
+
,
t
)
=
w
(
x
,
0
−
,
t
)
,
D
y
(
x
,
0
+
,
t
)
=
D
y
(
x
,
0
−
,
t
)
(33)
σ
y
z
(
x
,
0
+
,
t
)
=
σ
y
z
(
x
,
0
−
,
t
)
=
−
τ
∞
H
(
t
)
D
y
(
x
,
0
+
,
t
)
=
D
y
(
x
,
0
−
,
t
)
=
D
*
−
D
∞
H
(
t
)
(34)
where
D
*
is the electrical displacement in the crack plane. Take the Laplace transform of Equation (34), we obtain
σ
y
z
*
(
x
,
0
+
,
p
)
=
σ
y
z
*
(
x
,
0
−
,
p
)
=
−
τ
∞
p
,
D
y
*
(
x
,
0
+
,
p
)
=
D
y
*
(
x
,
0
−
,
p
)
=
−
D
∞
−
D
*
p
(35)
5.1. Impermeable Crack (D-P)
For impermeable cracks, the electric potential in the crack cavity can be reasonably neglected. This is because the dielectric constant of piezoelectric materials is usually three to four orders of magnitude higher than that of air or vacuum. Therefore, the physical quantities appearing in Equation (35)
D
*
meet the following conditions [
23
]:
D
*
=
0
(36)
According to the aforementioned derivation process, the stress-intensity factor represented by
a
1
and
b
1
can be expressed as
lim
x
→
a
1
−
σ
y
z
*
(
x
,
0
,
p
)
=
σ
∞
C
1
*
p
B
*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
+
σ
∞
C
2
*
p
B
*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(37)
lim
x
→
b
1
+
σ
y
z
*
(
x
,
0
,
p
)
=
σ
∞
C
1
*
p
B
*
lim
c
→
1
+
−
h
1
(
1
)
c
2
−
1
+
σ
∞
C
2
*
p
B
*
lim
c
→
1
+
−
h
2
(
1
)
c
2
−
1
(38)
lim
x
→
a
1
−
D
y
*
(
x
,
0
,
p
)
=
D
∞
C
3
*
p
B
*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
+
D
∞
C
4
*
p
B
*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(39)
lim
x
→
b
1
+
D
y
*
(
x
,
0
,
p
)
=
D
∞
C
3
*
p
B
*
lim
c
→
1
+
−
h
1
(
1
)
c
2
−
1
+
D
∞
C
4
*
p
B
*
lim
c
→
1
+
−
h
2
(
1
)
c
2
−
1
(40)
In terms of the definitions of stress and electric displacement, the corresponding factors at the crack fronts can be formally derived as:
K
I
I
I
(
a
1
)
K
I
V
(
a
1
)
=
2
π
(
a
1
−
x
)
x
→
a
1
−
0
σ
y
z
(
x
,
0
,
t
)
D
y
(
x
,
0
,
t
)
(41)
K
I
I
I
(
b
1
)
K
I
V
(
b
1
)
=
2
π
(
x
−
b
1
)
x
→
b
1
+
0
σ
y
z
(
x
,
0
,
t
)
D
y
(
x
,
0
,
t
)
(42)
According to Equations (41) and (42) and combined with the virtual-crack-closure technique, the corresponding energy release rate G can be obtained [
24
].
5.2. Permeable Crack
For the penetrating cracks, when the material in the text degenerates into an isotropic material, the coupling degree of the mechanical and electrical loads can be expressed as
c
55
(
k
)
(
∂
2
w
∂
X
2
+
∂
2
w
∂
Y
2
)
+
e
15
(
k
)
(
∂
2
ϕ
∂
X
2
+
∂
2
ϕ
∂
Y
2
)
=
ρ
∂
2
w
∂
t
2
(43)
e
15
(
k
)
(
∂
2
w
∂
X
2
+
∂
2
w
∂
Y
2
)
−
ε
11
(
k
)
(
∂
2
ϕ
∂
X
2
+
∂
2
ϕ
∂
Y
2
)
=
0
(44)
e
15
D
∞
ε
11
τ
∞
=
λ
(45)
Note: In the previous calculation,
c
55
(
k
)
=
c
44
(
k
)
,
e
15
(
k
)
=
e
24
(
k
)
and
ε
11
(
k
)
=
ε
22
(
k
)
.
For the permeable crack model, the continuity of the electric potential at the upper and lower surfaces of the crack is the core electrical boundary condition, which means that the electric potential of the media on both sides of the crack is equal everywhere on the crack surface
h
ϕ
(
x
)
=
0
(46)
According to the aforementioned derivation process, the stress-intensity factor represented by
a
1
and
b
1
can be expressed as
lim
x
→
a
1
−
σ
y
z
*
(
x
,
0
,
p
)
=
σ
∞
C
1
*
p
B
*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
+
σ
∞
C
2
*
p
B
*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(47)
lim
x
→
b
1
+
σ
y
z
*
(
x
,
0
,
p
)
=
σ
∞
C
1
*
p
B
*
lim
c
→
1
+
−
h
1
(
1
)
c
2
−
1
+
σ
∞
C
2
*
p
B
*
lim
c
→
1
+
−
h
2
(
1
)
c
2
−
1
(48)
lim
x
→
a
1
−
D
y
*
(
x
,
0
,
p
)
=
λ
ε
11
τ
∞
C
3
*
e
15
p
B
*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
+
λ
ε
11
τ
∞
C
4
*
e
15
p
B
*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(49)
lim
x
→
b
1
+
D
y
*
(
x
,
0
,
p
)
=
λ
ε
11
τ
∞
C
3
*
e
15
p
B
*
lim
c
→
1
+
−
h
1
(
1
)
c
2
−
1
+
λ
ε
11
τ
∞
C
4
*
e
15
p
B
*
lim
c
→
1
+
−
h
2
(
1
)
c
2
−
1
(50)
6. Numerical Examples
To verify the above method and analyze its application, we simulated some numerical examples. In the following examples,
a
represents the half-length of the crack.
6.1. The Cracks Change in Accordance with the Influence of Elastic Constants
The variation in the dimensionless function with different elastic constant ratios is shown in
Figure 2
a,b. When the elastic constant ratio
ς
c
=
c
44
(
2
)
/
c
44
(
1
)
was different, the piezoelectric constant ratio and dielectric constant ratio were both 1. The crack changes under different thickness ratios of materials 1 and 2 (
d
1
/
d
2
=
0.0125
and
d
1
/
d
2
=
0.01
) were analyzed respectively in the cases of elastic constant ratio
ς
c
=
2
(
Figure 2
a) and
ς
c
=
1.5
(
Figure 2
b).
As can be seen from the figure, the non-dimensional function decreases continuously as
d
1
/
a
increases and gradually approaches a stable value. Moreover, when the elastic modulus is large, the stable value of the non-dimensional function is also larger. At the same time, it can be observed that when the thickness of material 2 is fixed, the larger the thickness of material 1, the smaller the influence on the non-dimensional function. That is, the smaller the thickness of material 1, the greater the stress intensity generated. Therefore, in application, if we make material 1 relatively thinner, it will have higher practical application value.
6.2. The Variation in Crack Under Impact Time and Crack Length
It can be clearly seen from the two-dimensional and three-dimensional curves presented in
Figure 3
a,b that the dimensionless dynamic stress intensity factor at the crack tip of the piezoelectric bimaterial interface shows obvious transient dynamic behaviors under impact load. Over time, the dimensionless function first rises rapidly, reaching a significant peak within a short period of time, then begins to gradually decline, presenting small oscillations within a certain range, and finally converges gradually to the corresponding stress intensity factor value under the static load.
This change process fully demonstrates the comprehensive influence of stress-wave propagation, reflection, superposition, and energy dissipation under impact load on the mechanical field at the crack tip and also reflects the significant differences between dynamic fracture behavior and quasi-static fracture behavior. The appearance of the peak indicates that in the initial stage of impact, the stress wave focuses strongly at the crack tip, resulting in a significant enhancement of the local mechanical and electric field coupling effect, while the oscillation convergence in the later stage indicates that the system gradually stabilizes and the dynamic effect gradually decays.
At the same time, the numerical results also show that, while keeping the thickness of material 1 unchanged, the crack length has a significant positive regulatory effect on the dimensionless function. As the crack length increases, the dimensionless dynamic stress intensity factor shows a significant increasing trend overall, which means that the expansion of the crack size will further intensify the stress concentration at the interface and increase the risk of dynamic expansion and unstable fracture of the crack.
These laws not only reveal the intrinsic relationship between geometric parameters, load forms and dynamic fracture behavior, but also provide important theoretical basis and data support for the strength design, damage identification and safety assessment of piezoelectric composite material structures under complex dynamic conditions such as impact and vibration.
6.3. Impermeable and Permeable Cracks
The variation in the normalized energy-release rate under far-field stress loading is shown in
Figure 4
. The dynamic nondimensional energy-release rate can be expressed as
G
/
G
r
(
G
r
is the dynamic energy-release rate for a stationary crack).
Figure 4
shows that when the far-field shear stress
τ
∞
increases from zero, the normalized energy release rate
G
/
G
r
always shows an upward trend. However, when
τ
∞
decreases from zero (i.e., in a negative shear stress state),
G
/
G
r
first decreases and then increases. This phenomenon demonstrates that the effects of positive and negative shear stresses on crack propagation have significant asymmetry: positive shear stress can continuously increase the energy release rate, thereby promoting crack propagation, while negative shear stress will reduce the energy release rate within a certain small range, manifesting as inhibition or hindrance of crack development. However, as the absolute value of negative shear stress further increases, the energy release rate will rise again. Thus, stress loading can either promote crack propagation or have interference or blocking effects, the specific manifestation depending on the direction and magnitude of the far-field stress. The above change pattern is completely consistent with the theoretical analysis results established by Pak [
23
], further verifying the rationality of the calculation model and the shear fracture theory presented in this paper.
In order to conduct a more in-depth study on the dynamic crack propagation under the coupling impact of mechanics and electricity,
Figure 5
presents the evolution of the normalized dynamic strength factor under different electromagnetic–mechanical coupling coefficients λ. It is clearly observable that as the parameter ct/(2a) related to time increases, the normalized dynamic strength factor
F
0
first rises rapidly, reaching a significant peak, then gradually decreases and eventually stabilizes at a constant value. This characteristic trend, including the sharp rise, peak response, and final steady-state behavior, is highly consistent with the dynamic fracture response studied by Wang [
25
].
6.4. Response-Surface Statistical Analysis
This study employed Response Surface Methodology (RSM), a statistical technique that integrates experimental design, modeling, and optimization techniques, aiming to systematically explore and quantify the influence of multiple controllable factors on a specific response variable, and ultimately determine the conditions for obtaining the optimal response. In this study, we specifically focused on the influence pattern of three key factors—Factor A (
d
1
/
a
), Factor B (
c
p
t
/
a
), and Factor C (
c
44
(
2
)
/
c
44
(
1
)
)—on the Stress Intensity Factor R1 (SIF) response variable.
First, to efficiently and comprehensively obtain data points for modeling, we adopted the Box–Behnken Design (BBD) experimental design method provided by Design-Expert software V22. BBD is a commonly used response surface design method, particularly suitable for situations where it is necessary to estimate the interaction and quadratic effects between factors.
Next, we used the statistical module of Design-Expert to perform a quadratic polynomial regression fit on the data. To assess the reliability and explanatory power of the established model, we conducted rigorous statistical analyses. The model-fitting goodness was measured by the coefficient of determination R
2
(R-squared) and the adjusted coefficient of determination Adjusted R
2
. The R
2
value represents the proportion of the total variation in the response variable explained by the model; the closer it is to 1, the better the model fit. Adjusted R
2
considers the number of independent variables in the model and adjusts R
2
, providing a more realistic reflection of the model’s actual predictive ability. At the same time, we performed an Analysis of Variance (ANOVA). Analysis of Variance (ANOVA) decomposes the contribution of independent variables to the total variation to assess the impact of each factor on the response variable (Stress Intensity Factor R1). The significance of each factor on the response value is mainly achieved by comparing the variance between groups with that within groups, thereby evaluating the accuracy and reliability of the model. The significance of the model is tested by analyzing the sources of error, and the F-value and
p
-value are combined for joint judgment. The larger the F value and the smaller the
p
value, the more significant the model is. When the
p
value is greater than 0.05, the model is not significant. When the
p
-value is less than 0.05, the model is significant, especially when the
p
-value is less than 0.0001, the model shows extremely high significance. According to the data in
Table 1
, the
p
value is less than 0.0001, so the model is significant. The ANOVA results not only validated the statistical significance of the entire model (through F-test, typically requiring a
p
-value < 0.05) but also tested the significance of individual coefficients (linear, quadratic, interaction terms) (through
t
-test, typically requiring a
p
-value < 0.05 or 0.01), helping us understand which factors and their interactions significantly affect SIF. Residual analysis was also conducted to further validate the model’s predictive ability and the reasonableness of its assumptions.
Finally, based on the established and validated quadratic response surface model, we utilized the optimization tools built into the Design-Expert software to perform optimization calculations on the combination of Factors A, B, and C. The optimization objective was set to maximize the Stress Intensity Factor (SIF) response value R1. The optimization process considered all significant factor effects (including linear, quadratic, and interaction effects) and may have set operational constraints on the factors (such as factors not exceeding their physical or process-allowable ranges). The optimization algorithm in Design-Expert searched the factor space to find the optimal combination of Factors A, B, and C that maximizes SIF under the constraint conditions. This optimal combination not only provides clear guidance for practical applications but also intuitively demonstrates the core value of RSM—namely, finding optimal operating conditions from complex factor interactions through mathematical modeling and optimization.
The results of the analysis of variance (as shown in
Table 1
and
Table 2
) show that the F value of the regression model is 706.4, and the
p
values of each model are all less than 0.05, indicating the significance of the model. Among them, the
p
values of A-A, B-B, A
2
, B
2
, and C
2
are all less than 0.0001, indicating that the model is extremely significant and highly reliable, and the experimental design is reasonable. Therefore, the results of the analysis of variance indicate that the fitting model is accurate and statistically significant.
According to the data, a second-order response surface mathematical model was obtained through analysis. The second-order response surface regression model regarding the three variables A, B, and C, as well as R1, is shown in Equation (51).
R
1
=
−
12.2865
+
0.8313
A
+
1.5650
B
+
14.0871
C
+
0.1817
AB
+
0.1818
AC
−
0.0691
BC
−
0.8666
A
2
−
0.4544
B
2
−
4.07760
C
2
(51)
Figure 6
depicts the relationship between the predicted values and the actual values. The horizontal axis represents the actual values, and the vertical axis represents the predicted values. Under ideal conditions, if the predicted value is equal to the actual value, it lies on the black straight line. Different data points are represented by squares of different colors. By observing the image, it can be found that the vast majority of data points are concentrated near the black straight line, which directly indicates that the predicted values are close to the actual values. This also shows that the established prediction model has good accuracy and reliability, can predict the relevant values relatively accurately, and can effectively reflect the true relationship between the variables.
It can be seen from
Table 2
that the interaction terms AC and BC are statistically insignificant. This does not mean that there are no interaction effects at the physical level, but reflects that the strengths of these interactions are weaker compared with the main effects of factors A, B, and C. During the experiments, the coupled electromechanical effects corresponding to AC and BC do exist objectively, but their influences on the response variables (e.g., interfacial electromechanical coupling strength, stress–electric coupling distribution at the crack tip, etc.) are overshadowed by the dominant main effects. For example, A (thickness of Material 1/crack length) directly determines the constraint effect of Material 1 on crack propagation, and C (ratio of elastic constants of Material 2 to Material 1) mainly affects the mechanical matching across the interface. Their independent effects on the electromechanical coupling are more significant, whereas their synergistic regulation only induces marginal changes. In addition, we have verified the residual distribution of the model and found no systematic deviation, indicating that the insignificance of AC and BC does not result from improper model specification, but truly reflects the weak intensity of these interactions.
From the perspective of physical mechanism, the coupled electromechanical effects at the interface are essentially a complex process involving multi-factor synergy. According to the variable definitions, the synergistic pathways of AC and BC are inherently constrained by the geometrical parameters and mechanical properties of the materials, limiting their interactive contributions to a small range. Specifically, for the interaction term BC, the variation in crack propagation length is directly related to the stress concentration at the crack tip, while the elastic constant ratio determines the stress-transfer efficiency across the interface. During crack propagation, the regulation of stress transfer by the elastic constant ratio is overshadowed by the stress-concentration effect at the crack tip, producing only weak local effects that cannot be manifested as statistically significant changes in the macroscopic response variables. For the interaction term AC, their synergistic effect theoretically affects the distribution of interfacial electromechanical coupling. However, when the thickness of Material 1 is sufficient to form a stable load-bearing structure, the contribution of the interaction to the response variable becomes weak. Due to its small amplitude, it does not reach the threshold of statistical significance in the macroscopic response model, thus showing statistical insignificance.
In summary, the statistical insignificance of the interaction terms AC and BC is a reasonable reflection of the coupled electromechanical effects at the interface under the constraints of specific geometrical parameters and mechanical properties, rather than a contradiction with physical reality. The key to reconciling them is to clarify that “statistical significance” depends on the detectability of effects within the experimental design and modeling framework, while “physical existence” refers to the inherent synergistic behavior of the factors. These two concepts are complementary rather than mutually exclusive.
Through response-surface analysis of the parameters, A 3D response-surface map (as shown in
Figure 7
) and a contour map (as shown in
Figure 7
) were obtained. As shown in
Figure 7
, in the three-dimensional response-surface plot, the horizontal axis represents parameters A and B respectively, and the vertical axis represents the corresponding response values. The shape and gradient changes in the response surface intuitively reflect the variation pattern of the response values under the combined effect of the two factors: the red area corresponds to higher response values, the blue area corresponds to lower response values, and the degree of color gradient clearly shows the increase or decrease range of the response values.
In the contour plot, the same pattern is also presented, that is, the color that is closer to red indicates a higher response value, and the color that is closer to blue indicates a lower response value. Overall, the response-surface plot and the contour plot can intuitively, comprehensively and quantitatively reveal the spatial distribution characteristics, interaction influence rules and extreme distribution areas of the response values with respect to parameters A and B, and can provide reliable visual evidence and theoretical support for the subsequent identification of key parameters, analysis of the influencing mechanism and selection of the optimal parameter range.
In the three-dimensional response-surface plot, the horizontal axis is A and C, and the vertical axis R1 represents the obtained response values (as shown in
Figure 8
). This surface is used to reflect the continuous pattern of response values changing with the two independent variables. The steeper the surface, the stronger the influence of the corresponding independent variable on the response value in that area. If the surface is flat, the influence of the variable on the response value is relatively weak. For the contour plot, the color represents the gradual change (green-yellow-red), indicating high or low response values. The red area has the highest response value, while the green area has the lowest response value. Points on the same curve have the same response value. The denser the curve, the more intense the change in the response value in that area. Within the 1.5 range, the contours are elliptical, indicating a significant interaction between the two independent variables and the change in response values.
Figure 9
shows the three-dimensional response-surface plot and contour plot corresponding to the joint variation in independent variables B and C. From the figure, it can be seen that the response value shows a continuous and smooth trend as variables B and C change. The color in the figure gradually transitions from blue to red, intuitively reflecting the distribution pattern of the response value from low to high. The response surface as a whole presents a morphological characteristic of rapid increase first and then tending to be stable. This indicates that when B and C increase to a certain range, the growth rate of the response value gradually slows down and tends to stabilize.
Further analysis reveals that when the value of variable B is within the range of 1.5 to 2, the response value R1 can reach the optimal range of 1.8 to 2. At the same time, from the contour plot, a significant stretched elliptical distribution feature can be clearly observed. This feature directly indicates that there is a significant interaction between independent variables B and C, and their influence on the response value is not independent of each other but shows a clear coupling effect.
Depicted in
Figure 10
are the response values of A, B, and C under the maximum stress intensity factor condition. From the figure, it can be seen that the optimal values of A, B and C are 0.4045, 1.6797 and 1.9035 respectively, and the stress intensity factor reaches its maximum value of 1.3375 at these values. Additionally, it is found that the interaction between B and C is the most significant. The combination of B and C has a more significant impact on the response value than the combinations of other factors.
Figure 10
shows the parameter combination when the stress-intensity factor reaches its global maximum: the corresponding values of variables A, B, and C are 0.4045, 1.6797, and 1.9035, respectively, at which the peak stress intensity factor is 1.3375. Meanwhile, the response surface analysis reveals that the interaction effect between factors B and C is statistically significant (
p
< 0.05), and their coupled influence on the stress intensity factor is significantly stronger than that of other factor combinations (see
Table 2
for the significance analysis of interaction terms).
It should be clearly defined that the term “optimal parameters” here refers only to the mathematical extreme solution from the response-surface model, specifically the parameter point that maximizes the stress-intensity factor, rather than the optimal design scheme based on structural safety criteria. As a key parameter characterizing the stress-field intensity at the crack tip, the magnitude of the stress-intensity factor directly determines the crack initiation threshold and propagation rate. A higher stress-intensity factor indicates a more severe stress concentration at the crack tip, a greater probability of brittle fracture or fatigue failure, and lower structural safety. Therefore, the parameter combination shown in
Figure 9
corresponds to the most unfavorable service condition of the structure, whose core value is to provide a boundary reference for dangerous working conditions in engineering design, rather than a recommended safe parameter range.
In practical structural safety design and fracture control, the high-risk parameter combination shown in
Figure 9
should be avoided. We have supplemented the distinction between mathematically optimal extreme values and engineering safety-optimized solutions in the revised manuscript to prevent confusion of academic concepts, ensure the consistency between research conclusions and engineering practice, and thus provide a more valuable reference for structural safety design.
7. Conclusions
This paper conducts a systematic study on the dynamic propagation of orthogonal anisotropic interface cracks in piezoelectric bimaterial structures. By introducing Laplace transform and Fourier transform to convert the governing equations and combining the Chebyshev collocation method to discretize the original problem into algebraic equations, the dynamic variation laws of the stress intensity factor at the crack tip and the electric displacement intensity factor are obtained after numerical inverse transformation. Based on this, the influence mechanism of key parameters on the fracture behavior is analyzed by using the response surface method, and the following main conclusions are obtained:
(1)
The elastic parameters have a significant regulatory effect on the dynamic propagation of interface cracks. Within a certain range, the larger the elastic parameters, the higher the dimensionless dynamic stress-intensity factor at the crack tip, indicating that the elastic properties of the material directly affect the degree of stress concentration and the driving ability of the crack.
(2)
The impact-load application time significantly changes the dynamic response characteristics of the crack. With the passage of time, the dimensionless stress-intensity factor shows a pattern of first increasing, reaching a peak, and then oscillating and converging, eventually approaching the corresponding value under static load, reflecting the combined effect of stress wave propagation, reflection, and energy dissipation.
(3)
The crack length is an important geometric factor affecting the fracture risk. The larger the crack length, the more significant the stress concentration at the tip, and the corresponding dimensionless dynamic stress-intensity factor increases. The structural fracture failure risk is higher.
(4)
For non-permeable cracks, the dual effect of forward shear stress promoting crack propagation and negative shear stress first inhibiting and then driving the crack propagation has verified the rationality of the shear-fracture calculation model established in this paper. For the dimensionless dynamic strength factors corresponding to different force–electric coupling coefficients λ, they all show a typical characteristic of rapid increase—reaching a peak—slow decrease—approaching stability as
c
p
t
/
(
2
a
)
changes. This further proves that the model in this paper can accurately describe the transient fracture-mechanics behavior of piezoelectric materials under coupled loads.
(5)
The constructed prediction model has high accuracy and strong reliability, with a determination coefficient R
2
of 0.9886. The real values and predicted values are highly consistent, providing an efficient and accurate data-driven method for the rapid prediction of the stress-intensity factor at the crack tip.
(6)
The key parameter combinations that maximize the stress-intensity factor were obtained through response surface analysis: A = 0.4045, B = 1.6797, C = 1.9035, corresponding to a peak stress intensity factor of 1.3375. Significance analysis indicates that there is a significant interaction between factors B and C (
p
< 0.05), and their coupling effect is stronger than other parameter combinations.
(7)
The “optimal parameter combination” obtained in this paper is only a mathematical extremum point. Even the parameter combination that maximizes the stress-intensity factor corresponds to the most dangerous service condition of the structure, rather than an engineering safety optimization scheme. This result can be used to identify high-risk parameter intervals, providing important theoretical references for the fracture prevention, safety design, and boundary determination of dangerous conditions of engineering structures.
Overall, the analytical–numerical combined analysis framework established in this paper not only enriches the dynamic fracture theory of piezoelectric bimaterial interface cracks but also provides feasible methods for the safety assessment, parameter optimization, and failure prevention of intelligent materials and structures. It has certain theoretical value and practical significance for promoting the reliable application of piezoelectric structures in engineering fields. In future research, we will incorporate the thickness ratio (d
1
/d
2
) as a key parameter into our analysis and conduct systematic quantitative studies using appropriate numerical simulation software. This extension aims to more comprehensively reveal the influence law of the thickness ratio on the dynamic fracture behavior of piezoelectric bimaterial interface cracks, further improving the completeness and engineering applicability of the proposed theoretical framework.
Author Contributions
Conceptualization, Y.Z. and J.L.; methodology, J.L.; software, J.L.; validation, Y.Z., J.L. and X.L.; formal analysis, X.L.; investigation, J.M.; resources, Y.Z. and J.M.; data curation, Y.Z.; writing—original draft preparation, Y.Z.; writing—review and editing, X.L.; visualization, J.M.; supervision, J.L.; project administration, Y.Z.; funding acquisition, J.M. and J.L. All authors have read and agreed to the published version of the manuscript.
Funding
This research was funded by National Natural Science Foundation of China grant number 12301591, and also by the Shanxi Provincial Key Research and Development Project grant number 202102090301027. The APC was funded by the Taiyuan University of Science and Technology doctoral research start-up fund grant number 20252008.
Data Availability Statement
The original contributions presented in this study are included in the article, further inquiries can be directed to the corresponding author.
Conflicts of Interest
The authors declare no conflicts of interest.
Appendix A
ϖ
1
(
k
)
=
(
−
ε
22
(
k
)
c
44
(
k
)
(
(
c
55
(
k
)
c
44
(
k
)
−
v
2
c
k
¯
2
)
s
2
+
2
v
p
s
i
c
k
¯
2
+
p
2
c
k
¯
2
)
−
2
e
24
(
k
)
e
15
(
k
)
s
2
c
44
(
k
)
c
44
(
k
)
−
ε
11
(
k
)
s
2
c
44
(
k
)
)
/
(
(
ε
22
(
k
)
c
44
(
k
)
+
(
e
24
(
k
)
c
44
(
k
)
)
2
)
,
ϖ
2
(
k
)
=
(
ε
11
(
k
)
c
44
(
k
)
(
(
c
55
(
k
)
c
44
(
k
)
−
v
2
c
k
¯
2
)
s
4
+
2
v
p
s
3
i
c
k
¯
2
+
p
2
s
2
c
k
¯
2
)
+
e
15
(
k
)
e
15
(
k
)
s
4
c
44
(
k
)
c
44
(
k
)
)
/
(
(
ε
22
(
k
)
c
44
(
k
)
+
(
e
24
(
k
)
c
44
(
k
)
)
2
)
,
ϑ
1
(
k
)
=
−
72
ϖ
1
(
k
)
ϖ
2
(
k
)
−
2
ϖ
1
(
k
)
3
54
+
(
72
ϖ
1
(
k
)
ϖ
2
(
k
)
−
2
ϖ
1
(
k
)
3
)
2
2916
−
(
12
ϖ
2
(
k
)
+
ϖ
1
(
k
)
2
)
3
729
3
,
ϑ
1
(
k
)
=
−
72
ϖ
1
(
k
)
ϖ
2
(
k
)
−
2
ϖ
1
(
k
)
3
54
−
(
72
ϖ
1
(
k
)
ϖ
2
(
k
)
−
2
ϖ
1
(
k
)
3
54
)
2
−
(
12
ϖ
2
(
k
)
+
ϖ
1
(
k
)
2
9
)
3
3
,
γ
1
(
k
)
=
−
−
2
ϖ
1
(
k
)
/
3
+
ϑ
1
(
k
)
+
ϑ
2
(
k
)
+
−
4
ϖ
1
(
k
)
/
3
−
ϑ
1
(
k
)
−
ϑ
2
(
k
)
2
,
γ
2
(
k
)
=
−
−
2
ϖ
1
(
k
)
/
3
+
ϑ
1
(
k
)
+
ϑ
2
(
k
)
−
−
4
ϖ
1
(
k
)
/
3
−
ϑ
1
(
k
)
−
ϑ
2
(
k
)
2
,
γ
3
(
k
)
=
−
2
ϖ
1
(
k
)
/
3
+
ϑ
1
(
k
)
+
ϑ
2
(
k
)
+
−
4
ϖ
1
(
k
)
/
3
−
ϑ
1
(
k
)
−
ϑ
2
(
k
)
2
,
γ
4
(
k
)
=
−
2
ϖ
1
(
k
)
/
3
+
ϑ
1
(
k
)
+
ϑ
2
(
k
)
−
−
4
ϖ
1
(
k
)
/
3
−
ϑ
1
(
k
)
−
ϑ
2
(
k
)
2
,
A
*
=
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
(
−
χ
11
(
1
)
χ
22
(
2
)
+
χ
12
(
2
)
χ
21
(
1
)
+
χ
12
(
1
)
χ
22
(
2
)
−
χ
12
(
2
)
χ
22
(
1
)
)
+
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
(
−
χ
11
(
2
)
χ
21
(
1
)
+
χ
11
(
1
)
χ
21
(
2
)
−
χ
12
(
1
)
χ
21
(
2
)
+
χ
11
(
2
)
χ
22
(
1
)
)
,
A
1
=
−
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
χ
22
(
2
)
+
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
χ
21
(
2
)
,
B
1
=
−
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
χ
11
(
2
)
+
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
χ
12
(
2
)
,
A
2
=
−
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
χ
21
(
2
)
+
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
χ
22
(
2
)
,
B
2
=
−
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
χ
12
(
2
)
+
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
χ
11
(
2
)
,
A
1
*
=
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
(
χ
21
(
1
)
−
χ
22
(
1
)
)
,
B
1
*
=
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
(
χ
12
(
1
)
−
χ
11
(
1
)
)
,
A
2
*
=
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
(
χ
22
(
1
)
−
χ
21
(
1
)
)
,
B
2
*
=
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
(
χ
11
(
1
)
−
χ
12
(
1
)
)
,
C
1
*
=
(
(
c
44
(
1
)
χ
11
(
1
)
+
e
24
(
1
)
χ
21
(
1
)
)
γ
1
(
1
)
−
(
c
44
(
1
)
χ
12
(
1
)
+
e
24
(
1
)
χ
22
(
1
)
)
γ
2
(
1
)
)
(
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
χ
21
(
2
)
−
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
χ
22
(
2
)
)
,
C
2
*
=
(
(
c
44
(
1
)
χ
11
(
1
)
+
e
24
(
1
)
χ
21
(
1
)
)
γ
1
(
1
)
−
(
c
44
(
1
)
χ
12
(
1
)
+
e
24
(
1
)
χ
22
(
1
)
)
γ
2
(
1
)
)
(
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
χ
12
(
2
)
−
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
χ
11
(
2
)
)
,
C
3
*
=
(
(
e
15
(
1
)
χ
11
(
1
)
−
ε
22
(
1
)
χ
21
(
1
)
)
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
−
(
e
15
(
1
)
χ
12
(
1
)
−
ε
22
(
1
)
χ
22
(
1
)
)
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
)
γ
1
(
1
)
γ
2
(
2
)
(
χ
21
(
1
)
−
χ
22
(
1
)
)
,
C
4
*
=
(
(
e
15
(
1
)
χ
11
(
1
)
−
ε
22
(
1
)
χ
21
(
1
)
)
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
−
(
e
15
(
1
)
χ
12
(
1
)
−
ε
22
(
1
)
χ
22
(
1
)
)
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
)
γ
1
(
1
)
γ
2
(
2
)
(
χ
12
(
1
)
−
χ
11
(
1
)
)
,
B
*
=
(
c
44
(
2
)
χ
11
(
2
)
+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
(
−
χ
11
(
1
)
χ
22
(
2
)
+
χ
12
(
2
)
χ
21
(
1
)
+
χ
12
(
1
)
χ
22
(
2
)
−
χ
12
(
2
)
χ
22
(
1
)
)
+
(
c
44
(
2
)
χ
12
(
2
)
+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
(
−
χ
11
(
2
)
χ
21
(
1
)
+
χ
11
(
1
)
χ
21
(
2
)
−
χ
12
(
1
)
χ
21
(
2
)
+
χ
11
(
2
)
χ
22
(
1
)
)
.
References
Singh, V.; Jangid, K.; Bui, T.Q. A study of strain and electric field gradient effects on two collinear cracks in an arbitrary polarized piezoelectric layer.
Eur. J. Mech.-A/Solids
2025
,
114
, 105723. [
Google Scholar
] [
CrossRef
]
Mu, X.; Zhu, Z.; Zhang, L.; Gao, Y. The Green’s functions of two-dimensional piezoelectric quasicrystal semi-infinite crack and finite crack.
Appl. Math. Model.
2025
,
137
, 115732. [
Google Scholar
] [
CrossRef
]
Li, Y.; Yan, S.; Zhao, M.; Ren, J. Fracture Analysis of Planar Cracks in 3D Thermal Piezoelectric Semiconductors.
Int. J. Mech. Sci.
2024
,
273
, 109212. [
Google Scholar
] [
CrossRef
]
Tan, Y.; Rao, W.; Wan, K.; Peng, K.; Zhao, J.; Li, X. Phase-field model for fatigue crack growth in piezoelectrics: Energetically consistent boundary condition.
Int. J. Solids Struct.
2025
,
316
, 113378. [
Google Scholar
] [
CrossRef
]
Nourazar, M.; Yang, W.; Chen, Z. Fracture analysis of a curved crack in a piezoelectric plane under general thermal loading.
Eng. Fract. Mech.
2023
,
284
, 109208. [
Google Scholar
] [
CrossRef
]
Cao, T.; Zhang, Y.; Qin, T. Research on interaction of two straight cracks embedded in different spaces of piezoelectric biomaterial.
Eng. Anal. Bound. Elem.
2025
,
179
, 106401. [
Google Scholar
] [
CrossRef
]
Yu, H.; Zhu, S.; Ma, H.; Wang, J. Interface crack analysis of piezoelectric laminates considering initial strain.
Int. J. Mech. Sci.
2024
,
271
, 109104. [
Google Scholar
] [
CrossRef
]
Zhu, S.; Yu, H.; Guo, L. Analysis of an interfacial crack between two nonhomogeneous piezoelectric materials using a new domain-independent interaction integral.
Compos. Struct.
2024
,
331
, 117873. [
Google Scholar
] [
CrossRef
]
Zhang, Y.; Li, J.; Xie, X. SH-wave scattering by the interface crack of piezoelectric ceramic polymer composites.
J. Eng. Res.
2024
,
12
, 259–267. [
Google Scholar
] [
CrossRef
]
Zhao, S.; Govorukha, V.; Sheveleva, A.; Loboda, V. On energy release rate for an electrically permeable interface crack between two different 1D hexagonal piezoelectric QCs.
Procedia Struct. Integr.
2023
,
51
, 69–75. [
Google Scholar
] [
CrossRef
]
Loboda, V.; Sheveleva, A.; Komarov, O.; Chapelle, F.; Lapusta, Y. Arbitrary number of electrically permeable cracks on the interface between two one-dimensional piezoelectric quasicrystals with piezoelectric effect.
Eng. Fract. Mech.
2022
,
276
, 108878. [
Google Scholar
] [
CrossRef
]
Sheveleva, A.; Loboda, V.; Lapusta, Y. A conductive crack and a remote electrode at the interface between two piezoelectric materials.
Appl. Math. Model.
2020
,
87
, 287–299. [
Google Scholar
] [
CrossRef
]
Sun, W.; Qu, W.; Gu, Y. Arbitrary-order GFDM for crack analysis of homogeneous and bi-material piezoelectric problems.
Theor. Appl. Fract. Mech.
2025
,
139
, 105107. [
Google Scholar
] [
CrossRef
]
Singh, R. Inspection of dynamic fracture behavior of multiple interfacial cracks emanating from circular holes in functionally graded piezoelectric bi-materials.
Appl. Math. Model.
2025
,
143
, 116041. [
Google Scholar
] [
CrossRef
]
Yan, Z.; Wen, C.; Feng, W.; Zhang, C. Dynamic growth properties of an internal crack in piezoelectric plates based on the improved mechanical energy release rate criterion.
Theor. Appl. Fract. Mech.
2025
,
140
, 105131. [
Google Scholar
] [
CrossRef
]
Zhu, S.; Yu, H.; Wang, Z. Interfacial dynamic impermeable crack analysis in dissimilar piezoelectric materials by a new interaction integral.
Compos. Struct.
2025
,
352
, 118668. [
Google Scholar
] [
CrossRef
]
Liu, H.; Zhu, S. Dynamic analysis of two collinear permeable Mode-I cracks in piezoelectric materials based on non-local piezoelectricity theory.
Multidiscip. Model. Mater. Struct.
2019
,
15
, 1274–1293. [
Google Scholar
] [
CrossRef
]
Liu, H.-T. Dynamic non-local theory solution to a permeable mode-I crack in a piezoelectric medium.
Eng. Fract. Mech.
2017
,
179
, 43–58. [
Google Scholar
] [
CrossRef
]
Liu, H.-T.; Wu, J.-G.; Li, T.-J. Dynamic analytical solution of a limited-permeable mode-I crack in piezoelectric materials based on the non-local theory.
Wave Motion
2019
,
90
, 82–98. [
Google Scholar
] [
CrossRef
]
Afshar, H.; Bagheri, R. Several embedded cracks in a functionally graded piezoelectric strip under dynamic loading.
Comput. Math. Appl.
2018
,
76
, 534–550. [
Google Scholar
] [
CrossRef
]
Wünsche, M.; Sladek, J.; Sladek, V.; Zhang, C.; García-Sánchez, F.; Sáez, A. Dynamic crack analysis in piezoelectric solids under time-harmonic loadings with a symmetric Galerkin boundary element method.
Eng. Anal. Bound. Elem.
2017
,
84
, 141–153. [
Google Scholar
] [
CrossRef
]
Yu, T.; Bui, T.Q.; Liu, P.; Zhang, C.; Hirose, S. Interfacial dynamic impermeable cracks analysis in dissimilar piezoelectric materials under coupled electromechanical loading with the extended finite element method.
Int. J. Solids Struct.
2015
,
67–68
, 205–218. [
Google Scholar
] [
CrossRef
]
Pak, Y.E. Crack extension force in a piezoelectric material.
J. Appl. Mech.
1990
,
57
, 647–653. [
Google Scholar
] [
CrossRef
]
Wang, X.-Y.; Yu, S.-W. Transient response of a crack in piezoelectric strip subjected to the mechanical and electrical impacts:mode-III problem.
Int. J. Solids Struct.
2000
,
37
, 5795–5808. [
Google Scholar
] [
CrossRef
]
Wang, B.-L.; Han, J.-C.; Du, S.-Y. Dynamic response for non-homogeneous piezoelectric medium with multiple cracks.
Eng. Fract. Mech.
1998
,
61
, 607–617. [
Google Scholar
] [
CrossRef
]
Figure 1.
Structure of the piezoelectric bimaterials model with cracks (The green part in the picture indicates the direction of the crack movement).
Figure 1.
Structure of the piezoelectric bimaterials model with cracks (The green part in the picture indicates the direction of the crack movement).
Figure 2.
The variation of
K
with respect to
d
1
/
a
((
a
):
ς
c
=
2
; (
b
):
ς
c
=
1.5
).
Figure 2.
The variation of
K
with respect to
d
1
/
a
((
a
):
ς
c
=
2
; (
b
):
ς
c
=
1.5
).
Figure 3.
(
a
) The 2D plot that
K
varies with
c
p
t
/
a
; (
b
) The three-dimensional graph that
K
varies with
c
p
t
/
a
.
Figure 3.
(
a
) The 2D plot that
K
varies with
c
p
t
/
a
; (
b
) The three-dimensional graph that
K
varies with
c
p
t
/
a
.
Figure 4.
The variation in the normalized energy release rate under far-field stress loading.
Figure 4.
The variation in the normalized energy release rate under far-field stress loading.
Figure 5.
The variation in the normalized intensity factor under different
λ
values.
Figure 5.
The variation in the normalized intensity factor under different
λ
values.
Figure 6.
Comparison chart of real values and predicted values.
Figure 6.
Comparison chart of real values and predicted values.
Figure 7.
The three-dimensional graph and contour map of the interaction between A and B in terms of the influence on the stress intensity factor.
Figure 7.
The three-dimensional graph and contour map of the interaction between A and B in terms of the influence on the stress intensity factor.
Figure 8.
The three-dimensional graph and contour map of the interaction between A and C in terms of the influence on the stress-intensity factor.
Figure 8.
The three-dimensional graph and contour map of the interaction between A and C in terms of the influence on the stress-intensity factor.
Figure 9.
The three-dimensional graph and contour map of the interaction between B and C in terms of the influence on the stress intensity factor.
Figure 9.
The three-dimensional graph and contour map of the interaction between B and C in terms of the influence on the stress intensity factor.
Figure 10.
The response values of A, B and C when the stress-intensity factor reaches its maximum value (The red and grey dots represent the corresponding data values).
Figure 10.
The response values of A, B and C when the stress-intensity factor reaches its maximum value (The red and grey dots represent the corresponding data values).
Table 1.
Variance analysis of predictive models.
Table 1.
Variance analysis of predictive models.
Standard Error
Correlation Coefficient
Correction Coefficient
F
p
Significant
0.0285
0.9975
0.9886
706.4
<0.0001
significant
Table 2.
Variance analysis of predictive model of response value.
Table 2.
Variance analysis of predictive model of response value.
Source
Sum of Square
Df
Mean Square
F-Value
p
-Value
Significant
A-A
0.1013
1
0.1013
125.03
<0.0001
significant
B-B
3.42
1
3.42
4214.47
<0.0001
significant
C-C
0.0092
1
0.0092
11.32
0.0120
AB
0.0399
1
0.0399
49.29
0.0002
AC
0.0025
1
0.0025
3.09
0.1224
BC
0.0012
1
0.0012
1.47
0.2642
A
2
0.2894
1
0.2894
357.08
<0.0001
significant
B
2
0.8693
1
0.8693
1072.74
<0.0001
significant
C
2
0.2735
1
0.2735
337.47
<0.0001
significant
Residual
0.0057
7
0.0008
Lack of Fit
0.0034
3
0.0011
2.06
0.2480
not significant
Pure Error
0.0022
4
0.0006
Cor Total
5.16
16
Disclaimer/Publisher’s Note:
The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
© 2026 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the
Creative Commons Attribution (CC BY) license
. |
| Markdown | Next Article in Journal
[Effect of Heat Input on the Hydrogen Embrittlement Sensitivity of CGHAZ of X60 Pipeline Steel](https://www.mdpi.com/1996-1944/19/5/961)
Previous Article in Journal
[A Numerical Study on Optimization of Shape and Dimensions for Cold-Extruded Blank of Copper Pin-Type Heat Dissipation Substrates](https://www.mdpi.com/1996-1944/19/5/962)
## Journals
[Active Journals](https://www.mdpi.com/about/journals) [Find a Journal](https://www.mdpi.com/about/journalfinder) [Journal Proposal](https://www.mdpi.com/about/journals/proposal) [Proceedings Series](https://www.mdpi.com/about/proceedings)
[Topics](https://www.mdpi.com/topics)
## Information
[For Authors](https://www.mdpi.com/authors) [For Reviewers](https://www.mdpi.com/reviewers) [For Editors](https://www.mdpi.com/editors) [For Librarians](https://www.mdpi.com/librarians) [For Publishers](https://www.mdpi.com/publishing_services) [For Societies](https://www.mdpi.com/societies) [For Conference Organizers](https://www.mdpi.com/conference_organizers)
[Open Access Policy](https://www.mdpi.com/openaccess) [Institutional Open Access Program](https://www.mdpi.com/ioap) [Special Issues Guidelines](https://www.mdpi.com/special_issues_guidelines) [Editorial Process](https://www.mdpi.com/editorial_process) [Research and Publication Ethics](https://www.mdpi.com/ethics) [Article Processing Charges](https://www.mdpi.com/apc) [Awards](https://www.mdpi.com/awards) [Testimonials](https://www.mdpi.com/testimonials)
[Author Services](https://www.mdpi.com/authors/english)
## Initiatives
[Sciforum](https://sciforum.net/) [MDPI Books](https://www.mdpi.com/books) [Preprints.org](https://www.preprints.org/) [Scilit](https://www.scilit.com/) [SciProfiles](https://sciprofiles.com/) [Encyclopedia](https://encyclopedia.pub/) [JAMS](https://jams.pub/) [Proceedings Series](https://www.mdpi.com/about/proceedings)
## About
[Overview](https://www.mdpi.com/about) [Contact](https://www.mdpi.com/about/contact) [Careers](https://careers.mdpi.com/) [News](https://www.mdpi.com/about/announcements) [Press](https://www.mdpi.com/about/press) [Blog](http://blog.mdpi.com/)
[Sign In / Sign Up](https://www.mdpi.com/user/login)
## Notice
[*clear*]()
## Notice
You are accessing a machine-readable page. In order to be human-readable, please install an RSS reader.
[Continue]() [Cancel]()
[*clear*]()
All articles published by MDPI are made immediately available worldwide under an open access license. No special permission is required to reuse all or part of the article published by MDPI, including figures and tables. For articles published under an open access Creative Common CC BY license, any part of the article may be reused without permission provided that the original article is clearly cited. For more information, please refer to <https://www.mdpi.com/openaccess>.
Feature papers represent the most advanced research with significant potential for high impact in the field. A Feature Paper should be a substantial original Article that involves several techniques or approaches, provides an outlook for future research directions and describes possible research applications.
Feature papers are submitted upon individual invitation or recommendation by the scientific editors and must receive positive feedback from the reviewers.
Editor’s Choice articles are based on recommendations by the scientific editors of MDPI journals from around the world. Editors select a small number of articles recently published in the journal that they believe will be particularly interesting to readers, or important in the respective research area. The aim is to provide a snapshot of some of the most exciting work published in the various research areas of the journal.
Original Submission Date Received: .
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
[](https://www.mdpi.com/)
[*clear*](https://www.mdpi.com/1996-1944/19/5/964) [*zoom\_out\_map*](https://www.mdpi.com/toggle_desktop_layout_cookie "Toggle desktop layout") [*search*](https://www.mdpi.com/1996-1944/19/5/964) [*menu*](https://www.mdpi.com/1996-1944/19/5/964 "MDPI main page")
[](https://www.mdpi.com/)
- [Journals](https://www.mdpi.com/about/journals)
- - [Active Journals](https://www.mdpi.com/about/journals)
- [Find a Journal](https://www.mdpi.com/about/journalfinder)
- [Journal Proposal](https://www.mdpi.com/about/journals/proposal)
- [Proceedings Series](https://www.mdpi.com/about/proceedings)
- [Topics](https://www.mdpi.com/topics)
- [Information](https://www.mdpi.com/authors)
- - [For Authors](https://www.mdpi.com/authors)
- [For Reviewers](https://www.mdpi.com/reviewers)
- [For Editors](https://www.mdpi.com/editors)
- [For Librarians](https://www.mdpi.com/librarians)
- [For Publishers](https://www.mdpi.com/publishing_services)
- [For Societies](https://www.mdpi.com/societies)
- [For Conference Organizers](https://www.mdpi.com/conference_organizers)
- [Open Access Policy](https://www.mdpi.com/openaccess)
- [Institutional Open Access Program](https://www.mdpi.com/ioap)
- [Special Issues Guidelines](https://www.mdpi.com/special_issues_guidelines)
- [Editorial Process](https://www.mdpi.com/editorial_process)
- [Research and Publication Ethics](https://www.mdpi.com/ethics)
- [Article Processing Charges](https://www.mdpi.com/apc)
- [Awards](https://www.mdpi.com/awards)
- [Testimonials](https://www.mdpi.com/testimonials)
- [Author Services](https://www.mdpi.com/authors/english)
- [Initiatives](https://www.mdpi.com/about/initiatives)
- - [Sciforum](https://sciforum.net/)
- [MDPI Books](https://www.mdpi.com/books)
- [Preprints.org](https://www.preprints.org/)
- [Scilit](https://www.scilit.com/)
- [SciProfiles](https://sciprofiles.com/)
- [Encyclopedia](https://encyclopedia.pub/)
- [JAMS](https://jams.pub/)
- [Proceedings Series](https://www.mdpi.com/about/proceedings)
- [About](https://www.mdpi.com/about)
- - [Overview](https://www.mdpi.com/about)
- [Contact](https://www.mdpi.com/about/contact)
- [Careers](https://careers.mdpi.com/)
- [News](https://www.mdpi.com/about/announcements)
- [Press](https://www.mdpi.com/about/press)
- [Blog](http://blog.mdpi.com/)
[Sign In / Sign Up](https://www.mdpi.com/user/login) [Submit](https://www.mdpi.com/1996-1944/19/5/%20%20%20%20%20%20%20%20%20%20%20%20https://susy.mdpi.com/user/manuscripts/upload?journal=materials%0A%20%20%20%20)
*error\_outline* You can access [the new MDPI.com website here](https://www.mdpi.com/redirect/new_site?return=/1996-1944/19/5/964). Explore and share your feedback with us. [*close*](https://www.mdpi.com/1996-1944/19/5/964)
[Journals](https://www.mdpi.com/about/journals)
[Materials](https://www.mdpi.com/journal/materials)
[Volume 19](https://www.mdpi.com/1996-1944/19)
[Issue 5](https://www.mdpi.com/1996-1944/19/5)
[10\.3390/ma19050964](https://www.mdpi.com/1996-1944/19/5/964)
[](https://www.mdpi.com/journal/materials)
[Submit to this Journal](https://susy.mdpi.com/user/manuscripts/upload?form%5Bjournal_id%5D%3D14) [Review for this Journal](https://susy.mdpi.com/volunteer/journals/review) [Propose a Special Issue](https://www.mdpi.com/journalproposal/sendproposalspecialissue/materials)
[► ▼ Article Menu](https://www.mdpi.com/1996-1944/19/5/964)
## Article Menu
- [Academic Editor](https://www.mdpi.com/1996-1944/19/5/964#academic_editors)
[Dong-Joo Kim](https://sciprofiles.com/profile/511359?utm_source=mdpi.com&utm_medium=website&utm_campaign=avatar_name)
- [Recommended Articles](https://www.mdpi.com/1996-1944/19/5/964)
- [Related Info Link](https://www.mdpi.com/1996-1944/19/5/964#related)
- [Google Scholar](http://scholar.google.com/scholar?q=Study%20on%20Interfacial%20Crack%20of%20Piezoelectric%20Bimaterials%20Under%20Dynamic%20Loading)
- [More by Authors Links](https://www.mdpi.com/1996-1944/19/5/964#authors)
- [on DOAJ]()
- [Zhang, Y.](http://doaj.org/search/articles?source=%7B%22query%22%3A%7B%22query_string%22%3A%7B%22query%22%3A%22%5C%22Yani%20Zhang%5C%22%22%2C%22default_operator%22%3A%22AND%22%2C%22default_field%22%3A%22bibjson.author.name%22%7D%7D%7D)
- [Li, J.](http://doaj.org/search/articles?source=%7B%22query%22%3A%7B%22query_string%22%3A%7B%22query%22%3A%22%5C%22Junlin%20Li%5C%22%22%2C%22default_operator%22%3A%22AND%22%2C%22default_field%22%3A%22bibjson.author.name%22%7D%7D%7D)
- [Li, X.](http://doaj.org/search/articles?source=%7B%22query%22%3A%7B%22query_string%22%3A%7B%22query%22%3A%22%5C%22Xiangyu%20Li%5C%22%22%2C%22default_operator%22%3A%22AND%22%2C%22default_field%22%3A%22bibjson.author.name%22%7D%7D%7D)
- [Ma, J.](http://doaj.org/search/articles?source=%7B%22query%22%3A%7B%22query_string%22%3A%7B%22query%22%3A%22%5C%22Junye%20Ma%5C%22%22%2C%22default_operator%22%3A%22AND%22%2C%22default_field%22%3A%22bibjson.author.name%22%7D%7D%7D)
- [on Google Scholar]()
- [Zhang, Y.](http://scholar.google.com/scholar?q=Yani%20Zhang)
- [Li, J.](http://scholar.google.com/scholar?q=Junlin%20Li)
- [Li, X.](http://scholar.google.com/scholar?q=Xiangyu%20Li)
- [Ma, J.](http://scholar.google.com/scholar?q=Junye%20Ma)
- [on PubMed]()
- [Zhang, Y.](http://www.pubmed.gov/?cmd=Search&term=Yani%20Zhang)
- [Li, J.](http://www.pubmed.gov/?cmd=Search&term=Junlin%20Li)
- [Li, X.](http://www.pubmed.gov/?cmd=Search&term=Xiangyu%20Li)
- [Ma, J.](http://www.pubmed.gov/?cmd=Search&term=Junye%20Ma)
[Article Views](https://www.mdpi.com/1996-1944/19/5/964#metrics)
[Citations \-](https://www.mdpi.com/1996-1944/19/5/964#metrics)
- [Table of Contents](https://www.mdpi.com/1996-1944/19/5/964#table_of_contents)
Altmetric [*share* Share](https://www.mdpi.com/1996-1944/19/5/964 "Share") [*announcement* Help](https://www.mdpi.com/1996-1944/19/5/964 "Help") [*format\_quote* Cite]()
## Need Help?
### Support
Find support for a specific problem in the support section of our website.
[Get Support](https://www.mdpi.com/about/contactform)
### Feedback
Please let us know what you think of our products and services.
[Give Feedback](https://www.mdpi.com/feedback/send)
### Information
Visit our dedicated information section to learn more about MDPI.
[Get Information](https://www.mdpi.com/authors)
[*clear*]()
## JSmol Viewer
[*clear*]()
*first\_page*
[Download PDF](https://www.mdpi.com/1996-1944/19/5/964/pdf?version=1772595422)
*settings*
[Order Article Reprints](https://www.mdpi.com/1996-1944/19/5/964/reprints)
Font Type:
*Arial* *Georgia* *Verdana*
Font Size:
Aa Aa Aa
Line Spacing:
** ** **
Column Width:
** ** **
Background:
Open AccessArticle
# Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading
by
Yani Zhang
Yani Zhang
[SciProfiles](https://sciprofiles.com/profile/author/eE9UcHRPdmNIaVpXRnZoN1F2bEIzUUJWaksyNHppeHFZMzFtWktNSkp2dz0=?utm_source=mdpi.com&utm_medium=website&utm_campaign=avatar_name) [Scilit](https://scilit.com/scholars?q=Yani%20Zhang) [Preprints.org](https://www.preprints.org/search?condition_blocks=[{%22value%22:%22Yani+Zhang%22,%22type%22:%22author%22,%22operator%22:null}]&sort_field=relevance&sort_dir=desc&page=1&exact_match=true) [Google Scholar](https://scholar.google.com/scholar?q=Yani+Zhang)
1,\*,
Junlin Li
Junlin Li
[SciProfiles](https://sciprofiles.com/profile/2248912?utm_source=mdpi.com&utm_medium=website&utm_campaign=avatar_name) [Scilit](https://scilit.com/scholars?q=Junlin%20Li) [Preprints.org](https://www.preprints.org/search?condition_blocks=[{%22value%22:%22Junlin+Li%22,%22type%22:%22author%22,%22operator%22:null}]&sort_field=relevance&sort_dir=desc&page=1&exact_match=true) [Google Scholar](https://scholar.google.com/scholar?q=Junlin+Li)
1[](https://orcid.org/0000-0003-3049-8949),
Xiangyu Li
Xiangyu Li
[SciProfiles](https://sciprofiles.com/profile/author/a2hoS3c4TTIxcmNtV3hpVUlZc0x6ajRHL2l4MVVkZmhFOE5hM2hLa0J5ND0=?utm_source=mdpi.com&utm_medium=website&utm_campaign=avatar_name) [Scilit](https://scilit.com/scholars?q=Xiangyu%20Li) [Preprints.org](https://www.preprints.org/search?condition_blocks=[{%22value%22:%22Xiangyu+Li%22,%22type%22:%22author%22,%22operator%22:null}]&sort_field=relevance&sort_dir=desc&page=1&exact_match=true) [Google Scholar](https://scholar.google.com/scholar?q=Xiangyu+Li)
2 and
Junye Ma
Junye Ma
[SciProfiles](https://sciprofiles.com/profile/5071399?utm_source=mdpi.com&utm_medium=website&utm_campaign=avatar_name) [Scilit](https://scilit.com/scholars?q=Junye%20Ma) [Preprints.org](https://www.preprints.org/search?condition_blocks=[{%22value%22:%22Junye+Ma%22,%22type%22:%22author%22,%22operator%22:null}]&sort_field=relevance&sort_dir=desc&page=1&exact_match=true) [Google Scholar](https://scholar.google.com/scholar?q=Junye+Ma)
1
1
School of Applied Science, Taiyuan University of Science and Technology, Taiyuan 030024, China
2
School of Materials Science and Engineering, Taiyuan University of Science and Technology, Taiyuan 030024, China
\*
Author to whom correspondence should be addressed.
*Materials* **2026**, *19*(5), 964; https://doi.org/10.3390/ma19050964 (registering DOI)
Submission received: 8 January 2026 / Revised: 19 February 2026 / Accepted: 24 February 2026 / Published: 2 March 2026
(This article belongs to the Section [Mechanics of Materials](https://www.mdpi.com/journal/materials/sections/mechanics_materials))
[Download *keyboard\_arrow\_down*]()
[Download PDF](https://www.mdpi.com/1996-1944/19/5/964/pdf?version=1772595422)
[Download XML](https://www.mdpi.com/1996-1944/19/5/964)
[Browse Figures](https://www.mdpi.com/1996-1944/19/5/964)
[Versions Notes](https://www.mdpi.com/1996-1944/19/5/964/notes)
## Abstract
To meet the requirements of effectiveness and strength in actual engineering, based on the dynamic fracture characteristics, the dynamic propagation of orthogonal anisotropic interface cracks in piezoelectric bimaterials was analyzed. By performing Laplace transformation and Fourier transformation on the governing equations, the problem was transformed into a singular integral equation. Using the Chebyshev point method and Laplace inversion, the stress and electric displacement intensity factors at the crack tip of the orthogonal anisotropic interface were obtained. The results show that the crack length affects the dimensionless function. The longer the crack, the larger the dimensionless function. Under certain conditions, the smaller the elastic parameters, the smaller the dimensionless dynamic stress intensity factor. At the same time, the impact time also affects the dynamic crack propagation. With the passage of time, the dimensionless function first increases, then reaches a peak, and finally oscillates and converges to the static value. On this basis, the response surface method was used for analysis and prediction. The R2 value of the random forest model is 0.9886, which indicates that the model has high predictive accuracy. When the optimal values of A (
d
1
/
a
), B (
c
p
t
/
a
) and C (
c
44
(
2
)
/
c
44
(
1
)
) are 0.4045, 1.6797 and 1.9035 respectively, the stress intensity reaches its maximum value of 1.3375.
Keywords:
[integral transformation](https://www.mdpi.com/search?q=integral+transformation); [piezoelectric bimaterials](https://www.mdpi.com/search?q=piezoelectric+bimaterials); [dynamic fracture](https://www.mdpi.com/search?q=dynamic+fracture); [interfacial crack](https://www.mdpi.com/search?q=interfacial+crack); [non-dimensional function](https://www.mdpi.com/search?q=non-dimensional+function)
## 1\. Introduction
Piezoelectric materials have been widely used in advanced equipment, instruments and materials due to their electromechanical coupling properties. However, during manufacturing or usage, defects such as cracks are inevitable, which may lead to structural damage or device failure. Therefore, the crack problem of piezoelectric materials has attracted considerable attention \[[1](https://www.mdpi.com/1996-1944/19/5/964#B1-materials-19-00964),[2](https://www.mdpi.com/1996-1944/19/5/964#B2-materials-19-00964),[3](https://www.mdpi.com/1996-1944/19/5/964#B3-materials-19-00964),[4](https://www.mdpi.com/1996-1944/19/5/964#B4-materials-19-00964),[5](https://www.mdpi.com/1996-1944/19/5/964#B5-materials-19-00964)\]. In recent years, piezoelectric bimaterials have been widely applied in devices. The sudden change in materials on both sides of the interface may cause stress concentration, which is one of the important reasons for material failure \[[6](https://www.mdpi.com/1996-1944/19/5/964#B6-materials-19-00964),[7](https://www.mdpi.com/1996-1944/19/5/964#B7-materials-19-00964),[8](https://www.mdpi.com/1996-1944/19/5/964#B8-materials-19-00964)\]. Li et al. \[[9](https://www.mdpi.com/1996-1944/19/5/964#B9-materials-19-00964)\] analyzed the SH-wave scattering by interfacial cracks in piezoelectric ceramic-polymer composites. Their core objective was to establish a dynamic response analysis framework for cracks under a steady-state wave field, and they preliminarily investigated the effects of frequency, incident angle, and material properties on crack propagation. This work provides a fundamental wave-based analysis model for subsequent dynamic fracture research; however, it does not address the transient response under impact loading or the coupling behavior in complex anisotropic material systems.
Research on the behavior of interface crack propagation has achieved fruitful results. Loboda et al. conducted studies on the conductive interface cracks between one-dimensional hexagonal piezoelectric quasicrystals and ordinary piezoelectric materials, analyzing the coupling effect between the cracks and the far-field electrodes, and clarifying the energy-evolution law of piezoelectric hetero-interface cracks, providing theoretical references for the analysis of related issues \[[10](https://www.mdpi.com/1996-1944/19/5/964#B10-materials-19-00964),[11](https://www.mdpi.com/1996-1944/19/5/964#B11-materials-19-00964),[12](https://www.mdpi.com/1996-1944/19/5/964#B12-materials-19-00964)\]. Sun et al. \[[13](https://www.mdpi.com/1996-1944/19/5/964#B13-materials-19-00964)\] proposed the arbitrary-order generalized finite difference method (GFDM), providing a new numerical solution path for homogeneous and dual-material piezoelectric crack problems. The crack propagation behavior of multiple interfaces caused by circular holes under dynamic loading was studied, and the influence law of functional gradient characteristics on the dynamic response of cracks was clarified, confirming that the gradient distribution of material parameters would significantly change the stress-electric field concentration effect at the crack tip \[[14](https://www.mdpi.com/1996-1944/19/5/964#B14-materials-19-00964)\]. Based on the improved mechanical energy release rate criterion, the dynamic expansion characteristics of cracks in piezoelectric plates were revealed \[[15](https://www.mdpi.com/1996-1944/19/5/964#B15-materials-19-00964)\]. In terms of innovative analysis methods, a new interaction integral method was proposed, achieving precise analysis of dynamic impermeable interface cracks in heterogeneous piezoelectric materials, effectively improving the accuracy of field quantity calculation at the crack tip \[[16](https://www.mdpi.com/1996-1944/19/5/964#B16-materials-19-00964)\]. Liu et al. \[[17](https://www.mdpi.com/1996-1944/19/5/964#B17-materials-19-00964),[18](https://www.mdpi.com/1996-1944/19/5/964#B18-materials-19-00964),[19](https://www.mdpi.com/1996-1944/19/5/964#B19-materials-19-00964)\] based on the non-local piezoelectric theory, systematically studied the dynamic response of I-type penetrating cracks in piezoelectric materials, derived the dynamic analytical solutions for finite-penetrating and complete-penetrating cracks, and improved the theoretical system of dynamic fracture in piezoelectric materials. Afshar et al. \[[20](https://www.mdpi.com/1996-1944/19/5/964#B20-materials-19-00964)\] studied the dynamic fracture behavior of multiple intruded cracks in functional gradient piezoelectric strips \[[20](https://www.mdpi.com/1996-1944/19/5/964#B20-materials-19-00964)\]. Wünsche et al. \[[21](https://www.mdpi.com/1996-1944/19/5/964#B21-materials-19-00964)\] used the symmetrical Galerkin boundary element method to establish an analysis model for dynamic cracks in piezoelectric bodies under harmonic loading, providing an efficient solution method for piezoelectric fracture problems under harmonic loading. Yu et al. \[[22](https://www.mdpi.com/1996-1944/19/5/964#B22-materials-19-00964)\] based on the extended finite element method (XFEM), established a numerical solution framework for the dynamic crack of heterogeneous piezoelectric materials under force–electric coupling loading, providing an efficient calculation tool for interface fracture problems in complex configurations.
Through in-depth analysis of the existing research, we can find that the electromechanical coupling constitutive relationship of orthotropic anisotropic materials is more complex, and the direction dependence of their physical parameters leads to the inability to directly transplant the existing theoretical models for isotropic or functionally graded materials; at the same time, although some studies involve dynamic loads, they are mainly concentrated on harmonic loads or do not fully capture the transient response characteristics under impact loads; and the existing research mostly focuses on specific types of conductive cracks, non-conductive cracks, etc., and mainly focuses on the analysis of the influence of a single factor, without deeply exploring the coupling effect and quantitative influence laws of multiple parameters such as elastic parameters, crack length, and impact time, which is difficult to meet the systematic research requirements of this article on the dynamic expansion mechanism of interface cracks in orthotropic anisotropic piezoelectric dual materials.
In summary, although the existing research has laid a foundation for the crack problem of dual-layer materials, it has obvious limitations in adapting to the specific research object, load conditions, and analysis goals of this article.
The above research has constructed a research framework for the crack problem of piezoelectric materials from aspects such as material system, crack form, analysis method, and load type. However, the theory and methods still have obvious limitations in applying to the orthogonal anisotropic piezoelectric dual-material interface crack problem under dynamic loads. On one hand, the existing research on dual-material interface cracks mostly focuses on static/sinusoidal loads, isotropic/functional gradient piezoelectric media, or do not consider the orthogonal anisotropy characteristics of the materials. The electromechanical coupling constitutive relationship of orthogonal anisotropic piezoelectric dual-material is more complex, and the transient response laws under dynamic impact loads are significantly different from those of conventional media. The existing static theories and analysis methods under sinusoidal loads cannot be directly transplanted. On the other hand, most dynamic-crack research focuses on internal cracks of a single piezoelectric material or does not deeply explore the coupling effects of multiple parameters and quantitative influence rules. There is no complete analytical–numerical coupled solution system for the dynamic-fracture problem of orthogonal anisotropic piezoelectric dual-material interfaces, and the research on the quantitative influence of key parameters such as elastic parameters, crack length, and impact time and the interaction between parameters is still scarce. In addition, most existing studies adopt a single analytical or numerical method for crack analysis. The research of integrating the response-surface method into dynamic fracture analysis to achieve rapid prediction of multiple parameters and identification of optimal working conditions is still in the exploratory stage.
During manufacturing and usage, cracks may occur at the interface due to manufacturing level and other factors. Meanwhile, under dynamic loads, the material thickness, elastic constants, etc., of the piezoelectric bimaterials will affect the dynamic propagation of interface cracks, which may lead to device failure and losses. However, there are few articles on this topic. Moreover, the dynamic propagation of interface cracks in piezoelectric bimaterials has limitations in software simulation and experiments. Therefore, precise quantitative analysis of the dynamic fracture problem of interface cracks in piezoelectric bimaterials is of great significance for structural design and safety.
The motivation for this study stems from the urgent need to address the dynamic fracture behavior of interface cracks in piezoelectric bimaterials. This is a topic that has received limited attention in existing literature, especially concerning the complex scenario involving orthogonally anisotropic materials under dynamic loads. Piezoelectric materials play a critical role in advanced equipment, instruments, and sensors due to their unique electromechanical coupling properties. However, their brittle nature and inevitable defects, such as cracks, during manufacturing and use can lead to structural damage or even device failure. While research on static fracture problems has a certain foundation, many piezoelectric components in practical applications are often subjected to dynamic loads, making the study of dynamic fracture behavior particularly important and challenging. Currently, precise quantitative analysis of dynamic interface cracks in orthogonally anisotropic piezoelectric bimaterials, whether in theoretical analysis, numerical simulation, or experimental verification, faces significant deficiencies.
Our primary contribution lies in the innovative application and integration of advanced analytical and computational methods to deeply explore this complex problem. Specifically, we employed the powerful mathematical tools of Laplace transformation and Fourier transformation, successfully converting the governing equations that describe the dynamic expansion of orthogonally anisotropic interface cracks into a complex singular integral equation. Subsequently, we applied a high-precision Chebyshev point method for numerical discretization of this integral equation, combined with efficient Laplace inversion techniques, to accurately solve for the stress-intensity factors and electric-displacement-intensity factors at the crack tip. This combined analytical and numerical approach provides a solid foundation for understanding dynamic-fracture mechanisms. Furthermore, to further enhance analysis efficiency and predictive capability, we proactively integrated the Response Surface Methodology (RSM) and specifically developed and applied a high-precision Random Forest (Random Forest) machine learning model (with a determination coefficient R2 as high as 0.9886, demonstrating the model’s excellent goodness of fit and predictive accuracy) for effective analysis and prediction of complex multi-parameter influences.
This multi-level comprehensive method offers significant added value. It not only precisely reveals the quantitative influence laws of key factors such as crack length, elastic parameters (e.g., modulus), and impact time on dynamic fracture behavior—for example, revealing the positive correlation between crack length and the dimensionless function, as well as the inverse relationship between elastic parameters and the dimensionless dynamic stress intensity factor under certain conditions—but also captures the dynamic evolution process induced by impact time, where the dynamic crack response first increases, reaches a peak, and then oscillates and converges to the static value. Moreover, it ultimately enables the identification of optimal conditions affecting stress intensity. This in-depth and quantitative understanding greatly enhances the prediction accuracy of dynamic crack propagation mechanisms, providing a key theoretical basis for evaluating and optimizing the dynamic reliability of piezoelectric devices, thereby significantly improving the operational reliability and structural safety of piezoelectric equipment in practical engineering applications. By directly addressing the limitations of current software simulations in handling dynamic fracture problems of complex anisotropic materials, as well as the challenges of experimental research in terms of cost, scale, and dynamic loading control, this study aims to fill a critical gap in the field and provide important scientific support and technical guidance for the future design of more robust, safer, and more reliable piezoelectric smart structures.
## 2\. Establishment of the Model
As shown in [Figure 1](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f001), the thicknesses of the two materials are
d
1
and
d
2
, respectively. The crack is from
a
1
to
b
1
. With the interface direction as the X-axis and the thickness direction as the Y-axis. Let us assume that the polarization direction of the piezoelectric dual materials is the Z-axis, the structure is orthotropic, and the symmetry axes of the materials are perfectly aligned at the interface.
The governing equation for the crack in the piezoelectric bimaterial orthotropic interface is:
c
55
(
k
)
∂
2
w
∂
X
2
\+
c
44
(
k
)
∂
2
w
∂
Y
2
\+
e
15
(
k
)
∂
2
ϕ
∂
X
2
\+
e
24
(
k
)
∂
2
ϕ
∂
Y
2
\=
ρ
∂
2
w
∂
t
2
(1)
e
15
(
k
)
∂
2
w
∂
X
2
\+
e
24
(
k
)
∂
2
w
∂
Y
2
−
ε
11
(
k
)
∂
2
ϕ
∂
X
2
−
ε
22
(
k
)
∂
2
ϕ
∂
Y
2
\=
0
(2)
In the above equations,
w
\=
w
(
x
,
y
,
t
)
denotes the mechanical displacements,
φ
(
x
,
y
,
t
)
represents the electric potentials,
ρ
is the mass density,
c
55
(
k
)
and
c
44
(
k
)
are the elastic stiffness constants,
e
15
(
k
)
and
e
24
(
k
)
stand for the piezoelectric constants, and
ε
11
(
k
)
and
ε
22
(
k
)
are the dielectric constants. Here,
k
\=
1
,
2
correspond to piezoelectric material 1 and piezoelectric material 2, respectively.
Introduce the transformation
x
\=
X
−
v
t
,
y
\=
Y
, then Equations (1) and (2) becomes:
c
55
(
k
)
c
44
(
k
)
−
v
2
c
k
¯
2
1
e
15
(
k
)
c
44
(
k
)
e
24
(
k
)
c
44
(
k
)
e
15
(
k
)
c
44
(
k
)
e
24
(
k
)
c
44
(
k
)
−
ε
11
(
k
)
c
44
(
k
)
−
ε
22
(
k
)
c
44
(
k
)
∂
2
w
∂
x
2
∂
2
w
∂
y
2
∂
2
ϕ
∂
x
2
∂
2
ϕ
∂
y
2
\=
Ψ
(
k
)
0
(3)
Perform the Laplace transform on time
t
and the Fourier transform on
x
. The solution of Equation (3) can be expressed as
\[
w
\*
x
,
y
,
p
ϕ
\*
x
,
y
,
p
\]
\=
1
2
π
∫
−
∞
\+
∞
\[
χ
11
(
k
)
χ
12
(
k
)
χ
21
(
k
)
χ
22
(
k
)
\]
\[
A
11
(
k
)
e
(
3
−
2
k
)
s
γ
1
(
k
)
y
A
12
(
k
)
e
(
3
−
2
k
)
s
γ
2
(
k
)
y
\]
e
−
i
s
x
d
s
(4)
where
χ
1
η
(
k
)
\=
e
24
(
k
)
c
44
(
k
)
γ
η
(
k
)
2
−
s
2
e
15
(
k
)
c
44
(
k
)
,
χ
2
η
(
k
)
\=
γ
η
(
k
)
2
−
s
2
\[
(
c
55
(
k
)
c
44
(
k
)
−
v
2
c
k
¯
2
)
\+
2
v
p
i
c
k
¯
2
s
\+
p
2
c
k
¯
2
s
2
\]
,
η
\=
1
,
2
.
The value of
γ
η
(
k
)
can be found in [Appendix A](https://www.mdpi.com/1996-1944/19/5/964#app1-materials-19-00964).
By analyzing the model, it can be concluded that the boundary conditions that the cracks must satisfy are
σ
y
z
(
x
,
d
1
,
t
)
\=
σ
y
z
(
x
,
−
d
2
,
t
)
\=
0
(5)
w
(
x
,
0
\+
,
t
)
\=
w
(
x
,
0
−
,
t
)
,
D
y
(
x
,
0
\+
,
t
)
\=
D
y
(
x
,
0
−
,
t
)
(6)
σ
y
z
(
x
,
0
\+
,
t
)
\=
σ
y
z
(
x
,
0
−
,
t
)
\=
−
τ
0
H
(
t
)
D
y
(
x
,
0
\+
,
t
)
\=
D
y
(
x
,
0
−
,
t
)
\=
−
D
0
H
(
t
)
(7)
where
H
(
t
)
\=
0
,
t
\<
0
1
,
t
\>
0
.
By analyzing the model, it can be concluded that the boundary conditions that the cracks must satisfy are
σ
y
z
\*
(
x
,
d
1
,
p
)
\=
σ
y
z
\*
(
x
,
−
d
2
,
p
)
\=
0
(8)
w
\*
(
x
,
0
\+
,
p
)
\=
w
\*
(
x
,
0
−
,
p
)
,
D
y
\*
(
x
,
0
\+
,
p
)
\=
D
y
\*
(
x
,
0
−
,
p
)
(9)
σ
y
z
\*
(
x
,
0
\+
,
p
)
\=
σ
y
z
\*
(
x
,
0
−
,
p
)
\=
−
τ
0
p
,
D
y
\*
(
x
,
0
\+
,
p
)
\=
D
y
\*
(
x
,
0
−
,
p
)
\=
−
D
0
p
(10)
The electric displacement corresponding to the stress in Equation (4) is
σ
y
z
\*
D
y
\*
\=
1
2
π
∫
−
∞
\+
∞
(
c
44
(
k
)
e
24
(
k
)
e
15
(
k
)
−
ε
22
(
k
)
χ
11
(
k
)
χ
12
(
k
)
χ
21
(
k
)
χ
22
(
k
)
γ
1
(
k
)
e
(
3
−
2
k
)
s
γ
1
(
k
)
y
A
11
(
k
)
γ
2
(
k
)
e
(
3
−
2
k
)
s
γ
2
(
k
)
y
A
12
(
k
)
(
3
−
2
k
)
s
)
e
−
i
s
x
d
s
(11)
## 3\. Solution of the Model
In the process of crack propagation, the expansion of the crack is affected by a variety of complex factors. To more accurately describe the crack propagation behavior, introducing a density function is an effective method. Therefore, the following density function is introduced:
h
w
(
x
)
\=
d
d
x
\[
w
(
x
,
0
\+
,
p
)
−
w
(
x
,
0
−
,
p
)
\]
(12)
h
ϕ
(
x
)
\=
d
d
x
\[
ϕ
(
x
,
0
\+
,
p
)
−
ϕ
(
x
,
0
−
,
p
)
\]
(13)
Since displacement and potential are continuous, we can obtain:
∫
a
1
b
1
h
w
(
x
)
d
x
\=
0
,
x
∉
(
a
1
,
b
1
)
(14)
∫
a
1
b
1
h
ϕ
(
x
)
d
x
\=
0
,
x
∉
(
a
1
,
b
1
)
(15)
Based on the boundary conditions (8–10) and Equations (12) and (13), the following conclusion can be obtained:
χ
11
(
1
)
χ
12
(
1
)
χ
21
(
1
)
χ
22
(
1
)
A
11
B
11
−
χ
11
(
2
)
χ
12
(
2
)
χ
21
(
2
)
χ
22
(
2
)
A
12
B
12
\=
i
s
∫
a
1
b
1
h
w
(
s
)
h
ϕ
(
s
)
e
i
s
f
d
f
\=
H
1
H
2
(16)
A
11
(
1
)
A
12
(
1
)
A
11
(
2
)
A
12
(
2
)
\=
1
A
\*
A
1
B
1
A
2
A
2
A
1
\*
B
1
\*
A
2
\*
B
2
\*
H
1
H
2
(17)
The values of
A
1
,
B
1
,
A
2
,
B
2
,
A
1
\*
,
B
1
\*
,
A
2
\*
,
B
2
\*
,
A
\*
can be found in [Appendix A](https://www.mdpi.com/1996-1944/19/5/964#app1-materials-19-00964).
Substituting Equation (17) into Equation (11) yields
σ
y
z
\*
(
x
,
0
\+
,
p
)
D
y
\*
(
x
,
0
\+
,
p
)
\=
lim
y
→
0
\+
1
2
π
∫
−
∞
\+
∞
s
e
s
γ
y
A
\*
C
1
C
2
C
3
C
4
H
1
H
2
e
−
i
s
x
d
s
(18)
where
C
1
\=
∑
i
\=
1
2
(
c
44
(
1
)
χ
1
i
(
1
)
\+
e
24
(
1
)
χ
2
i
(
1
)
)
A
i
γ
i
(
1
)
,
C
2
\=
∑
i
\=
1
2
(
c
44
(
1
)
χ
1
i
(
1
)
\+
e
24
(
1
)
χ
2
i
(
1
)
)
B
i
γ
i
(
1
)
,
C
3
\=
∑
i
\=
1
2
(
e
15
(
1
)
χ
1
i
(
1
)
−
ε
22
(
1
)
χ
2
i
(
1
)
)
A
i
\*
γ
i
(
1
)
,
C
4
\=
∑
i
\=
1
2
(
e
15
(
1
)
χ
1
i
(
1
)
−
ε
22
(
1
)
χ
2
i
(
1
)
)
B
i
\*
γ
i
(
1
)
,
e
s
γ
y
\=
∑
i
\=
1
2
e
s
γ
i
(
1
)
y
On the crack surface, there is
∫
a
1
b
1
(
C
1
\*
C
2
\*
C
3
\*
C
4
\*
1
B
\*
(
f
−
x
)
\+
k
1
(
x
,
f
)
k
2
(
x
,
f
)
k
3
(
x
,
f
)
k
4
(
x
,
f
)
)
h
w
(
f
)
h
ϕ
(
f
)
d
f
\=
−
π
p
σ
0
D
0
(19)
where
k
1
(
x
,
f
)
\=
∫
0
∞
(
C
1
A
\*
−
C
1
\*
B
\*
)
sin
(
s
(
f
−
x
)
)
d
s
,
k
2
(
x
,
f
)
\=
∫
0
∞
(
C
2
A
\*
−
C
2
\*
B
\*
)
sin
(
s
(
f
−
x
)
)
d
s
,
k
3
(
x
,
f
)
\=
∫
0
∞
(
C
3
A
\*
−
C
3
\*
B
\*
)
sin
(
s
(
f
−
x
)
)
d
s
,
k
4
(
x
,
f
)
\=
∫
0
∞
(
C
4
A
\*
−
C
4
\*
B
\*
)
sin
(
s
(
f
−
x
)
)
d
s
.
The value of
C
1
\*
,
C
2
\*
,
C
3
\*
,
C
4
\*
,
B
\*
is referred to in [Appendix A](https://www.mdpi.com/1996-1944/19/5/964#app1-materials-19-00964).
To simplify the calculation, we will transform the normalized quantity into the standard form and introduce the dislocation density function
h
w
(
x
)
\=
q
1
(
x
)
,
h
ϕ
(
x
)
\=
q
2
(
x
)
,
k
l
(
x
,
f
)
\=
P
l
(
r
,
c
)
,
(
l
\=
1
,
2
,
3
,
4
)
,
m
o
\=
b
1
−
a
1
2
,
n
o
\=
b
1
\+
a
1
2
,
f
\=
x
−
n
0
m
0
,
c
\=
f
−
n
0
m
0
,
(20)
q
1
(
c
)
\=
h
1
(
c
)
1
−
c
2
π
σ
0
p
,
q
2
(
c
)
\=
h
2
(
c
)
1
−
c
2
π
D
0
p
.
According to Equation (20) and the Chebyshev placement method, Equation (19) can be transformed into
1
N
∑
Q
\=
0
N
Φ
Q
(
C
1
\*
C
2
\*
C
3
\*
C
4
\*
Ξ
\+
m
0
P
1
(
r
u
,
c
Q
)
P
2
(
r
u
,
c
Q
)
P
3
(
r
u
,
c
Q
)
P
4
(
r
u
,
c
Q
)
)
h
1
(
c
Q
)
h
2
(
c
Q
)
\=
−
1
−
1
(21)
where
Ξ
\=
m
0
B
\*
(
m
0
c
Q
\+
n
0
−
m
0
r
u
−
n
0
)
,
∑
Q
\=
0
N
Φ
Q
h
1
(
c
Q
)
\=
0
,
∑
Q
\=
0
N
Φ
Q
h
2
(
c
Q
)
\=
0
,
u
\=
1
,
2
,
…
N
,
u
\=
1
,
2
,
…
N
,
Φ
0
\=
Φ
N
\=
1
2
,
Φ
1
\=
…
\=
Φ
N
−
1
\=
1
,
r
u
\=
cos
(
2
u
−
1
)
π
2
N
,
c
Q
\=
cos
Q
π
N
(N represents a node).
## 4\. Intensity Factor
According to Equation (21), it can be obtained that
lim
x
→
a
1
−
σ
y
z
\*
(
x
,
0
,
p
)
\=
σ
0
C
1
\*
p
B
\*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
\+
σ
0
C
2
\*
p
B
\*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(22)
lim
x
→
b
1
\+
σ
y
z
\*
(
x
,
0
,
p
)
\=
σ
0
C
1
\*
p
B
\*
lim
c
→
1
\+
−
h
1
(
1
)
c
2
−
1
\+
σ
0
C
2
\*
p
B
\*
lim
c
→
1
\+
−
h
2
(
1
)
c
2
−
1
(23)
lim
x
→
a
1
−
D
y
\*
(
x
,
0
,
p
)
\=
D
0
C
3
\*
p
B
\*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
\+
D
0
C
4
\*
p
B
\*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(24)
lim
x
→
b
1
\+
D
y
\*
(
x
,
0
,
p
)
\=
D
0
C
3
\*
p
B
\*
lim
c
→
1
\+
−
h
1
(
1
)
c
2
−
1
\+
D
0
C
4
\*
p
B
\*
lim
c
→
1
\+
−
h
2
(
1
)
c
2
−
1
(25)
Using Miller and Guy, assuming that the Laplace transform of
f
(
t
)
is
f
\*
(
p
)
,
p
can be given by the following discrete points
p
\=
δ
(
β
\+
1
\+
k
)
(26)
where
δ
\>
0
,
β
\>
−
1
,
k
\=
0
,
1
,
2
,
⋯
, Then it can be concluded that
f
(
t
)
\=
∑
T
\=
0
N
B
T
P
T
(
0
,
β
)
2
e
−
δ
t
−
1
,
δ
\>
0
,
β
\>
−
1
(27)
Here,
P
T
(
0
,
β
)
is a Jacobi polynomial, and the coefficients
B
T
are determined by the following system of equations
δ
f
\*
(
δ
(
β
\+
1
\+
k
)
)
\=
∑
m
\=
0
k
k
(
k
−
1
)
⋯
k
−
(
m
−
1
)
(
k
\+
β
\+
1
)
(
k
\+
β
\+
2
)
⋯
(
k
\+
β
\+
1
\+
m
)
B
T
(28)
Based on the definitions of stress and electric displacement, the stress intensity factors
K
I
I
I
and the electric displacement intensity factors
K
I
V
at the crack tips
x
\=
a
1
and
x
\=
b
1
can be obtained as
K
I
I
I
(
a
1
)
K
I
V
(
a
1
)
\=
2
π
(
a
1
−
x
)
x
→
a
1
−
0
σ
y
z
(
x
,
0
,
t
)
D
y
(
x
,
0
,
t
)
(29)
K
I
I
I
(
b
1
)
K
I
V
(
b
1
)
\=
2
π
(
x
−
b
1
)
x
→
b
1
\+
0
σ
y
z
(
x
,
0
,
t
)
D
y
(
x
,
0
,
t
)
(30)
Here, a non-dimensional function is defined
K
\=
K
I
I
I
(
t
,
v
)
/
K
I
I
I
(
0
,
v
)
(31)
Here,
K
I
I
I
(
0
,
v
)
represents the stress intensity factor under static load conditions.
## 5\. Far-Field Forces and Electric Impact Loading
For the crack problem in piezoelectric media, as a cutting-edge branch of multi-physics field-coupled fracture mechanics, its core difficulty lies in accurately describing the complex boundary conditions at the crack tip and its surface. These conditions involve not only mechanical stress and strain but are also closely related to the electric field and electric displacement, forming a highly nonlinear boundary-value problem.
The classical strict electro-mechanical boundary conditions require that the normal displacement components at the surface of the piezoelectric medium (as the defect interface) and within the material must be continuous in the vicinity of the crack tip and at the crack interface, and the potential or electric field may also need to satisfy specific continuity or jump conditions. This strict condition stems from the principle of physical continuity and is currently the most accurate theoretical description. However, as stated in the original text, this greatly limits the scope of obtaining analytical solutions or closed-form solutions, making it extremely difficult to deeply understand and predict the fracture behavior of piezoelectric materials under complex conditions, severely restricting the application of related theories in engineering practice.
To overcome this challenge, the engineering and academic communities have developed a practical and effective approximation method. The core idea is to simplify the electric field based on the “transparency” of the crack surface. Currently, the most widely accepted and applied approach is to classify crack-boundary conditions into two categories: impermeable cracks and permeable cracks. The impermeable crack model, also known as the D-P boundary condition, is the most commonly used simplified model. It is based on a key assumption: the crack surface is completely insulated and does not allow any free charges, so the free charge density on the crack surface is zero. According to Gauss’s law, this means that the normal electric displacement component on the crack surface must be zero. Additionally, this model typically assumes that the potential on the crack surface remains constant. This boundary condition significantly simplifies the mathematical processing as it directly eliminates the influence of the crack surface as an electric field boundary, making the crack surface a “free surface” or an “equipotential surface”.
Although from a physical perspective this might be overly idealized (completely preventing the electric field from penetrating), in many engineering applications, especially when the crack propagation direction is not consistent with the electric field direction or when the crack size is much smaller than the characteristic wavelength, the D-P condition can provide results with reasonable accuracy, becoming a common benchmark for fracture mechanics analysis (such as calculating stress-intensity factors, electric-displacement-intensity factors) and structural-safety assessment. The permeable-crack model takes the opposite extreme, allowing the electric field to pass through the crack surface. In the permeable model, it is assumed that the potentials on the upper and lower surfaces of the crack are equal (i.e., the potential is continuous), or more generally, it allows for a certain jump in the potential on both sides of the crack surface, but the key is that the normal electric displacement components on both sides of the crack surface remain continuous. This means that the crack surface is no longer completely insulated but allows current to pass through. The permeable model typically better reflects the actual distribution of the electric field at the crack, especially in strong electric-field effects or specific geometric structures. However, its mathematical processing is usually more complex, and analytical solutions are more difficult to obtain.
By introducing these approximate boundary conditions, although the accuracy of the model is sacrificed to some extent, it significantly reduces the difficulty of solving the piezoelectric crack problem, enabling the use of analytical methods, semi-analytical methods to analyze fracture behavior, calculate stress intensity factors to be possible, and assess the safety of structures under complex load coupling. These approximate models provide important theoretical tools and engineering methods for understanding and predicting the failure mechanism of piezoelectric intelligent structures (such as piezoelectric sensors, detectors, and actuators).
In this study, we investigate the crack-growth problem subjected to stress and electric displacement impact loading at infinity, denoted by τ∞H(t) and D∞H(t), respectively. The boundary conditions described in Equations (5) and (7) are revised and transformed into the following forms
σ
y
z
(
x
,
d
1
,
t
)
\=
σ
y
z
(
x
,
−
d
2
,
t
)
\=
0
(32)
w
(
x
,
0
\+
,
t
)
\=
w
(
x
,
0
−
,
t
)
,
D
y
(
x
,
0
\+
,
t
)
\=
D
y
(
x
,
0
−
,
t
)
(33)
σ
y
z
(
x
,
0
\+
,
t
)
\=
σ
y
z
(
x
,
0
−
,
t
)
\=
−
τ
∞
H
(
t
)
D
y
(
x
,
0
\+
,
t
)
\=
D
y
(
x
,
0
−
,
t
)
\=
D
\*
−
D
∞
H
(
t
)
(34)
where
D
\*
is the electrical displacement in the crack plane. Take the Laplace transform of Equation (34), we obtain
σ
y
z
\*
(
x
,
0
\+
,
p
)
\=
σ
y
z
\*
(
x
,
0
−
,
p
)
\=
−
τ
∞
p
,
D
y
\*
(
x
,
0
\+
,
p
)
\=
D
y
\*
(
x
,
0
−
,
p
)
\=
−
D
∞
−
D
\*
p
(35)
### 5\.1. Impermeable Crack (D-P)
For impermeable cracks, the electric potential in the crack cavity can be reasonably neglected. This is because the dielectric constant of piezoelectric materials is usually three to four orders of magnitude higher than that of air or vacuum. Therefore, the physical quantities appearing in Equation (35)
D
\*
meet the following conditions \[[23](https://www.mdpi.com/1996-1944/19/5/964#B23-materials-19-00964)\]:
D
\*
\=
0
(36)
According to the aforementioned derivation process, the stress-intensity factor represented by
a
1
and
b
1
can be expressed as
lim
x
→
a
1
−
σ
y
z
\*
(
x
,
0
,
p
)
\=
σ
∞
C
1
\*
p
B
\*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
\+
σ
∞
C
2
\*
p
B
\*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(37)
lim
x
→
b
1
\+
σ
y
z
\*
(
x
,
0
,
p
)
\=
σ
∞
C
1
\*
p
B
\*
lim
c
→
1
\+
−
h
1
(
1
)
c
2
−
1
\+
σ
∞
C
2
\*
p
B
\*
lim
c
→
1
\+
−
h
2
(
1
)
c
2
−
1
(38)
lim
x
→
a
1
−
D
y
\*
(
x
,
0
,
p
)
\=
D
∞
C
3
\*
p
B
\*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
\+
D
∞
C
4
\*
p
B
\*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(39)
lim
x
→
b
1
\+
D
y
\*
(
x
,
0
,
p
)
\=
D
∞
C
3
\*
p
B
\*
lim
c
→
1
\+
−
h
1
(
1
)
c
2
−
1
\+
D
∞
C
4
\*
p
B
\*
lim
c
→
1
\+
−
h
2
(
1
)
c
2
−
1
(40)
In terms of the definitions of stress and electric displacement, the corresponding factors at the crack fronts can be formally derived as:
K
I
I
I
(
a
1
)
K
I
V
(
a
1
)
\=
2
π
(
a
1
−
x
)
x
→
a
1
−
0
σ
y
z
(
x
,
0
,
t
)
D
y
(
x
,
0
,
t
)
(41)
K
I
I
I
(
b
1
)
K
I
V
(
b
1
)
\=
2
π
(
x
−
b
1
)
x
→
b
1
\+
0
σ
y
z
(
x
,
0
,
t
)
D
y
(
x
,
0
,
t
)
(42)
According to Equations (41) and (42) and combined with the virtual-crack-closure technique, the corresponding energy release rate G can be obtained \[[24](https://www.mdpi.com/1996-1944/19/5/964#B24-materials-19-00964)\].
### 5\.2. Permeable Crack
For the penetrating cracks, when the material in the text degenerates into an isotropic material, the coupling degree of the mechanical and electrical loads can be expressed as
c
55
(
k
)
(
∂
2
w
∂
X
2
\+
∂
2
w
∂
Y
2
)
\+
e
15
(
k
)
(
∂
2
ϕ
∂
X
2
\+
∂
2
ϕ
∂
Y
2
)
\=
ρ
∂
2
w
∂
t
2
(43)
e
15
(
k
)
(
∂
2
w
∂
X
2
\+
∂
2
w
∂
Y
2
)
−
ε
11
(
k
)
(
∂
2
ϕ
∂
X
2
\+
∂
2
ϕ
∂
Y
2
)
\=
0
(44)
e
15
D
∞
ε
11
τ
∞
\=
λ
(45)
Note: In the previous calculation,
c
55
(
k
)
\=
c
44
(
k
)
,
e
15
(
k
)
\=
e
24
(
k
)
and
ε
11
(
k
)
\=
ε
22
(
k
)
.
For the permeable crack model, the continuity of the electric potential at the upper and lower surfaces of the crack is the core electrical boundary condition, which means that the electric potential of the media on both sides of the crack is equal everywhere on the crack surface
h
ϕ
(
x
)
\=
0
(46)
According to the aforementioned derivation process, the stress-intensity factor represented by
a
1
and
b
1
can be expressed as
lim
x
→
a
1
−
σ
y
z
\*
(
x
,
0
,
p
)
\=
σ
∞
C
1
\*
p
B
\*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
\+
σ
∞
C
2
\*
p
B
\*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(47)
lim
x
→
b
1
\+
σ
y
z
\*
(
x
,
0
,
p
)
\=
σ
∞
C
1
\*
p
B
\*
lim
c
→
1
\+
−
h
1
(
1
)
c
2
−
1
\+
σ
∞
C
2
\*
p
B
\*
lim
c
→
1
\+
−
h
2
(
1
)
c
2
−
1
(48)
lim
x
→
a
1
−
D
y
\*
(
x
,
0
,
p
)
\=
λ
ε
11
τ
∞
C
3
\*
e
15
p
B
\*
lim
c
→
−
1
−
h
1
(
−
1
)
c
2
−
1
\+
λ
ε
11
τ
∞
C
4
\*
e
15
p
B
\*
lim
c
→
−
1
−
h
2
(
−
1
)
c
2
−
1
(49)
lim
x
→
b
1
\+
D
y
\*
(
x
,
0
,
p
)
\=
λ
ε
11
τ
∞
C
3
\*
e
15
p
B
\*
lim
c
→
1
\+
−
h
1
(
1
)
c
2
−
1
\+
λ
ε
11
τ
∞
C
4
\*
e
15
p
B
\*
lim
c
→
1
\+
−
h
2
(
1
)
c
2
−
1
(50)
## 6\. Numerical Examples
To verify the above method and analyze its application, we simulated some numerical examples. In the following examples,
a
represents the half-length of the crack.
### 6\.1. The Cracks Change in Accordance with the Influence of Elastic Constants
The variation in the dimensionless function with different elastic constant ratios is shown in [Figure 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f002)a,b. When the elastic constant ratio
ς
c
\=
c
44
(
2
)
/
c
44
(
1
)
was different, the piezoelectric constant ratio and dielectric constant ratio were both 1. The crack changes under different thickness ratios of materials 1 and 2 (
d
1
/
d
2
\=
0\.0125
and
d
1
/
d
2
\=
0\.01
) were analyzed respectively in the cases of elastic constant ratio
ς
c
\=
2
([Figure 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f002)a) and
ς
c
\=
1\.5
([Figure 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f002)b).
As can be seen from the figure, the non-dimensional function decreases continuously as
d
1
/
a
increases and gradually approaches a stable value. Moreover, when the elastic modulus is large, the stable value of the non-dimensional function is also larger. At the same time, it can be observed that when the thickness of material 2 is fixed, the larger the thickness of material 1, the smaller the influence on the non-dimensional function. That is, the smaller the thickness of material 1, the greater the stress intensity generated. Therefore, in application, if we make material 1 relatively thinner, it will have higher practical application value.
### 6\.2. The Variation in Crack Under Impact Time and Crack Length
It can be clearly seen from the two-dimensional and three-dimensional curves presented in [Figure 3](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f003)a,b that the dimensionless dynamic stress intensity factor at the crack tip of the piezoelectric bimaterial interface shows obvious transient dynamic behaviors under impact load. Over time, the dimensionless function first rises rapidly, reaching a significant peak within a short period of time, then begins to gradually decline, presenting small oscillations within a certain range, and finally converges gradually to the corresponding stress intensity factor value under the static load.
This change process fully demonstrates the comprehensive influence of stress-wave propagation, reflection, superposition, and energy dissipation under impact load on the mechanical field at the crack tip and also reflects the significant differences between dynamic fracture behavior and quasi-static fracture behavior. The appearance of the peak indicates that in the initial stage of impact, the stress wave focuses strongly at the crack tip, resulting in a significant enhancement of the local mechanical and electric field coupling effect, while the oscillation convergence in the later stage indicates that the system gradually stabilizes and the dynamic effect gradually decays.
At the same time, the numerical results also show that, while keeping the thickness of material 1 unchanged, the crack length has a significant positive regulatory effect on the dimensionless function. As the crack length increases, the dimensionless dynamic stress intensity factor shows a significant increasing trend overall, which means that the expansion of the crack size will further intensify the stress concentration at the interface and increase the risk of dynamic expansion and unstable fracture of the crack.
These laws not only reveal the intrinsic relationship between geometric parameters, load forms and dynamic fracture behavior, but also provide important theoretical basis and data support for the strength design, damage identification and safety assessment of piezoelectric composite material structures under complex dynamic conditions such as impact and vibration.
### 6\.3. Impermeable and Permeable Cracks
The variation in the normalized energy-release rate under far-field stress loading is shown in [Figure 4](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f004). The dynamic nondimensional energy-release rate can be expressed as
G
/
G
r
(
G
r
is the dynamic energy-release rate for a stationary crack). [Figure 4](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f004) shows that when the far-field shear stress
τ
∞
increases from zero, the normalized energy release rate
G
/
G
r
always shows an upward trend. However, when
τ
∞
decreases from zero (i.e., in a negative shear stress state),
G
/
G
r
first decreases and then increases. This phenomenon demonstrates that the effects of positive and negative shear stresses on crack propagation have significant asymmetry: positive shear stress can continuously increase the energy release rate, thereby promoting crack propagation, while negative shear stress will reduce the energy release rate within a certain small range, manifesting as inhibition or hindrance of crack development. However, as the absolute value of negative shear stress further increases, the energy release rate will rise again. Thus, stress loading can either promote crack propagation or have interference or blocking effects, the specific manifestation depending on the direction and magnitude of the far-field stress. The above change pattern is completely consistent with the theoretical analysis results established by Pak \[[23](https://www.mdpi.com/1996-1944/19/5/964#B23-materials-19-00964)\], further verifying the rationality of the calculation model and the shear fracture theory presented in this paper.
In order to conduct a more in-depth study on the dynamic crack propagation under the coupling impact of mechanics and electricity, [Figure 5](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f005) presents the evolution of the normalized dynamic strength factor under different electromagnetic–mechanical coupling coefficients λ. It is clearly observable that as the parameter ct/(2a) related to time increases, the normalized dynamic strength factor
F
0
first rises rapidly, reaching a significant peak, then gradually decreases and eventually stabilizes at a constant value. This characteristic trend, including the sharp rise, peak response, and final steady-state behavior, is highly consistent with the dynamic fracture response studied by Wang \[[25](https://www.mdpi.com/1996-1944/19/5/964#B25-materials-19-00964)\].
### 6\.4. Response-Surface Statistical Analysis
This study employed Response Surface Methodology (RSM), a statistical technique that integrates experimental design, modeling, and optimization techniques, aiming to systematically explore and quantify the influence of multiple controllable factors on a specific response variable, and ultimately determine the conditions for obtaining the optimal response. In this study, we specifically focused on the influence pattern of three key factors—Factor A (
d
1
/
a
), Factor B (
c
p
t
/
a
), and Factor C (
c
44
(
2
)
/
c
44
(
1
)
)—on the Stress Intensity Factor R1 (SIF) response variable.
First, to efficiently and comprehensively obtain data points for modeling, we adopted the Box–Behnken Design (BBD) experimental design method provided by Design-Expert software V22. BBD is a commonly used response surface design method, particularly suitable for situations where it is necessary to estimate the interaction and quadratic effects between factors.
Next, we used the statistical module of Design-Expert to perform a quadratic polynomial regression fit on the data. To assess the reliability and explanatory power of the established model, we conducted rigorous statistical analyses. The model-fitting goodness was measured by the coefficient of determination R2 (R-squared) and the adjusted coefficient of determination Adjusted R2. The R2 value represents the proportion of the total variation in the response variable explained by the model; the closer it is to 1, the better the model fit. Adjusted R2 considers the number of independent variables in the model and adjusts R2, providing a more realistic reflection of the model’s actual predictive ability. At the same time, we performed an Analysis of Variance (ANOVA). Analysis of Variance (ANOVA) decomposes the contribution of independent variables to the total variation to assess the impact of each factor on the response variable (Stress Intensity Factor R1). The significance of each factor on the response value is mainly achieved by comparing the variance between groups with that within groups, thereby evaluating the accuracy and reliability of the model. The significance of the model is tested by analyzing the sources of error, and the F-value and p\-value are combined for joint judgment. The larger the F value and the smaller the p value, the more significant the model is. When the p value is greater than 0.05, the model is not significant. When the p\-value is less than 0.05, the model is significant, especially when the p\-value is less than 0.0001, the model shows extremely high significance. According to the data in [Table 1](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t001), the p value is less than 0.0001, so the model is significant. The ANOVA results not only validated the statistical significance of the entire model (through F-test, typically requiring a p\-value \< 0.05) but also tested the significance of individual coefficients (linear, quadratic, interaction terms) (through t\-test, typically requiring a p\-value \< 0.05 or 0.01), helping us understand which factors and their interactions significantly affect SIF. Residual analysis was also conducted to further validate the model’s predictive ability and the reasonableness of its assumptions.
Finally, based on the established and validated quadratic response surface model, we utilized the optimization tools built into the Design-Expert software to perform optimization calculations on the combination of Factors A, B, and C. The optimization objective was set to maximize the Stress Intensity Factor (SIF) response value R1. The optimization process considered all significant factor effects (including linear, quadratic, and interaction effects) and may have set operational constraints on the factors (such as factors not exceeding their physical or process-allowable ranges). The optimization algorithm in Design-Expert searched the factor space to find the optimal combination of Factors A, B, and C that maximizes SIF under the constraint conditions. This optimal combination not only provides clear guidance for practical applications but also intuitively demonstrates the core value of RSM—namely, finding optimal operating conditions from complex factor interactions through mathematical modeling and optimization.
The results of the analysis of variance (as shown in [Table 1](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t001) and [Table 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t002)) show that the F value of the regression model is 706.4, and the p values of each model are all less than 0.05, indicating the significance of the model. Among them, the p values of A-A, B-B, A2, B2, and C2 are all less than 0.0001, indicating that the model is extremely significant and highly reliable, and the experimental design is reasonable. Therefore, the results of the analysis of variance indicate that the fitting model is accurate and statistically significant.
According to the data, a second-order response surface mathematical model was obtained through analysis. The second-order response surface regression model regarding the three variables A, B, and C, as well as R1, is shown in Equation (51).
R
1
\=
−
12\.2865
\+
0\.8313
A
\+
1\.5650
B
\+
14\.0871
C
\+
0\.1817
AB
\+
0\.1818
AC
−
0\.0691
BC
−
0\.8666
A
2
−
0\.4544
B
2
−
4\.07760
C
2
(51)
[Figure 6](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f006) depicts the relationship between the predicted values and the actual values. The horizontal axis represents the actual values, and the vertical axis represents the predicted values. Under ideal conditions, if the predicted value is equal to the actual value, it lies on the black straight line. Different data points are represented by squares of different colors. By observing the image, it can be found that the vast majority of data points are concentrated near the black straight line, which directly indicates that the predicted values are close to the actual values. This also shows that the established prediction model has good accuracy and reliability, can predict the relevant values relatively accurately, and can effectively reflect the true relationship between the variables.
It can be seen from [Table 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t002) that the interaction terms AC and BC are statistically insignificant. This does not mean that there are no interaction effects at the physical level, but reflects that the strengths of these interactions are weaker compared with the main effects of factors A, B, and C. During the experiments, the coupled electromechanical effects corresponding to AC and BC do exist objectively, but their influences on the response variables (e.g., interfacial electromechanical coupling strength, stress–electric coupling distribution at the crack tip, etc.) are overshadowed by the dominant main effects. For example, A (thickness of Material 1/crack length) directly determines the constraint effect of Material 1 on crack propagation, and C (ratio of elastic constants of Material 2 to Material 1) mainly affects the mechanical matching across the interface. Their independent effects on the electromechanical coupling are more significant, whereas their synergistic regulation only induces marginal changes. In addition, we have verified the residual distribution of the model and found no systematic deviation, indicating that the insignificance of AC and BC does not result from improper model specification, but truly reflects the weak intensity of these interactions.
From the perspective of physical mechanism, the coupled electromechanical effects at the interface are essentially a complex process involving multi-factor synergy. According to the variable definitions, the synergistic pathways of AC and BC are inherently constrained by the geometrical parameters and mechanical properties of the materials, limiting their interactive contributions to a small range. Specifically, for the interaction term BC, the variation in crack propagation length is directly related to the stress concentration at the crack tip, while the elastic constant ratio determines the stress-transfer efficiency across the interface. During crack propagation, the regulation of stress transfer by the elastic constant ratio is overshadowed by the stress-concentration effect at the crack tip, producing only weak local effects that cannot be manifested as statistically significant changes in the macroscopic response variables. For the interaction term AC, their synergistic effect theoretically affects the distribution of interfacial electromechanical coupling. However, when the thickness of Material 1 is sufficient to form a stable load-bearing structure, the contribution of the interaction to the response variable becomes weak. Due to its small amplitude, it does not reach the threshold of statistical significance in the macroscopic response model, thus showing statistical insignificance.
In summary, the statistical insignificance of the interaction terms AC and BC is a reasonable reflection of the coupled electromechanical effects at the interface under the constraints of specific geometrical parameters and mechanical properties, rather than a contradiction with physical reality. The key to reconciling them is to clarify that “statistical significance” depends on the detectability of effects within the experimental design and modeling framework, while “physical existence” refers to the inherent synergistic behavior of the factors. These two concepts are complementary rather than mutually exclusive.
Through response-surface analysis of the parameters, A 3D response-surface map (as shown in [Figure 7](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f007)) and a contour map (as shown in [Figure 7](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f007)) were obtained. As shown in [Figure 7](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f007), in the three-dimensional response-surface plot, the horizontal axis represents parameters A and B respectively, and the vertical axis represents the corresponding response values. The shape and gradient changes in the response surface intuitively reflect the variation pattern of the response values under the combined effect of the two factors: the red area corresponds to higher response values, the blue area corresponds to lower response values, and the degree of color gradient clearly shows the increase or decrease range of the response values.
In the contour plot, the same pattern is also presented, that is, the color that is closer to red indicates a higher response value, and the color that is closer to blue indicates a lower response value. Overall, the response-surface plot and the contour plot can intuitively, comprehensively and quantitatively reveal the spatial distribution characteristics, interaction influence rules and extreme distribution areas of the response values with respect to parameters A and B, and can provide reliable visual evidence and theoretical support for the subsequent identification of key parameters, analysis of the influencing mechanism and selection of the optimal parameter range.
In the three-dimensional response-surface plot, the horizontal axis is A and C, and the vertical axis R1 represents the obtained response values (as shown in [Figure 8](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f008)). This surface is used to reflect the continuous pattern of response values changing with the two independent variables. The steeper the surface, the stronger the influence of the corresponding independent variable on the response value in that area. If the surface is flat, the influence of the variable on the response value is relatively weak. For the contour plot, the color represents the gradual change (green-yellow-red), indicating high or low response values. The red area has the highest response value, while the green area has the lowest response value. Points on the same curve have the same response value. The denser the curve, the more intense the change in the response value in that area. Within the 1.5 range, the contours are elliptical, indicating a significant interaction between the two independent variables and the change in response values.
[Figure 9](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f009) shows the three-dimensional response-surface plot and contour plot corresponding to the joint variation in independent variables B and C. From the figure, it can be seen that the response value shows a continuous and smooth trend as variables B and C change. The color in the figure gradually transitions from blue to red, intuitively reflecting the distribution pattern of the response value from low to high. The response surface as a whole presents a morphological characteristic of rapid increase first and then tending to be stable. This indicates that when B and C increase to a certain range, the growth rate of the response value gradually slows down and tends to stabilize.
Further analysis reveals that when the value of variable B is within the range of 1.5 to 2, the response value R1 can reach the optimal range of 1.8 to 2. At the same time, from the contour plot, a significant stretched elliptical distribution feature can be clearly observed. This feature directly indicates that there is a significant interaction between independent variables B and C, and their influence on the response value is not independent of each other but shows a clear coupling effect.
Depicted in [Figure 10](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f010) are the response values of A, B, and C under the maximum stress intensity factor condition. From the figure, it can be seen that the optimal values of A, B and C are 0.4045, 1.6797 and 1.9035 respectively, and the stress intensity factor reaches its maximum value of 1.3375 at these values. Additionally, it is found that the interaction between B and C is the most significant. The combination of B and C has a more significant impact on the response value than the combinations of other factors.
[Figure 10](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f010) shows the parameter combination when the stress-intensity factor reaches its global maximum: the corresponding values of variables A, B, and C are 0.4045, 1.6797, and 1.9035, respectively, at which the peak stress intensity factor is 1.3375. Meanwhile, the response surface analysis reveals that the interaction effect between factors B and C is statistically significant (p \< 0.05), and their coupled influence on the stress intensity factor is significantly stronger than that of other factor combinations (see [Table 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t002) for the significance analysis of interaction terms).
It should be clearly defined that the term “optimal parameters” here refers only to the mathematical extreme solution from the response-surface model, specifically the parameter point that maximizes the stress-intensity factor, rather than the optimal design scheme based on structural safety criteria. As a key parameter characterizing the stress-field intensity at the crack tip, the magnitude of the stress-intensity factor directly determines the crack initiation threshold and propagation rate. A higher stress-intensity factor indicates a more severe stress concentration at the crack tip, a greater probability of brittle fracture or fatigue failure, and lower structural safety. Therefore, the parameter combination shown in [Figure 9](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f009) corresponds to the most unfavorable service condition of the structure, whose core value is to provide a boundary reference for dangerous working conditions in engineering design, rather than a recommended safe parameter range.
In practical structural safety design and fracture control, the high-risk parameter combination shown in [Figure 9](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f009) should be avoided. We have supplemented the distinction between mathematically optimal extreme values and engineering safety-optimized solutions in the revised manuscript to prevent confusion of academic concepts, ensure the consistency between research conclusions and engineering practice, and thus provide a more valuable reference for structural safety design.
## 7\. Conclusions
This paper conducts a systematic study on the dynamic propagation of orthogonal anisotropic interface cracks in piezoelectric bimaterial structures. By introducing Laplace transform and Fourier transform to convert the governing equations and combining the Chebyshev collocation method to discretize the original problem into algebraic equations, the dynamic variation laws of the stress intensity factor at the crack tip and the electric displacement intensity factor are obtained after numerical inverse transformation. Based on this, the influence mechanism of key parameters on the fracture behavior is analyzed by using the response surface method, and the following main conclusions are obtained:
(1)
The elastic parameters have a significant regulatory effect on the dynamic propagation of interface cracks. Within a certain range, the larger the elastic parameters, the higher the dimensionless dynamic stress-intensity factor at the crack tip, indicating that the elastic properties of the material directly affect the degree of stress concentration and the driving ability of the crack.
(2)
The impact-load application time significantly changes the dynamic response characteristics of the crack. With the passage of time, the dimensionless stress-intensity factor shows a pattern of first increasing, reaching a peak, and then oscillating and converging, eventually approaching the corresponding value under static load, reflecting the combined effect of stress wave propagation, reflection, and energy dissipation.
(3)
The crack length is an important geometric factor affecting the fracture risk. The larger the crack length, the more significant the stress concentration at the tip, and the corresponding dimensionless dynamic stress-intensity factor increases. The structural fracture failure risk is higher.
(4)
For non-permeable cracks, the dual effect of forward shear stress promoting crack propagation and negative shear stress first inhibiting and then driving the crack propagation has verified the rationality of the shear-fracture calculation model established in this paper. For the dimensionless dynamic strength factors corresponding to different force–electric coupling coefficients λ, they all show a typical characteristic of rapid increase—reaching a peak—slow decrease—approaching stability as
c
p
t
/
(
2
a
)
changes. This further proves that the model in this paper can accurately describe the transient fracture-mechanics behavior of piezoelectric materials under coupled loads.
(5)
The constructed prediction model has high accuracy and strong reliability, with a determination coefficient R2 of 0.9886. The real values and predicted values are highly consistent, providing an efficient and accurate data-driven method for the rapid prediction of the stress-intensity factor at the crack tip.
(6)
The key parameter combinations that maximize the stress-intensity factor were obtained through response surface analysis: A = 0.4045, B = 1.6797, C = 1.9035, corresponding to a peak stress intensity factor of 1.3375. Significance analysis indicates that there is a significant interaction between factors B and C (p \< 0.05), and their coupling effect is stronger than other parameter combinations.
(7)
The “optimal parameter combination” obtained in this paper is only a mathematical extremum point. Even the parameter combination that maximizes the stress-intensity factor corresponds to the most dangerous service condition of the structure, rather than an engineering safety optimization scheme. This result can be used to identify high-risk parameter intervals, providing important theoretical references for the fracture prevention, safety design, and boundary determination of dangerous conditions of engineering structures.
Overall, the analytical–numerical combined analysis framework established in this paper not only enriches the dynamic fracture theory of piezoelectric bimaterial interface cracks but also provides feasible methods for the safety assessment, parameter optimization, and failure prevention of intelligent materials and structures. It has certain theoretical value and practical significance for promoting the reliable application of piezoelectric structures in engineering fields. In future research, we will incorporate the thickness ratio (d1/d2) as a key parameter into our analysis and conduct systematic quantitative studies using appropriate numerical simulation software. This extension aims to more comprehensively reveal the influence law of the thickness ratio on the dynamic fracture behavior of piezoelectric bimaterial interface cracks, further improving the completeness and engineering applicability of the proposed theoretical framework.
## Author Contributions
Conceptualization, Y.Z. and J.L.; methodology, J.L.; software, J.L.; validation, Y.Z., J.L. and X.L.; formal analysis, X.L.; investigation, J.M.; resources, Y.Z. and J.M.; data curation, Y.Z.; writing—original draft preparation, Y.Z.; writing—review and editing, X.L.; visualization, J.M.; supervision, J.L.; project administration, Y.Z.; funding acquisition, J.M. and J.L. All authors have read and agreed to the published version of the manuscript.
## Funding
This research was funded by National Natural Science Foundation of China grant number 12301591, and also by the Shanxi Provincial Key Research and Development Project grant number 202102090301027. The APC was funded by the Taiyuan University of Science and Technology doctoral research start-up fund grant number 20252008.
## Data Availability Statement
The original contributions presented in this study are included in the article, further inquiries can be directed to the corresponding author.
## Conflicts of Interest
The authors declare no conflicts of interest.
## Appendix A
ϖ
1
(
k
)
\=
(
−
ε
22
(
k
)
c
44
(
k
)
(
(
c
55
(
k
)
c
44
(
k
)
−
v
2
c
k
¯
2
)
s
2
\+
2
v
p
s
i
c
k
¯
2
\+
p
2
c
k
¯
2
)
−
2
e
24
(
k
)
e
15
(
k
)
s
2
c
44
(
k
)
c
44
(
k
)
−
ε
11
(
k
)
s
2
c
44
(
k
)
)
/
(
(
ε
22
(
k
)
c
44
(
k
)
\+
(
e
24
(
k
)
c
44
(
k
)
)
2
)
,
ϖ
2
(
k
)
\=
(
ε
11
(
k
)
c
44
(
k
)
(
(
c
55
(
k
)
c
44
(
k
)
−
v
2
c
k
¯
2
)
s
4
\+
2
v
p
s
3
i
c
k
¯
2
\+
p
2
s
2
c
k
¯
2
)
\+
e
15
(
k
)
e
15
(
k
)
s
4
c
44
(
k
)
c
44
(
k
)
)
/
(
(
ε
22
(
k
)
c
44
(
k
)
\+
(
e
24
(
k
)
c
44
(
k
)
)
2
)
,
ϑ
1
(
k
)
\=
−
72
ϖ
1
(
k
)
ϖ
2
(
k
)
−
2
ϖ
1
(
k
)
3
54
\+
(
72
ϖ
1
(
k
)
ϖ
2
(
k
)
−
2
ϖ
1
(
k
)
3
)
2
2916
−
(
12
ϖ
2
(
k
)
\+
ϖ
1
(
k
)
2
)
3
729
3
,
ϑ
1
(
k
)
\=
−
72
ϖ
1
(
k
)
ϖ
2
(
k
)
−
2
ϖ
1
(
k
)
3
54
−
(
72
ϖ
1
(
k
)
ϖ
2
(
k
)
−
2
ϖ
1
(
k
)
3
54
)
2
−
(
12
ϖ
2
(
k
)
\+
ϖ
1
(
k
)
2
9
)
3
3
,
γ
1
(
k
)
\=
−
−
2
ϖ
1
(
k
)
/
3
\+
ϑ
1
(
k
)
\+
ϑ
2
(
k
)
\+
−
4
ϖ
1
(
k
)
/
3
−
ϑ
1
(
k
)
−
ϑ
2
(
k
)
2
,
γ
2
(
k
)
\=
−
−
2
ϖ
1
(
k
)
/
3
\+
ϑ
1
(
k
)
\+
ϑ
2
(
k
)
−
−
4
ϖ
1
(
k
)
/
3
−
ϑ
1
(
k
)
−
ϑ
2
(
k
)
2
,
γ
3
(
k
)
\=
−
2
ϖ
1
(
k
)
/
3
\+
ϑ
1
(
k
)
\+
ϑ
2
(
k
)
\+
−
4
ϖ
1
(
k
)
/
3
−
ϑ
1
(
k
)
−
ϑ
2
(
k
)
2
,
γ
4
(
k
)
\=
−
2
ϖ
1
(
k
)
/
3
\+
ϑ
1
(
k
)
\+
ϑ
2
(
k
)
−
−
4
ϖ
1
(
k
)
/
3
−
ϑ
1
(
k
)
−
ϑ
2
(
k
)
2
,
A
\*
\=
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
(
−
χ
11
(
1
)
χ
22
(
2
)
\+
χ
12
(
2
)
χ
21
(
1
)
\+
χ
12
(
1
)
χ
22
(
2
)
−
χ
12
(
2
)
χ
22
(
1
)
)
\+
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
(
−
χ
11
(
2
)
χ
21
(
1
)
\+
χ
11
(
1
)
χ
21
(
2
)
−
χ
12
(
1
)
χ
21
(
2
)
\+
χ
11
(
2
)
χ
22
(
1
)
)
,
A
1
\=
−
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
χ
22
(
2
)
\+
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
χ
21
(
2
)
,
B
1
\=
−
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
χ
11
(
2
)
\+
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
χ
12
(
2
)
,
A
2
\=
−
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
χ
21
(
2
)
\+
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
χ
22
(
2
)
,
B
2
\=
−
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
χ
12
(
2
)
\+
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
χ
11
(
2
)
,
A
1
\*
\=
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
(
χ
21
(
1
)
−
χ
22
(
1
)
)
,
B
1
\*
\=
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
e
s
γ
2
(
2
)
d
2
(
χ
12
(
1
)
−
χ
11
(
1
)
)
,
A
2
\*
\=
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
(
χ
22
(
1
)
−
χ
21
(
1
)
)
,
B
2
\*
\=
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
e
s
γ
1
(
2
)
d
2
(
χ
11
(
1
)
−
χ
12
(
1
)
)
,
C
1
\*
\=
(
(
c
44
(
1
)
χ
11
(
1
)
\+
e
24
(
1
)
χ
21
(
1
)
)
γ
1
(
1
)
−
(
c
44
(
1
)
χ
12
(
1
)
\+
e
24
(
1
)
χ
22
(
1
)
)
γ
2
(
1
)
)
(
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
χ
21
(
2
)
−
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
χ
22
(
2
)
)
,
C
2
\*
\=
(
(
c
44
(
1
)
χ
11
(
1
)
\+
e
24
(
1
)
χ
21
(
1
)
)
γ
1
(
1
)
−
(
c
44
(
1
)
χ
12
(
1
)
\+
e
24
(
1
)
χ
22
(
1
)
)
γ
2
(
1
)
)
(
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
χ
12
(
2
)
−
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
χ
11
(
2
)
)
,
C
3
\*
\=
(
(
e
15
(
1
)
χ
11
(
1
)
−
ε
22
(
1
)
χ
21
(
1
)
)
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
−
(
e
15
(
1
)
χ
12
(
1
)
−
ε
22
(
1
)
χ
22
(
1
)
)
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
)
γ
1
(
1
)
γ
2
(
2
)
(
χ
21
(
1
)
−
χ
22
(
1
)
)
,
C
4
\*
\=
(
(
e
15
(
1
)
χ
11
(
1
)
−
ε
22
(
1
)
χ
21
(
1
)
)
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
−
(
e
15
(
1
)
χ
12
(
1
)
−
ε
22
(
1
)
χ
22
(
1
)
)
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
)
γ
1
(
1
)
γ
2
(
2
)
(
χ
12
(
1
)
−
χ
11
(
1
)
)
,
B
\*
\=
(
c
44
(
2
)
χ
11
(
2
)
\+
e
24
(
2
)
χ
21
(
2
)
)
γ
1
(
2
)
(
−
χ
11
(
1
)
χ
22
(
2
)
\+
χ
12
(
2
)
χ
21
(
1
)
\+
χ
12
(
1
)
χ
22
(
2
)
−
χ
12
(
2
)
χ
22
(
1
)
)
\+
(
c
44
(
2
)
χ
12
(
2
)
\+
e
24
(
2
)
χ
22
(
2
)
)
γ
2
(
2
)
(
−
χ
11
(
2
)
χ
21
(
1
)
\+
χ
11
(
1
)
χ
21
(
2
)
−
χ
12
(
1
)
χ
21
(
2
)
\+
χ
11
(
2
)
χ
22
(
1
)
)
.
## References
1. Singh, V.; Jangid, K.; Bui, T.Q. A study of strain and electric field gradient effects on two collinear cracks in an arbitrary polarized piezoelectric layer. Eur. J. Mech.-A/Solids **2025**, 114, 105723. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=A+study+of+strain+and+electric+field+gradient+effects+on+two+collinear+cracks+in+an+arbitrary+polarized+piezoelectric+layer&author=Singh,+V.&author=Jangid,+K.&author=Bui,+T.Q.&publication_year=2025&journal=Eur.+J.+Mech.-A/Solids&volume=114&pages=105723&doi=10.1016/j.euromechsol.2025.105723)\] \[[CrossRef](https://doi.org/10.1016/j.euromechsol.2025.105723)\]
2. Mu, X.; Zhu, Z.; Zhang, L.; Gao, Y. The Green’s functions of two-dimensional piezoelectric quasicrystal semi-infinite crack and finite crack. Appl. Math. Model. **2025**, 137, 115732. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=The+Green%E2%80%99s+functions+of+two-dimensional+piezoelectric+quasicrystal+semi-infinite+crack+and+finite+crack&author=Mu,+X.&author=Zhu,+Z.&author=Zhang,+L.&author=Gao,+Y.&publication_year=2025&journal=Appl.+Math.+Model.&volume=137&pages=115732&doi=10.1016/j.apm.2024.115732)\] \[[CrossRef](https://doi.org/10.1016/j.apm.2024.115732)\]
3. Li, Y.; Yan, S.; Zhao, M.; Ren, J. Fracture Analysis of Planar Cracks in 3D Thermal Piezoelectric Semiconductors. Int. J. Mech. Sci. **2024**, 273, 109212. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Fracture+Analysis+of+Planar+Cracks+in+3D+Thermal+Piezoelectric+Semiconductors&author=Li,+Y.&author=Yan,+S.&author=Zhao,+M.&author=Ren,+J.&publication_year=2024&journal=Int.+J.+Mech.+Sci.&volume=273&pages=109212&doi=10.1016/j.ijmecsci.2024.109212)\] \[[CrossRef](https://doi.org/10.1016/j.ijmecsci.2024.109212)\]
4. Tan, Y.; Rao, W.; Wan, K.; Peng, K.; Zhao, J.; Li, X. Phase-field model for fatigue crack growth in piezoelectrics: Energetically consistent boundary condition. Int. J. Solids Struct. **2025**, 316, 113378. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Phase-field+model+for+fatigue+crack+growth+in+piezoelectrics:+Energetically+consistent+boundary+condition&author=Tan,+Y.&author=Rao,+W.&author=Wan,+K.&author=Peng,+K.&author=Zhao,+J.&author=Li,+X.&publication_year=2025&journal=Int.+J.+Solids+Struct.&volume=316&pages=113378&doi=10.1016/j.ijsolstr.2025.113378)\] \[[CrossRef](https://doi.org/10.1016/j.ijsolstr.2025.113378)\]
5. Nourazar, M.; Yang, W.; Chen, Z. Fracture analysis of a curved crack in a piezoelectric plane under general thermal loading. Eng. Fract. Mech. **2023**, 284, 109208. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Fracture+analysis+of+a+curved+crack+in+a+piezoelectric+plane+under+general+thermal+loading&author=Nourazar,+M.&author=Yang,+W.&author=Chen,+Z.&publication_year=2023&journal=Eng.+Fract.+Mech.&volume=284&pages=109208&doi=10.1016/j.engfracmech.2023.109208)\] \[[CrossRef](https://doi.org/10.1016/j.engfracmech.2023.109208)\]
6. Cao, T.; Zhang, Y.; Qin, T. Research on interaction of two straight cracks embedded in different spaces of piezoelectric biomaterial. Eng. Anal. Bound. Elem. **2025**, 179, 106401. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Research+on+interaction+of+two+straight+cracks+embedded+in+different+spaces+of+piezoelectric+biomaterial&author=Cao,+T.&author=Zhang,+Y.&author=Qin,+T.&publication_year=2025&journal=Eng.+Anal.+Bound.+Elem.&volume=179&pages=106401&doi=10.1016/j.enganabound.2025.106401)\] \[[CrossRef](https://doi.org/10.1016/j.enganabound.2025.106401)\]
7. Yu, H.; Zhu, S.; Ma, H.; Wang, J. Interface crack analysis of piezoelectric laminates considering initial strain. Int. J. Mech. Sci. **2024**, 271, 109104. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Interface+crack+analysis+of+piezoelectric+laminates+considering+initial+strain&author=Yu,+H.&author=Zhu,+S.&author=Ma,+H.&author=Wang,+J.&publication_year=2024&journal=Int.+J.+Mech.+Sci.&volume=271&pages=109104&doi=10.1016/j.ijmecsci.2024.109104)\] \[[CrossRef](https://doi.org/10.1016/j.ijmecsci.2024.109104)\]
8. Zhu, S.; Yu, H.; Guo, L. Analysis of an interfacial crack between two nonhomogeneous piezoelectric materials using a new domain-independent interaction integral. Compos. Struct. **2024**, 331, 117873. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Analysis+of+an+interfacial+crack+between+two+nonhomogeneous+piezoelectric+materials+using+a+new+domain-independent+interaction+integral&author=Zhu,+S.&author=Yu,+H.&author=Guo,+L.&publication_year=2024&journal=Compos.+Struct.&volume=331&pages=117873&doi=10.1016/j.compstruct.2024.117873)\] \[[CrossRef](https://doi.org/10.1016/j.compstruct.2024.117873)\]
9. Zhang, Y.; Li, J.; Xie, X. SH-wave scattering by the interface crack of piezoelectric ceramic polymer composites. J. Eng. Res. **2024**, 12, 259–267. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=SH-wave+scattering+by+the+interface+crack+of+piezoelectric+ceramic+polymer+composites&author=Zhang,+Y.&author=Li,+J.&author=Xie,+X.&publication_year=2024&journal=J.+Eng.+Res.&volume=12&pages=259%E2%80%93267&doi=10.1016/j.jer.2023.08.003)\] \[[CrossRef](https://doi.org/10.1016/j.jer.2023.08.003)\]
10. Zhao, S.; Govorukha, V.; Sheveleva, A.; Loboda, V. On energy release rate for an electrically permeable interface crack between two different 1D hexagonal piezoelectric QCs. Procedia Struct. Integr. **2023**, 51, 69–75. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=On+energy+release+rate+for+an+electrically+permeable+interface+crack+between+two+different+1D+hexagonal+piezoelectric+QCs&author=Zhao,+S.&author=Govorukha,+V.&author=Sheveleva,+A.&author=Loboda,+V.&publication_year=2023&journal=Procedia+Struct.+Integr.&volume=51&pages=69%E2%80%9375&doi=10.1016/j.prostr.2023.10.069)\] \[[CrossRef](https://doi.org/10.1016/j.prostr.2023.10.069)\]
11. Loboda, V.; Sheveleva, A.; Komarov, O.; Chapelle, F.; Lapusta, Y. Arbitrary number of electrically permeable cracks on the interface between two one-dimensional piezoelectric quasicrystals with piezoelectric effect. Eng. Fract. Mech. **2022**, 276, 108878. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Arbitrary+number+of+electrically+permeable+cracks+on+the+interface+between+two+one-dimensional+piezoelectric+quasicrystals+with+piezoelectric+effect&author=Loboda,+V.&author=Sheveleva,+A.&author=Komarov,+O.&author=Chapelle,+F.&author=Lapusta,+Y.&publication_year=2022&journal=Eng.+Fract.+Mech.&volume=276&pages=108878&doi=10.1016/j.engfracmech.2022.108878)\] \[[CrossRef](https://doi.org/10.1016/j.engfracmech.2022.108878)\]
12. Sheveleva, A.; Loboda, V.; Lapusta, Y. A conductive crack and a remote electrode at the interface between two piezoelectric materials. Appl. Math. Model. **2020**, 87, 287–299. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=A+conductive+crack+and+a+remote+electrode+at+the+interface+between+two+piezoelectric+materials&author=Sheveleva,+A.&author=Loboda,+V.&author=Lapusta,+Y.&publication_year=2020&journal=Appl.+Math.+Model.&volume=87&pages=287%E2%80%93299&doi=10.1016/j.apm.2020.06.003)\] \[[CrossRef](https://doi.org/10.1016/j.apm.2020.06.003)\]
13. Sun, W.; Qu, W.; Gu, Y. Arbitrary-order GFDM for crack analysis of homogeneous and bi-material piezoelectric problems. Theor. Appl. Fract. Mech. **2025**, 139, 105107. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Arbitrary-order+GFDM+for+crack+analysis+of+homogeneous+and+bi-material+piezoelectric+problems&author=Sun,+W.&author=Qu,+W.&author=Gu,+Y.&publication_year=2025&journal=Theor.+Appl.+Fract.+Mech.&volume=139&pages=105107&doi=10.1016/j.tafmec.2025.105107)\] \[[CrossRef](https://doi.org/10.1016/j.tafmec.2025.105107)\]
14. Singh, R. Inspection of dynamic fracture behavior of multiple interfacial cracks emanating from circular holes in functionally graded piezoelectric bi-materials. Appl. Math. Model. **2025**, 143, 116041. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Inspection+of+dynamic+fracture+behavior+of+multiple+interfacial+cracks+emanating+from+circular+holes+in+functionally+graded+piezoelectric+bi-materials&author=Singh,+R.&publication_year=2025&journal=Appl.+Math.+Model.&volume=143&pages=116041&doi=10.1016/j.apm.2025.116041)\] \[[CrossRef](https://doi.org/10.1016/j.apm.2025.116041)\]
15. Yan, Z.; Wen, C.; Feng, W.; Zhang, C. Dynamic growth properties of an internal crack in piezoelectric plates based on the improved mechanical energy release rate criterion. Theor. Appl. Fract. Mech. **2025**, 140, 105131. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+growth+properties+of+an+internal+crack+in+piezoelectric+plates+based+on+the+improved+mechanical+energy+release+rate+criterion&author=Yan,+Z.&author=Wen,+C.&author=Feng,+W.&author=Zhang,+C.&publication_year=2025&journal=Theor.+Appl.+Fract.+Mech.&volume=140&pages=105131&doi=10.1016/j.tafmec.2025.105131)\] \[[CrossRef](https://doi.org/10.1016/j.tafmec.2025.105131)\]
16. Zhu, S.; Yu, H.; Wang, Z. Interfacial dynamic impermeable crack analysis in dissimilar piezoelectric materials by a new interaction integral. Compos. Struct. **2025**, 352, 118668. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Interfacial+dynamic+impermeable+crack+analysis+in+dissimilar+piezoelectric+materials+by+a+new+interaction+integral&author=Zhu,+S.&author=Yu,+H.&author=Wang,+Z.&publication_year=2025&journal=Compos.+Struct.&volume=352&pages=118668&doi=10.1016/j.compstruct.2024.118668)\] \[[CrossRef](https://doi.org/10.1016/j.compstruct.2024.118668)\]
17. Liu, H.; Zhu, S. Dynamic analysis of two collinear permeable Mode-I cracks in piezoelectric materials based on non-local piezoelectricity theory. Multidiscip. Model. Mater. Struct. **2019**, 15, 1274–1293. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+analysis+of+two+collinear+permeable+Mode-I+cracks+in+piezoelectric+materials+based+on+non-local+piezoelectricity+theory&author=Liu,+H.&author=Zhu,+S.&publication_year=2019&journal=Multidiscip.+Model.+Mater.+Struct.&volume=15&pages=1274%E2%80%931293&doi=10.1108/MMMS-12-2018-0215)\] \[[CrossRef](https://doi.org/10.1108/MMMS-12-2018-0215)\]
18. Liu, H.-T. Dynamic non-local theory solution to a permeable mode-I crack in a piezoelectric medium. Eng. Fract. Mech. **2017**, 179, 43–58. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+non-local+theory+solution+to+a+permeable+mode-I+crack+in+a+piezoelectric+medium&author=Liu,+H.-T.&publication_year=2017&journal=Eng.+Fract.+Mech.&volume=179&pages=43%E2%80%9358&doi=10.1016/j.engfracmech.2017.04.023)\] \[[CrossRef](https://doi.org/10.1016/j.engfracmech.2017.04.023)\]
19. Liu, H.-T.; Wu, J.-G.; Li, T.-J. Dynamic analytical solution of a limited-permeable mode-I crack in piezoelectric materials based on the non-local theory. Wave Motion **2019**, 90, 82–98. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+analytical+solution+of+a+limited-permeable+mode-I+crack+in+piezoelectric+materials+based+on+the+non-local+theory&author=Liu,+H.-T.&author=Wu,+J.-G.&author=Li,+T.-J.&publication_year=2019&journal=Wave+Motion&volume=90&pages=82%E2%80%9398&doi=10.1016/j.wavemoti.2019.05.003)\] \[[CrossRef](https://doi.org/10.1016/j.wavemoti.2019.05.003)\]
20. Afshar, H.; Bagheri, R. Several embedded cracks in a functionally graded piezoelectric strip under dynamic loading. Comput. Math. Appl. **2018**, 76, 534–550. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Several+embedded+cracks+in+a+functionally+graded+piezoelectric+strip+under+dynamic+loading&author=Afshar,+H.&author=Bagheri,+R.&publication_year=2018&journal=Comput.+Math.+Appl.&volume=76&pages=534%E2%80%93550&doi=10.1016/j.camwa.2018.04.035)\] \[[CrossRef](https://doi.org/10.1016/j.camwa.2018.04.035)\]
21. Wünsche, M.; Sladek, J.; Sladek, V.; Zhang, C.; García-Sánchez, F.; Sáez, A. Dynamic crack analysis in piezoelectric solids under time-harmonic loadings with a symmetric Galerkin boundary element method. Eng. Anal. Bound. Elem. **2017**, 84, 141–153. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+crack+analysis+in+piezoelectric+solids+under+time-harmonic+loadings+with+a+symmetric+Galerkin+boundary+element+method&author=W%C3%BCnsche,+M.&author=Sladek,+J.&author=Sladek,+V.&author=Zhang,+C.&author=Garc%C3%ADa-S%C3%A1nchez,+F.&author=S%C3%A1ez,+A.&publication_year=2017&journal=Eng.+Anal.+Bound.+Elem.&volume=84&pages=141%E2%80%93153&doi=10.1016/j.enganabound.2017.08.013)\] \[[CrossRef](https://doi.org/10.1016/j.enganabound.2017.08.013)\]
22. Yu, T.; Bui, T.Q.; Liu, P.; Zhang, C.; Hirose, S. Interfacial dynamic impermeable cracks analysis in dissimilar piezoelectric materials under coupled electromechanical loading with the extended finite element method. Int. J. Solids Struct. **2015**, 67–68, 205–218. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Interfacial+dynamic+impermeable+cracks+analysis+in+dissimilar+piezoelectric+materials+under+coupled+electromechanical+loading+with+the+extended+finite+element+method&author=Yu,+T.&author=Bui,+T.Q.&author=Liu,+P.&author=Zhang,+C.&author=Hirose,+S.&publication_year=2015&journal=Int.+J.+Solids+Struct.&volume=67%E2%80%9368&pages=205%E2%80%93218&doi=10.1016/j.ijsolstr.2015.03.037)\] \[[CrossRef](https://doi.org/10.1016/j.ijsolstr.2015.03.037)\]
23. Pak, Y.E. Crack extension force in a piezoelectric material. J. Appl. Mech. **1990**, 57, 647–653. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Crack+extension+force+in+a+piezoelectric+material&author=Pak,+Y.E.&publication_year=1990&journal=J.+Appl.+Mech.&volume=57&pages=647%E2%80%93653&doi=10.1115/1.2897071)\] \[[CrossRef](https://doi.org/10.1115/1.2897071)\]
24. Wang, X.-Y.; Yu, S.-W. Transient response of a crack in piezoelectric strip subjected to the mechanical and electrical impacts:mode-III problem. Int. J. Solids Struct. **2000**, 37, 5795–5808. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Transient+response+of+a+crack+in+piezoelectric+strip+subjected+to+the+mechanical+and+electrical+impacts:mode-III+problem&author=Wang,+X.-Y.&author=Yu,+S.-W.&publication_year=2000&journal=Int.+J.+Solids+Struct.&volume=37&pages=5795%E2%80%935808&doi=10.1016/S0020-7683\(99\)00268-1)\] \[[CrossRef](https://doi.org/10.1016/S0020-7683\(99\)00268-1)\]
25. Wang, B.-L.; Han, J.-C.; Du, S.-Y. Dynamic response for non-homogeneous piezoelectric medium with multiple cracks. Eng. Fract. Mech. **1998**, 61, 607–617. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+response+for+non-homogeneous+piezoelectric+medium+with+multiple+cracks&author=Wang,+B.-L.&author=Han,+J.-C.&author=Du,+S.-Y.&publication_year=1998&journal=Eng.+Fract.+Mech.&volume=61&pages=607%E2%80%93617&doi=10.1016/S0013-7944\(98\)00090-3)\] \[[CrossRef](https://doi.org/10.1016/S0013-7944\(98\)00090-3)\]
![Materials 19 00964 g001]()
**Figure 1.** Structure of the piezoelectric bimaterials model with cracks (The green part in the picture indicates the direction of the crack movement).
**Figure 1.** Structure of the piezoelectric bimaterials model with cracks (The green part in the picture indicates the direction of the crack movement).
![Materials 19 00964 g001]()
![Materials 19 00964 g002]()
**Figure 2.** The variation of
K
with respect to
d
1
/
a
((**a**):
ς
c
\=
2
; (**b**):
ς
c
\=
1\.5
).
**Figure 2.** The variation of
K
with respect to
d
1
/
a
((**a**):
ς
c
\=
2
; (**b**):
ς
c
\=
1\.5
).
![Materials 19 00964 g002]()
![Materials 19 00964 g003]()
**Figure 3.** (**a**) The 2D plot that
K
varies with
c
p
t
/
a
; (**b**) The three-dimensional graph that
K
varies with
c
p
t
/
a
.
**Figure 3.** (**a**) The 2D plot that
K
varies with
c
p
t
/
a
; (**b**) The three-dimensional graph that
K
varies with
c
p
t
/
a
.
![Materials 19 00964 g003]()
![Materials 19 00964 g004]()
**Figure 4.** The variation in the normalized energy release rate under far-field stress loading.
**Figure 4.** The variation in the normalized energy release rate under far-field stress loading.
![Materials 19 00964 g004]()
![Materials 19 00964 g005]()
**Figure 5.** The variation in the normalized intensity factor under different
λ
values.
**Figure 5.** The variation in the normalized intensity factor under different
λ
values.
![Materials 19 00964 g005]()
![Materials 19 00964 g006]()
**Figure 6.** Comparison chart of real values and predicted values.
**Figure 6.** Comparison chart of real values and predicted values.
![Materials 19 00964 g006]()
![Materials 19 00964 g007]()
**Figure 7.** The three-dimensional graph and contour map of the interaction between A and B in terms of the influence on the stress intensity factor.
**Figure 7.** The three-dimensional graph and contour map of the interaction between A and B in terms of the influence on the stress intensity factor.
![Materials 19 00964 g007]()
![Materials 19 00964 g008]()
**Figure 8.** The three-dimensional graph and contour map of the interaction between A and C in terms of the influence on the stress-intensity factor.
**Figure 8.** The three-dimensional graph and contour map of the interaction between A and C in terms of the influence on the stress-intensity factor.
![Materials 19 00964 g008]()
![Materials 19 00964 g009]()
**Figure 9.** The three-dimensional graph and contour map of the interaction between B and C in terms of the influence on the stress intensity factor.
**Figure 9.** The three-dimensional graph and contour map of the interaction between B and C in terms of the influence on the stress intensity factor.
![Materials 19 00964 g009]()
![Materials 19 00964 g010]()
**Figure 10.** The response values of A, B and C when the stress-intensity factor reaches its maximum value (The red and grey dots represent the corresponding data values).
**Figure 10.** The response values of A, B and C when the stress-intensity factor reaches its maximum value (The red and grey dots represent the corresponding data values).
![Materials 19 00964 g010]()
![]()
**Table 1.** Variance analysis of predictive models.
**Table 1.** Variance analysis of predictive models.
| Standard Error | Correlation Coefficient | Correction Coefficient | F | p | Significant |
|---|---|---|---|---|---|
| 0\.0285 | 0\.9975 | 0\.9886 | 706\.4 | \<0.0001 | significant |
![]()
**Table 2.** Variance analysis of predictive model of response value.
**Table 2.** Variance analysis of predictive model of response value.
| Source | Sum of Square | Df | Mean Square | F-Value | p\-Value | Significant |
|---|---|---|---|---|---|---|
| A-A | 0\.1013 | 1 | 0\.1013 | 125\.03 | \<0.0001 | significant |
| B-B | 3\.42 | 1 | 3\.42 | 4214\.47 | \<0.0001 | significant |
| C-C | 0\.0092 | 1 | 0\.0092 | 11\.32 | 0\.0120 | |
| AB | 0\.0399 | 1 | 0\.0399 | 49\.29 | 0\.0002 | |
| AC | 0\.0025 | 1 | 0\.0025 | 3\.09 | 0\.1224 | |
| BC | 0\.0012 | 1 | 0\.0012 | 1\.47 | 0\.2642 | |
| A2 | 0\.2894 | 1 | 0\.2894 | 357\.08 | \<0.0001 | significant |
| B2 | 0\.8693 | 1 | 0\.8693 | 1072\.74 | \<0.0001 | significant |
| C2 | 0\.2735 | 1 | 0\.2735 | 337\.47 | \<0.0001 | significant |
| Residual | 0\.0057 | 7 | 0\.0008 | | | |
| Lack of Fit | 0\.0034 | 3 | 0\.0011 | 2\.06 | 0\.2480 | not significant |
| Pure Error | 0\.0022 | 4 | 0\.0006 | | | |
| Cor Total | 5\.16 | 16 | | | | |
| | |
|---|---|
| | **Disclaimer/Publisher’s Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2026 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the [Creative Commons Attribution (CC BY) license](https://creativecommons.org/licenses/by/4.0/).
## Share and Cite
**MDPI and ACS Style**
Zhang, Y.; Li, J.; Li, X.; Ma, J. Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading. *Materials* **2026**, *19*, 964. https://doi.org/10.3390/ma19050964
**AMA Style**
Zhang Y, Li J, Li X, Ma J. Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading. *Materials*. 2026; 19(5):964. https://doi.org/10.3390/ma19050964
**Chicago/Turabian Style**
Zhang, Yani, Junlin Li, Xiangyu Li, and Junye Ma. 2026. "Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading" *Materials* 19, no. 5: 964. https://doi.org/10.3390/ma19050964
**APA Style**
Zhang, Y., Li, J., Li, X., & Ma, J. (2026). Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading. *Materials*, *19*(5), 964. https://doi.org/10.3390/ma19050964
Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details [here](https://www.mdpi.com/about/announcements/784).
## Article Metrics
No
No
### Article Access Statistics
For more information on the journal statistics, click [here](https://www.mdpi.com/journal/materials/stats).
Multiple requests from the same IP address are counted as one view.
[Zoom]() \| [Orient]() \| [As Lines]() \| [As Sticks]() \| [As Cartoon]() \| [As Surface]() \| [Previous Scene]() \| [Next Scene]()
## Cite
Export citation file: [BibTeX]() \| [EndNote]() \| [RIS]()
**MDPI and ACS Style**
Zhang, Y.; Li, J.; Li, X.; Ma, J. Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading. *Materials* **2026**, *19*, 964. https://doi.org/10.3390/ma19050964
**AMA Style**
Zhang Y, Li J, Li X, Ma J. Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading. *Materials*. 2026; 19(5):964. https://doi.org/10.3390/ma19050964
**Chicago/Turabian Style**
Zhang, Yani, Junlin Li, Xiangyu Li, and Junye Ma. 2026. "Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading" *Materials* 19, no. 5: 964. https://doi.org/10.3390/ma19050964
**APA Style**
Zhang, Y., Li, J., Li, X., & Ma, J. (2026). Study on Interfacial Crack of Piezoelectric Bimaterials Under Dynamic Loading. *Materials*, *19*(5), 964. https://doi.org/10.3390/ma19050964
Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details [here](https://www.mdpi.com/about/announcements/784).
[*clear*]()
*[Materials](https://www.mdpi.com/journal/materials)*, EISSN 1996-1944, Published by MDPI
[RSS](https://www.mdpi.com/rss/journal/materials) [Content Alert](https://www.mdpi.com/journal/materials/toc-alert)
### Further Information
[Article Processing Charges](https://www.mdpi.com/apc) [Pay an Invoice](https://www.mdpi.com/about/payment) [Open Access Policy](https://www.mdpi.com/openaccess) [Contact MDPI](https://www.mdpi.com/about/contact) [Jobs at MDPI](https://careers.mdpi.com/)
### Guidelines
[For Authors](https://www.mdpi.com/authors) [For Reviewers](https://www.mdpi.com/reviewers) [For Editors](https://www.mdpi.com/editors) [For Librarians](https://www.mdpi.com/librarians) [For Publishers](https://www.mdpi.com/publishing_services) [For Societies](https://www.mdpi.com/societies) [For Conference Organizers](https://www.mdpi.com/conference_organizers)
### MDPI Initiatives
[Sciforum](https://sciforum.net/) [MDPI Books](https://www.mdpi.com/books) [Preprints.org](https://www.preprints.org/) [Scilit](https://www.scilit.com/) [SciProfiles](https://sciprofiles.com/?utm_source=mpdi.com&utm_medium=bottom_menu&utm_campaign=initiative) [Encyclopedia](https://encyclopedia.pub/) [JAMS](https://jams.pub/) [Proceedings Series](https://www.mdpi.com/about/proceedings)
### Follow MDPI
[LinkedIn](https://www.linkedin.com/company/mdpi) [Facebook](https://www.facebook.com/MDPIOpenAccessPublishing) [X](https://x.com/MDPIOpenAccess)
© 1996-2026 MDPI (Basel, Switzerland) unless otherwise stated
[Disclaimer]()
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
[Terms and Conditions](https://www.mdpi.com/about/terms-and-conditions) [Privacy Policy](https://www.mdpi.com/about/privacy)
We use cookies on our website to ensure you get the best experience.
Read more about our cookies [here](https://www.mdpi.com/about/privacy).
[Accept](https://www.mdpi.com/accept_cookies)
## Share Link
[Copy]()
[*clear*]()
## Share
https://www.mdpi.com/3752464
[*clear*]()
[Back to TopTop](https://www.mdpi.com/1996-1944/19/5/964) |
| Readable Markdown | ## 1\. Introduction
Piezoelectric materials have been widely used in advanced equipment, instruments and materials due to their electromechanical coupling properties. However, during manufacturing or usage, defects such as cracks are inevitable, which may lead to structural damage or device failure. Therefore, the crack problem of piezoelectric materials has attracted considerable attention \[[1](https://www.mdpi.com/1996-1944/19/5/964#B1-materials-19-00964),[2](https://www.mdpi.com/1996-1944/19/5/964#B2-materials-19-00964),[3](https://www.mdpi.com/1996-1944/19/5/964#B3-materials-19-00964),[4](https://www.mdpi.com/1996-1944/19/5/964#B4-materials-19-00964),[5](https://www.mdpi.com/1996-1944/19/5/964#B5-materials-19-00964)\]. In recent years, piezoelectric bimaterials have been widely applied in devices. The sudden change in materials on both sides of the interface may cause stress concentration, which is one of the important reasons for material failure \[[6](https://www.mdpi.com/1996-1944/19/5/964#B6-materials-19-00964),[7](https://www.mdpi.com/1996-1944/19/5/964#B7-materials-19-00964),[8](https://www.mdpi.com/1996-1944/19/5/964#B8-materials-19-00964)\]. Li et al. \[[9](https://www.mdpi.com/1996-1944/19/5/964#B9-materials-19-00964)\] analyzed the SH-wave scattering by interfacial cracks in piezoelectric ceramic-polymer composites. Their core objective was to establish a dynamic response analysis framework for cracks under a steady-state wave field, and they preliminarily investigated the effects of frequency, incident angle, and material properties on crack propagation. This work provides a fundamental wave-based analysis model for subsequent dynamic fracture research; however, it does not address the transient response under impact loading or the coupling behavior in complex anisotropic material systems.
Research on the behavior of interface crack propagation has achieved fruitful results. Loboda et al. conducted studies on the conductive interface cracks between one-dimensional hexagonal piezoelectric quasicrystals and ordinary piezoelectric materials, analyzing the coupling effect between the cracks and the far-field electrodes, and clarifying the energy-evolution law of piezoelectric hetero-interface cracks, providing theoretical references for the analysis of related issues \[[10](https://www.mdpi.com/1996-1944/19/5/964#B10-materials-19-00964),[11](https://www.mdpi.com/1996-1944/19/5/964#B11-materials-19-00964),[12](https://www.mdpi.com/1996-1944/19/5/964#B12-materials-19-00964)\]. Sun et al. \[[13](https://www.mdpi.com/1996-1944/19/5/964#B13-materials-19-00964)\] proposed the arbitrary-order generalized finite difference method (GFDM), providing a new numerical solution path for homogeneous and dual-material piezoelectric crack problems. The crack propagation behavior of multiple interfaces caused by circular holes under dynamic loading was studied, and the influence law of functional gradient characteristics on the dynamic response of cracks was clarified, confirming that the gradient distribution of material parameters would significantly change the stress-electric field concentration effect at the crack tip \[[14](https://www.mdpi.com/1996-1944/19/5/964#B14-materials-19-00964)\]. Based on the improved mechanical energy release rate criterion, the dynamic expansion characteristics of cracks in piezoelectric plates were revealed \[[15](https://www.mdpi.com/1996-1944/19/5/964#B15-materials-19-00964)\]. In terms of innovative analysis methods, a new interaction integral method was proposed, achieving precise analysis of dynamic impermeable interface cracks in heterogeneous piezoelectric materials, effectively improving the accuracy of field quantity calculation at the crack tip \[[16](https://www.mdpi.com/1996-1944/19/5/964#B16-materials-19-00964)\]. Liu et al. \[[17](https://www.mdpi.com/1996-1944/19/5/964#B17-materials-19-00964),[18](https://www.mdpi.com/1996-1944/19/5/964#B18-materials-19-00964),[19](https://www.mdpi.com/1996-1944/19/5/964#B19-materials-19-00964)\] based on the non-local piezoelectric theory, systematically studied the dynamic response of I-type penetrating cracks in piezoelectric materials, derived the dynamic analytical solutions for finite-penetrating and complete-penetrating cracks, and improved the theoretical system of dynamic fracture in piezoelectric materials. Afshar et al. \[[20](https://www.mdpi.com/1996-1944/19/5/964#B20-materials-19-00964)\] studied the dynamic fracture behavior of multiple intruded cracks in functional gradient piezoelectric strips \[[20](https://www.mdpi.com/1996-1944/19/5/964#B20-materials-19-00964)\]. Wünsche et al. \[[21](https://www.mdpi.com/1996-1944/19/5/964#B21-materials-19-00964)\] used the symmetrical Galerkin boundary element method to establish an analysis model for dynamic cracks in piezoelectric bodies under harmonic loading, providing an efficient solution method for piezoelectric fracture problems under harmonic loading. Yu et al. \[[22](https://www.mdpi.com/1996-1944/19/5/964#B22-materials-19-00964)\] based on the extended finite element method (XFEM), established a numerical solution framework for the dynamic crack of heterogeneous piezoelectric materials under force–electric coupling loading, providing an efficient calculation tool for interface fracture problems in complex configurations.
Through in-depth analysis of the existing research, we can find that the electromechanical coupling constitutive relationship of orthotropic anisotropic materials is more complex, and the direction dependence of their physical parameters leads to the inability to directly transplant the existing theoretical models for isotropic or functionally graded materials; at the same time, although some studies involve dynamic loads, they are mainly concentrated on harmonic loads or do not fully capture the transient response characteristics under impact loads; and the existing research mostly focuses on specific types of conductive cracks, non-conductive cracks, etc., and mainly focuses on the analysis of the influence of a single factor, without deeply exploring the coupling effect and quantitative influence laws of multiple parameters such as elastic parameters, crack length, and impact time, which is difficult to meet the systematic research requirements of this article on the dynamic expansion mechanism of interface cracks in orthotropic anisotropic piezoelectric dual materials.
In summary, although the existing research has laid a foundation for the crack problem of dual-layer materials, it has obvious limitations in adapting to the specific research object, load conditions, and analysis goals of this article.
The above research has constructed a research framework for the crack problem of piezoelectric materials from aspects such as material system, crack form, analysis method, and load type. However, the theory and methods still have obvious limitations in applying to the orthogonal anisotropic piezoelectric dual-material interface crack problem under dynamic loads. On one hand, the existing research on dual-material interface cracks mostly focuses on static/sinusoidal loads, isotropic/functional gradient piezoelectric media, or do not consider the orthogonal anisotropy characteristics of the materials. The electromechanical coupling constitutive relationship of orthogonal anisotropic piezoelectric dual-material is more complex, and the transient response laws under dynamic impact loads are significantly different from those of conventional media. The existing static theories and analysis methods under sinusoidal loads cannot be directly transplanted. On the other hand, most dynamic-crack research focuses on internal cracks of a single piezoelectric material or does not deeply explore the coupling effects of multiple parameters and quantitative influence rules. There is no complete analytical–numerical coupled solution system for the dynamic-fracture problem of orthogonal anisotropic piezoelectric dual-material interfaces, and the research on the quantitative influence of key parameters such as elastic parameters, crack length, and impact time and the interaction between parameters is still scarce. In addition, most existing studies adopt a single analytical or numerical method for crack analysis. The research of integrating the response-surface method into dynamic fracture analysis to achieve rapid prediction of multiple parameters and identification of optimal working conditions is still in the exploratory stage.
During manufacturing and usage, cracks may occur at the interface due to manufacturing level and other factors. Meanwhile, under dynamic loads, the material thickness, elastic constants, etc., of the piezoelectric bimaterials will affect the dynamic propagation of interface cracks, which may lead to device failure and losses. However, there are few articles on this topic. Moreover, the dynamic propagation of interface cracks in piezoelectric bimaterials has limitations in software simulation and experiments. Therefore, precise quantitative analysis of the dynamic fracture problem of interface cracks in piezoelectric bimaterials is of great significance for structural design and safety.
The motivation for this study stems from the urgent need to address the dynamic fracture behavior of interface cracks in piezoelectric bimaterials. This is a topic that has received limited attention in existing literature, especially concerning the complex scenario involving orthogonally anisotropic materials under dynamic loads. Piezoelectric materials play a critical role in advanced equipment, instruments, and sensors due to their unique electromechanical coupling properties. However, their brittle nature and inevitable defects, such as cracks, during manufacturing and use can lead to structural damage or even device failure. While research on static fracture problems has a certain foundation, many piezoelectric components in practical applications are often subjected to dynamic loads, making the study of dynamic fracture behavior particularly important and challenging. Currently, precise quantitative analysis of dynamic interface cracks in orthogonally anisotropic piezoelectric bimaterials, whether in theoretical analysis, numerical simulation, or experimental verification, faces significant deficiencies.
Our primary contribution lies in the innovative application and integration of advanced analytical and computational methods to deeply explore this complex problem. Specifically, we employed the powerful mathematical tools of Laplace transformation and Fourier transformation, successfully converting the governing equations that describe the dynamic expansion of orthogonally anisotropic interface cracks into a complex singular integral equation. Subsequently, we applied a high-precision Chebyshev point method for numerical discretization of this integral equation, combined with efficient Laplace inversion techniques, to accurately solve for the stress-intensity factors and electric-displacement-intensity factors at the crack tip. This combined analytical and numerical approach provides a solid foundation for understanding dynamic-fracture mechanisms. Furthermore, to further enhance analysis efficiency and predictive capability, we proactively integrated the Response Surface Methodology (RSM) and specifically developed and applied a high-precision Random Forest (Random Forest) machine learning model (with a determination coefficient R2 as high as 0.9886, demonstrating the model’s excellent goodness of fit and predictive accuracy) for effective analysis and prediction of complex multi-parameter influences.
This multi-level comprehensive method offers significant added value. It not only precisely reveals the quantitative influence laws of key factors such as crack length, elastic parameters (e.g., modulus), and impact time on dynamic fracture behavior—for example, revealing the positive correlation between crack length and the dimensionless function, as well as the inverse relationship between elastic parameters and the dimensionless dynamic stress intensity factor under certain conditions—but also captures the dynamic evolution process induced by impact time, where the dynamic crack response first increases, reaches a peak, and then oscillates and converges to the static value. Moreover, it ultimately enables the identification of optimal conditions affecting stress intensity. This in-depth and quantitative understanding greatly enhances the prediction accuracy of dynamic crack propagation mechanisms, providing a key theoretical basis for evaluating and optimizing the dynamic reliability of piezoelectric devices, thereby significantly improving the operational reliability and structural safety of piezoelectric equipment in practical engineering applications. By directly addressing the limitations of current software simulations in handling dynamic fracture problems of complex anisotropic materials, as well as the challenges of experimental research in terms of cost, scale, and dynamic loading control, this study aims to fill a critical gap in the field and provide important scientific support and technical guidance for the future design of more robust, safer, and more reliable piezoelectric smart structures.
## 2\. Establishment of the Model
As shown in [Figure 1](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f001), the thicknesses of the two materials are d 1 and d 2, respectively. The crack is from a 1 to b 1. With the interface direction as the X-axis and the thickness direction as the Y-axis. Let us assume that the polarization direction of the piezoelectric dual materials is the Z-axis, the structure is orthotropic, and the symmetry axes of the materials are perfectly aligned at the interface.
The governing equation for the crack in the piezoelectric bimaterial orthotropic interface is:
c 55 ( k ) ∂ 2 w ∂ X 2 \+ c 44 ( k ) ∂ 2 w ∂ Y 2 \+ e 15 ( k ) ∂ 2 ϕ ∂ X 2 \+ e 24 ( k ) ∂ 2 ϕ ∂ Y 2 \= ρ ∂ 2 w ∂ t 2
(1)
e 15 ( k ) ∂ 2 w ∂ X 2 \+ e 24 ( k ) ∂ 2 w ∂ Y 2 − ε 11 ( k ) ∂ 2 ϕ ∂ X 2 − ε 22 ( k ) ∂ 2 ϕ ∂ Y 2 \= 0
(2)
In the above equations, w \= w ( x , y , t ) denotes the mechanical displacements, φ ( x , y , t ) represents the electric potentials, ρ is the mass density, c 55 ( k ) and c 44 ( k ) are the elastic stiffness constants, e 15 ( k ) and e 24 ( k ) stand for the piezoelectric constants, and ε 11 ( k ) and ε 22 ( k ) are the dielectric constants. Here, k \= 1 , 2 correspond to piezoelectric material 1 and piezoelectric material 2, respectively.
Introduce the transformation x \= X − v t, y \= Y, then Equations (1) and (2) becomes:
c 55 ( k ) c 44 ( k ) − v 2 c k ¯ 2 1 e 15 ( k ) c 44 ( k ) e 24 ( k ) c 44 ( k ) e 15 ( k ) c 44 ( k ) e 24 ( k ) c 44 ( k ) − ε 11 ( k ) c 44 ( k ) − ε 22 ( k ) c 44 ( k ) ∂ 2 w ∂ x 2 ∂ 2 w ∂ y 2 ∂ 2 ϕ ∂ x 2 ∂ 2 ϕ ∂ y 2 \= Ψ ( k ) 0
(3)
Perform the Laplace transform on time t and the Fourier transform on x. The solution of Equation (3) can be expressed as
\[ w \* x , y , p ϕ \* x , y , p \] \= 1 2 π ∫ − ∞ \+ ∞ \[ χ 11 ( k ) χ 12 ( k ) χ 21 ( k ) χ 22 ( k ) \] \[ A 11 ( k ) e ( 3 − 2 k ) s γ 1 ( k ) y A 12 ( k ) e ( 3 − 2 k ) s γ 2 ( k ) y \] e − i s x d s
(4)
where
χ 1 η ( k ) \= e 24 ( k ) c 44 ( k ) γ η ( k ) 2 − s 2 e 15 ( k ) c 44 ( k ) , χ 2 η ( k ) \= γ η ( k ) 2 − s 2 \[ ( c 55 ( k ) c 44 ( k ) − v 2 c k ¯ 2 ) \+ 2 v p i c k ¯ 2 s \+ p 2 c k ¯ 2 s 2 \] , η \= 1 , 2 .
The value of γ η ( k ) can be found in [Appendix A](https://www.mdpi.com/1996-1944/19/5/964#app1-materials-19-00964).
By analyzing the model, it can be concluded that the boundary conditions that the cracks must satisfy are
σ y z ( x , d 1 , t ) \= σ y z ( x , − d 2 , t ) \= 0
(5)
w ( x , 0 \+ , t ) \= w ( x , 0 − , t ) , D y ( x , 0 \+ , t ) \= D y ( x , 0 − , t )
(6)
σ y z ( x , 0 \+ , t ) \= σ y z ( x , 0 − , t ) \= − τ 0 H ( t ) D y ( x , 0 \+ , t ) \= D y ( x , 0 − , t ) \= − D 0 H ( t )
(7)
where H ( t ) \= 0 , t \< 0 1 , t \> 0.
By analyzing the model, it can be concluded that the boundary conditions that the cracks must satisfy are
σ y z \* ( x , d 1 , p ) \= σ y z \* ( x , − d 2 , p ) \= 0
(8)
w \* ( x , 0 \+ , p ) \= w \* ( x , 0 − , p ) , D y \* ( x , 0 \+ , p ) \= D y \* ( x , 0 − , p )
(9)
σ y z \* ( x , 0 \+ , p ) \= σ y z \* ( x , 0 − , p ) \= − τ 0 p , D y \* ( x , 0 \+ , p ) \= D y \* ( x , 0 − , p ) \= − D 0 p
(10)
The electric displacement corresponding to the stress in Equation (4) is
σ y z \* D y \* \= 1 2 π ∫ − ∞ \+ ∞ ( c 44 ( k ) e 24 ( k ) e 15 ( k ) − ε 22 ( k ) χ 11 ( k ) χ 12 ( k ) χ 21 ( k ) χ 22 ( k ) γ 1 ( k ) e ( 3 − 2 k ) s γ 1 ( k ) y A 11 ( k ) γ 2 ( k ) e ( 3 − 2 k ) s γ 2 ( k ) y A 12 ( k ) ( 3 − 2 k ) s ) e − i s x d s
(11)
## 3\. Solution of the Model
In the process of crack propagation, the expansion of the crack is affected by a variety of complex factors. To more accurately describe the crack propagation behavior, introducing a density function is an effective method. Therefore, the following density function is introduced:
h w ( x ) \= d d x \[ w ( x , 0 \+ , p ) − w ( x , 0 − , p ) \]
(12)
h ϕ ( x ) \= d d x \[ ϕ ( x , 0 \+ , p ) − ϕ ( x , 0 − , p ) \]
(13)
Since displacement and potential are continuous, we can obtain:
∫ a 1 b 1 h w ( x ) d x \= 0 , x ∉ ( a 1 , b 1 )
(14)
∫ a 1 b 1 h ϕ ( x ) d x \= 0 , x ∉ ( a 1 , b 1 )
(15)
Based on the boundary conditions (8–10) and Equations (12) and (13), the following conclusion can be obtained:
χ 11 ( 1 ) χ 12 ( 1 ) χ 21 ( 1 ) χ 22 ( 1 ) A 11 B 11 − χ 11 ( 2 ) χ 12 ( 2 ) χ 21 ( 2 ) χ 22 ( 2 ) A 12 B 12 \= i s ∫ a 1 b 1 h w ( s ) h ϕ ( s ) e i s f d f \= H 1 H 2
(16)
A 11 ( 1 ) A 12 ( 1 ) A 11 ( 2 ) A 12 ( 2 ) \= 1 A \* A 1 B 1 A 2 A 2 A 1 \* B 1 \* A 2 \* B 2 \* H 1 H 2
(17)
The values of A 1 , B 1 , A 2 , B 2 , A 1 \* , B 1 \* , A 2 \* , B 2 \* , A \* can be found in [Appendix A](https://www.mdpi.com/1996-1944/19/5/964#app1-materials-19-00964).
Substituting Equation (17) into Equation (11) yields
σ y z \* ( x , 0 \+ , p ) D y \* ( x , 0 \+ , p ) \= lim y → 0 \+ 1 2 π ∫ − ∞ \+ ∞ s e s γ y A \* C 1 C 2 C 3 C 4 H 1 H 2 e − i s x d s
(18)
where
C 1 \= ∑ i \= 1 2 ( c 44 ( 1 ) χ 1 i ( 1 ) \+ e 24 ( 1 ) χ 2 i ( 1 ) ) A i γ i ( 1 ) ,
C 2 \= ∑ i \= 1 2 ( c 44 ( 1 ) χ 1 i ( 1 ) \+ e 24 ( 1 ) χ 2 i ( 1 ) ) B i γ i ( 1 ) ,
C 3 \= ∑ i \= 1 2 ( e 15 ( 1 ) χ 1 i ( 1 ) − ε 22 ( 1 ) χ 2 i ( 1 ) ) A i \* γ i ( 1 ) ,
C 4 \= ∑ i \= 1 2 ( e 15 ( 1 ) χ 1 i ( 1 ) − ε 22 ( 1 ) χ 2 i ( 1 ) ) B i \* γ i ( 1 ) ,
e s γ y \= ∑ i \= 1 2 e s γ i ( 1 ) y
On the crack surface, there is
∫ a 1 b 1 ( C 1 \* C 2 \* C 3 \* C 4 \* 1 B \* ( f − x ) \+ k 1 ( x , f ) k 2 ( x , f ) k 3 ( x , f ) k 4 ( x , f ) ) h w ( f ) h ϕ ( f ) d f \= − π p σ 0 D 0
(19)
where
k 1 ( x , f ) \= ∫ 0 ∞ ( C 1 A \* − C 1 \* B \* ) sin ( s ( f − x ) ) d s ,
k 2 ( x , f ) \= ∫ 0 ∞ ( C 2 A \* − C 2 \* B \* ) sin ( s ( f − x ) ) d s ,
k 3 ( x , f ) \= ∫ 0 ∞ ( C 3 A \* − C 3 \* B \* ) sin ( s ( f − x ) ) d s ,
k 4 ( x , f ) \= ∫ 0 ∞ ( C 4 A \* − C 4 \* B \* ) sin ( s ( f − x ) ) d s .
The value of C 1 \* , C 2 \* , C 3 \* , C 4 \* , B \* is referred to in [Appendix A](https://www.mdpi.com/1996-1944/19/5/964#app1-materials-19-00964).
To simplify the calculation, we will transform the normalized quantity into the standard form and introduce the dislocation density function
h w ( x ) \= q 1 ( x ) , h ϕ ( x ) \= q 2 ( x ) , k l ( x , f ) \= P l ( r , c ) , ( l \= 1 , 2 , 3 , 4 ) ,
m o \= b 1 − a 1 2 , n o \= b 1 \+ a 1 2 , f \= x − n 0 m 0 , c \= f − n 0 m 0 ,
(20)
q 1 ( c ) \= h 1 ( c ) 1 − c 2 π σ 0 p , q 2 ( c ) \= h 2 ( c ) 1 − c 2 π D 0 p .
According to Equation (20) and the Chebyshev placement method, Equation (19) can be transformed into
1 N ∑ Q \= 0 N Φ Q ( C 1 \* C 2 \* C 3 \* C 4 \* Ξ \+ m 0 P 1 ( r u , c Q ) P 2 ( r u , c Q ) P 3 ( r u , c Q ) P 4 ( r u , c Q ) ) h 1 ( c Q ) h 2 ( c Q ) \= − 1 − 1
(21)
where Ξ \= m 0 B \* ( m 0 c Q \+ n 0 − m 0 r u − n 0 ), ∑ Q \= 0 N Φ Q h 1 ( c Q ) \= 0, ∑ Q \= 0 N Φ Q h 2 ( c Q ) \= 0, u \= 1 , 2 , … N, u \= 1 , 2 , … N, Φ 0 \= Φ N \= 1 2, Φ 1 \= … \= Φ N − 1 \= 1, r u \= cos ( 2 u − 1 ) π 2 N, c Q \= cos Q π N (N represents a node).
## 4\. Intensity Factor
According to Equation (21), it can be obtained that
lim x → a 1 − σ y z \* ( x , 0 , p ) \= σ 0 C 1 \* p B \* lim c → − 1 − h 1 ( − 1 ) c 2 − 1 \+ σ 0 C 2 \* p B \* lim c → − 1 − h 2 ( − 1 ) c 2 − 1
(22)
lim x → b 1 \+ σ y z \* ( x , 0 , p ) \= σ 0 C 1 \* p B \* lim c → 1 \+ − h 1 ( 1 ) c 2 − 1 \+ σ 0 C 2 \* p B \* lim c → 1 \+ − h 2 ( 1 ) c 2 − 1
(23)
lim x → a 1 − D y \* ( x , 0 , p ) \= D 0 C 3 \* p B \* lim c → − 1 − h 1 ( − 1 ) c 2 − 1 \+ D 0 C 4 \* p B \* lim c → − 1 − h 2 ( − 1 ) c 2 − 1
(24)
lim x → b 1 \+ D y \* ( x , 0 , p ) \= D 0 C 3 \* p B \* lim c → 1 \+ − h 1 ( 1 ) c 2 − 1 \+ D 0 C 4 \* p B \* lim c → 1 \+ − h 2 ( 1 ) c 2 − 1
(25)
Using Miller and Guy, assuming that the Laplace transform of f ( t ) is f \* ( p ), p can be given by the following discrete points
p \= δ ( β \+ 1 \+ k )
(26)
where δ \> 0, β \> − 1, k \= 0 , 1 , 2 , ⋯, Then it can be concluded that
f ( t ) \= ∑ T \= 0 N B T P T ( 0 , β ) 2 e − δ t − 1 , δ \> 0 , β \> − 1
(27)
Here, P T ( 0 , β ) is a Jacobi polynomial, and the coefficients B T are determined by the following system of equations
δ f \* ( δ ( β \+ 1 \+ k ) ) \= ∑ m \= 0 k k ( k − 1 ) ⋯ k − ( m − 1 ) ( k \+ β \+ 1 ) ( k \+ β \+ 2 ) ⋯ ( k \+ β \+ 1 \+ m ) B T
(28)
Based on the definitions of stress and electric displacement, the stress intensity factors K I I I and the electric displacement intensity factors K I V at the crack tips x \= a 1 and x \= b 1 can be obtained as
K I I I ( a 1 ) K I V ( a 1 ) \= 2 π ( a 1 − x ) x → a 1 − 0 σ y z ( x , 0 , t ) D y ( x , 0 , t )
(29)
K I I I ( b 1 ) K I V ( b 1 ) \= 2 π ( x − b 1 ) x → b 1 \+ 0 σ y z ( x , 0 , t ) D y ( x , 0 , t )
(30)
Here, a non-dimensional function is defined
K \= K I I I ( t , v ) / K I I I ( 0 , v )
(31)
Here, K I I I ( 0 , v ) represents the stress intensity factor under static load conditions.
## 5\. Far-Field Forces and Electric Impact Loading
For the crack problem in piezoelectric media, as a cutting-edge branch of multi-physics field-coupled fracture mechanics, its core difficulty lies in accurately describing the complex boundary conditions at the crack tip and its surface. These conditions involve not only mechanical stress and strain but are also closely related to the electric field and electric displacement, forming a highly nonlinear boundary-value problem.
The classical strict electro-mechanical boundary conditions require that the normal displacement components at the surface of the piezoelectric medium (as the defect interface) and within the material must be continuous in the vicinity of the crack tip and at the crack interface, and the potential or electric field may also need to satisfy specific continuity or jump conditions. This strict condition stems from the principle of physical continuity and is currently the most accurate theoretical description. However, as stated in the original text, this greatly limits the scope of obtaining analytical solutions or closed-form solutions, making it extremely difficult to deeply understand and predict the fracture behavior of piezoelectric materials under complex conditions, severely restricting the application of related theories in engineering practice.
To overcome this challenge, the engineering and academic communities have developed a practical and effective approximation method. The core idea is to simplify the electric field based on the “transparency” of the crack surface. Currently, the most widely accepted and applied approach is to classify crack-boundary conditions into two categories: impermeable cracks and permeable cracks. The impermeable crack model, also known as the D-P boundary condition, is the most commonly used simplified model. It is based on a key assumption: the crack surface is completely insulated and does not allow any free charges, so the free charge density on the crack surface is zero. According to Gauss’s law, this means that the normal electric displacement component on the crack surface must be zero. Additionally, this model typically assumes that the potential on the crack surface remains constant. This boundary condition significantly simplifies the mathematical processing as it directly eliminates the influence of the crack surface as an electric field boundary, making the crack surface a “free surface” or an “equipotential surface”.
Although from a physical perspective this might be overly idealized (completely preventing the electric field from penetrating), in many engineering applications, especially when the crack propagation direction is not consistent with the electric field direction or when the crack size is much smaller than the characteristic wavelength, the D-P condition can provide results with reasonable accuracy, becoming a common benchmark for fracture mechanics analysis (such as calculating stress-intensity factors, electric-displacement-intensity factors) and structural-safety assessment. The permeable-crack model takes the opposite extreme, allowing the electric field to pass through the crack surface. In the permeable model, it is assumed that the potentials on the upper and lower surfaces of the crack are equal (i.e., the potential is continuous), or more generally, it allows for a certain jump in the potential on both sides of the crack surface, but the key is that the normal electric displacement components on both sides of the crack surface remain continuous. This means that the crack surface is no longer completely insulated but allows current to pass through. The permeable model typically better reflects the actual distribution of the electric field at the crack, especially in strong electric-field effects or specific geometric structures. However, its mathematical processing is usually more complex, and analytical solutions are more difficult to obtain.
By introducing these approximate boundary conditions, although the accuracy of the model is sacrificed to some extent, it significantly reduces the difficulty of solving the piezoelectric crack problem, enabling the use of analytical methods, semi-analytical methods to analyze fracture behavior, calculate stress intensity factors to be possible, and assess the safety of structures under complex load coupling. These approximate models provide important theoretical tools and engineering methods for understanding and predicting the failure mechanism of piezoelectric intelligent structures (such as piezoelectric sensors, detectors, and actuators).
In this study, we investigate the crack-growth problem subjected to stress and electric displacement impact loading at infinity, denoted by τ∞H(t) and D∞H(t), respectively. The boundary conditions described in Equations (5) and (7) are revised and transformed into the following forms
σ y z ( x , d 1 , t ) \= σ y z ( x , − d 2 , t ) \= 0
(32)
w ( x , 0 \+ , t ) \= w ( x , 0 − , t ) , D y ( x , 0 \+ , t ) \= D y ( x , 0 − , t )
(33)
σ y z ( x , 0 \+ , t ) \= σ y z ( x , 0 − , t ) \= − τ ∞ H ( t ) D y ( x , 0 \+ , t ) \= D y ( x , 0 − , t ) \= D \* − D ∞ H ( t )
(34)
where D \* is the electrical displacement in the crack plane. Take the Laplace transform of Equation (34), we obtain
σ y z \* ( x , 0 \+ , p ) \= σ y z \* ( x , 0 − , p ) \= − τ ∞ p , D y \* ( x , 0 \+ , p ) \= D y \* ( x , 0 − , p ) \= − D ∞ − D \* p
(35)
### 5\.1. Impermeable Crack (D-P)
For impermeable cracks, the electric potential in the crack cavity can be reasonably neglected. This is because the dielectric constant of piezoelectric materials is usually three to four orders of magnitude higher than that of air or vacuum. Therefore, the physical quantities appearing in Equation (35) D \* meet the following conditions \[[23](https://www.mdpi.com/1996-1944/19/5/964#B23-materials-19-00964)\]:
D \* \= 0
(36)
According to the aforementioned derivation process, the stress-intensity factor represented by a 1 and b 1 can be expressed as
lim x → a 1 − σ y z \* ( x , 0 , p ) \= σ ∞ C 1 \* p B \* lim c → − 1 − h 1 ( − 1 ) c 2 − 1 \+ σ ∞ C 2 \* p B \* lim c → − 1 − h 2 ( − 1 ) c 2 − 1
(37)
lim x → b 1 \+ σ y z \* ( x , 0 , p ) \= σ ∞ C 1 \* p B \* lim c → 1 \+ − h 1 ( 1 ) c 2 − 1 \+ σ ∞ C 2 \* p B \* lim c → 1 \+ − h 2 ( 1 ) c 2 − 1
(38)
lim x → a 1 − D y \* ( x , 0 , p ) \= D ∞ C 3 \* p B \* lim c → − 1 − h 1 ( − 1 ) c 2 − 1 \+ D ∞ C 4 \* p B \* lim c → − 1 − h 2 ( − 1 ) c 2 − 1
(39)
lim x → b 1 \+ D y \* ( x , 0 , p ) \= D ∞ C 3 \* p B \* lim c → 1 \+ − h 1 ( 1 ) c 2 − 1 \+ D ∞ C 4 \* p B \* lim c → 1 \+ − h 2 ( 1 ) c 2 − 1
(40)
In terms of the definitions of stress and electric displacement, the corresponding factors at the crack fronts can be formally derived as:
K I I I ( a 1 ) K I V ( a 1 ) \= 2 π ( a 1 − x ) x → a 1 − 0 σ y z ( x , 0 , t ) D y ( x , 0 , t )
(41)
K I I I ( b 1 ) K I V ( b 1 ) \= 2 π ( x − b 1 ) x → b 1 \+ 0 σ y z ( x , 0 , t ) D y ( x , 0 , t )
(42)
According to Equations (41) and (42) and combined with the virtual-crack-closure technique, the corresponding energy release rate G can be obtained \[[24](https://www.mdpi.com/1996-1944/19/5/964#B24-materials-19-00964)\].
### 5\.2. Permeable Crack
For the penetrating cracks, when the material in the text degenerates into an isotropic material, the coupling degree of the mechanical and electrical loads can be expressed as
c 55 ( k ) ( ∂ 2 w ∂ X 2 \+ ∂ 2 w ∂ Y 2 ) \+ e 15 ( k ) ( ∂ 2 ϕ ∂ X 2 \+ ∂ 2 ϕ ∂ Y 2 ) \= ρ ∂ 2 w ∂ t 2
(43)
e 15 ( k ) ( ∂ 2 w ∂ X 2 \+ ∂ 2 w ∂ Y 2 ) − ε 11 ( k ) ( ∂ 2 ϕ ∂ X 2 \+ ∂ 2 ϕ ∂ Y 2 ) \= 0
(44)
e 15 D ∞ ε 11 τ ∞ \= λ
(45)
Note: In the previous calculation, c 55 ( k ) \= c 44 ( k ), e 15 ( k ) \= e 24 ( k ) and ε 11 ( k ) \= ε 22 ( k ).
For the permeable crack model, the continuity of the electric potential at the upper and lower surfaces of the crack is the core electrical boundary condition, which means that the electric potential of the media on both sides of the crack is equal everywhere on the crack surface
h ϕ ( x ) \= 0
(46)
According to the aforementioned derivation process, the stress-intensity factor represented by a 1 and b 1 can be expressed as
lim x → a 1 − σ y z \* ( x , 0 , p ) \= σ ∞ C 1 \* p B \* lim c → − 1 − h 1 ( − 1 ) c 2 − 1 \+ σ ∞ C 2 \* p B \* lim c → − 1 − h 2 ( − 1 ) c 2 − 1
(47)
lim x → b 1 \+ σ y z \* ( x , 0 , p ) \= σ ∞ C 1 \* p B \* lim c → 1 \+ − h 1 ( 1 ) c 2 − 1 \+ σ ∞ C 2 \* p B \* lim c → 1 \+ − h 2 ( 1 ) c 2 − 1
(48)
lim x → a 1 − D y \* ( x , 0 , p ) \= λ ε 11 τ ∞ C 3 \* e 15 p B \* lim c → − 1 − h 1 ( − 1 ) c 2 − 1 \+ λ ε 11 τ ∞ C 4 \* e 15 p B \* lim c → − 1 − h 2 ( − 1 ) c 2 − 1
(49)
lim x → b 1 \+ D y \* ( x , 0 , p ) \= λ ε 11 τ ∞ C 3 \* e 15 p B \* lim c → 1 \+ − h 1 ( 1 ) c 2 − 1 \+ λ ε 11 τ ∞ C 4 \* e 15 p B \* lim c → 1 \+ − h 2 ( 1 ) c 2 − 1
(50)
## 6\. Numerical Examples
To verify the above method and analyze its application, we simulated some numerical examples. In the following examples, a represents the half-length of the crack.
### 6\.1. The Cracks Change in Accordance with the Influence of Elastic Constants
The variation in the dimensionless function with different elastic constant ratios is shown in [Figure 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f002)a,b. When the elastic constant ratio ς c \= c 44 ( 2 ) / c 44 ( 1 ) was different, the piezoelectric constant ratio and dielectric constant ratio were both 1. The crack changes under different thickness ratios of materials 1 and 2 (d 1 / d 2 \= 0\.0125 and d 1 / d 2 \= 0\.01) were analyzed respectively in the cases of elastic constant ratio ς c \= 2 ([Figure 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f002)a) and ς c \= 1\.5 ([Figure 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f002)b).
As can be seen from the figure, the non-dimensional function decreases continuously as d 1 / a increases and gradually approaches a stable value. Moreover, when the elastic modulus is large, the stable value of the non-dimensional function is also larger. At the same time, it can be observed that when the thickness of material 2 is fixed, the larger the thickness of material 1, the smaller the influence on the non-dimensional function. That is, the smaller the thickness of material 1, the greater the stress intensity generated. Therefore, in application, if we make material 1 relatively thinner, it will have higher practical application value.
### 6\.2. The Variation in Crack Under Impact Time and Crack Length
It can be clearly seen from the two-dimensional and three-dimensional curves presented in [Figure 3](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f003)a,b that the dimensionless dynamic stress intensity factor at the crack tip of the piezoelectric bimaterial interface shows obvious transient dynamic behaviors under impact load. Over time, the dimensionless function first rises rapidly, reaching a significant peak within a short period of time, then begins to gradually decline, presenting small oscillations within a certain range, and finally converges gradually to the corresponding stress intensity factor value under the static load.
This change process fully demonstrates the comprehensive influence of stress-wave propagation, reflection, superposition, and energy dissipation under impact load on the mechanical field at the crack tip and also reflects the significant differences between dynamic fracture behavior and quasi-static fracture behavior. The appearance of the peak indicates that in the initial stage of impact, the stress wave focuses strongly at the crack tip, resulting in a significant enhancement of the local mechanical and electric field coupling effect, while the oscillation convergence in the later stage indicates that the system gradually stabilizes and the dynamic effect gradually decays.
At the same time, the numerical results also show that, while keeping the thickness of material 1 unchanged, the crack length has a significant positive regulatory effect on the dimensionless function. As the crack length increases, the dimensionless dynamic stress intensity factor shows a significant increasing trend overall, which means that the expansion of the crack size will further intensify the stress concentration at the interface and increase the risk of dynamic expansion and unstable fracture of the crack.
These laws not only reveal the intrinsic relationship between geometric parameters, load forms and dynamic fracture behavior, but also provide important theoretical basis and data support for the strength design, damage identification and safety assessment of piezoelectric composite material structures under complex dynamic conditions such as impact and vibration.
### 6\.3. Impermeable and Permeable Cracks
The variation in the normalized energy-release rate under far-field stress loading is shown in [Figure 4](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f004). The dynamic nondimensional energy-release rate can be expressed as G / G r (G r is the dynamic energy-release rate for a stationary crack). [Figure 4](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f004) shows that when the far-field shear stress τ ∞ increases from zero, the normalized energy release rate G / G r always shows an upward trend. However, when τ ∞ decreases from zero (i.e., in a negative shear stress state), G / G r first decreases and then increases. This phenomenon demonstrates that the effects of positive and negative shear stresses on crack propagation have significant asymmetry: positive shear stress can continuously increase the energy release rate, thereby promoting crack propagation, while negative shear stress will reduce the energy release rate within a certain small range, manifesting as inhibition or hindrance of crack development. However, as the absolute value of negative shear stress further increases, the energy release rate will rise again. Thus, stress loading can either promote crack propagation or have interference or blocking effects, the specific manifestation depending on the direction and magnitude of the far-field stress. The above change pattern is completely consistent with the theoretical analysis results established by Pak \[[23](https://www.mdpi.com/1996-1944/19/5/964#B23-materials-19-00964)\], further verifying the rationality of the calculation model and the shear fracture theory presented in this paper.
In order to conduct a more in-depth study on the dynamic crack propagation under the coupling impact of mechanics and electricity, [Figure 5](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f005) presents the evolution of the normalized dynamic strength factor under different electromagnetic–mechanical coupling coefficients λ. It is clearly observable that as the parameter ct/(2a) related to time increases, the normalized dynamic strength factor F 0 first rises rapidly, reaching a significant peak, then gradually decreases and eventually stabilizes at a constant value. This characteristic trend, including the sharp rise, peak response, and final steady-state behavior, is highly consistent with the dynamic fracture response studied by Wang \[[25](https://www.mdpi.com/1996-1944/19/5/964#B25-materials-19-00964)\].
### 6\.4. Response-Surface Statistical Analysis
This study employed Response Surface Methodology (RSM), a statistical technique that integrates experimental design, modeling, and optimization techniques, aiming to systematically explore and quantify the influence of multiple controllable factors on a specific response variable, and ultimately determine the conditions for obtaining the optimal response. In this study, we specifically focused on the influence pattern of three key factors—Factor A (d 1 / a), Factor B (c p t / a), and Factor C (c 44 ( 2 ) / c 44 ( 1 ))—on the Stress Intensity Factor R1 (SIF) response variable.
First, to efficiently and comprehensively obtain data points for modeling, we adopted the Box–Behnken Design (BBD) experimental design method provided by Design-Expert software V22. BBD is a commonly used response surface design method, particularly suitable for situations where it is necessary to estimate the interaction and quadratic effects between factors.
Next, we used the statistical module of Design-Expert to perform a quadratic polynomial regression fit on the data. To assess the reliability and explanatory power of the established model, we conducted rigorous statistical analyses. The model-fitting goodness was measured by the coefficient of determination R2 (R-squared) and the adjusted coefficient of determination Adjusted R2. The R2 value represents the proportion of the total variation in the response variable explained by the model; the closer it is to 1, the better the model fit. Adjusted R2 considers the number of independent variables in the model and adjusts R2, providing a more realistic reflection of the model’s actual predictive ability. At the same time, we performed an Analysis of Variance (ANOVA). Analysis of Variance (ANOVA) decomposes the contribution of independent variables to the total variation to assess the impact of each factor on the response variable (Stress Intensity Factor R1). The significance of each factor on the response value is mainly achieved by comparing the variance between groups with that within groups, thereby evaluating the accuracy and reliability of the model. The significance of the model is tested by analyzing the sources of error, and the F-value and p\-value are combined for joint judgment. The larger the F value and the smaller the p value, the more significant the model is. When the p value is greater than 0.05, the model is not significant. When the p\-value is less than 0.05, the model is significant, especially when the p\-value is less than 0.0001, the model shows extremely high significance. According to the data in [Table 1](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t001), the p value is less than 0.0001, so the model is significant. The ANOVA results not only validated the statistical significance of the entire model (through F-test, typically requiring a p\-value \< 0.05) but also tested the significance of individual coefficients (linear, quadratic, interaction terms) (through t\-test, typically requiring a p\-value \< 0.05 or 0.01), helping us understand which factors and their interactions significantly affect SIF. Residual analysis was also conducted to further validate the model’s predictive ability and the reasonableness of its assumptions.
Finally, based on the established and validated quadratic response surface model, we utilized the optimization tools built into the Design-Expert software to perform optimization calculations on the combination of Factors A, B, and C. The optimization objective was set to maximize the Stress Intensity Factor (SIF) response value R1. The optimization process considered all significant factor effects (including linear, quadratic, and interaction effects) and may have set operational constraints on the factors (such as factors not exceeding their physical or process-allowable ranges). The optimization algorithm in Design-Expert searched the factor space to find the optimal combination of Factors A, B, and C that maximizes SIF under the constraint conditions. This optimal combination not only provides clear guidance for practical applications but also intuitively demonstrates the core value of RSM—namely, finding optimal operating conditions from complex factor interactions through mathematical modeling and optimization.
The results of the analysis of variance (as shown in [Table 1](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t001) and [Table 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t002)) show that the F value of the regression model is 706.4, and the p values of each model are all less than 0.05, indicating the significance of the model. Among them, the p values of A-A, B-B, A2, B2, and C2 are all less than 0.0001, indicating that the model is extremely significant and highly reliable, and the experimental design is reasonable. Therefore, the results of the analysis of variance indicate that the fitting model is accurate and statistically significant.
According to the data, a second-order response surface mathematical model was obtained through analysis. The second-order response surface regression model regarding the three variables A, B, and C, as well as R1, is shown in Equation (51).
R 1 \= − 12\.2865 \+ 0\.8313 A \+ 1\.5650 B \+ 14\.0871 C \+ 0\.1817 AB \+ 0\.1818 AC − 0\.0691 BC − 0\.8666 A 2 − 0\.4544 B 2 − 4\.07760 C 2
(51)
[Figure 6](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f006) depicts the relationship between the predicted values and the actual values. The horizontal axis represents the actual values, and the vertical axis represents the predicted values. Under ideal conditions, if the predicted value is equal to the actual value, it lies on the black straight line. Different data points are represented by squares of different colors. By observing the image, it can be found that the vast majority of data points are concentrated near the black straight line, which directly indicates that the predicted values are close to the actual values. This also shows that the established prediction model has good accuracy and reliability, can predict the relevant values relatively accurately, and can effectively reflect the true relationship between the variables.
It can be seen from [Table 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t002) that the interaction terms AC and BC are statistically insignificant. This does not mean that there are no interaction effects at the physical level, but reflects that the strengths of these interactions are weaker compared with the main effects of factors A, B, and C. During the experiments, the coupled electromechanical effects corresponding to AC and BC do exist objectively, but their influences on the response variables (e.g., interfacial electromechanical coupling strength, stress–electric coupling distribution at the crack tip, etc.) are overshadowed by the dominant main effects. For example, A (thickness of Material 1/crack length) directly determines the constraint effect of Material 1 on crack propagation, and C (ratio of elastic constants of Material 2 to Material 1) mainly affects the mechanical matching across the interface. Their independent effects on the electromechanical coupling are more significant, whereas their synergistic regulation only induces marginal changes. In addition, we have verified the residual distribution of the model and found no systematic deviation, indicating that the insignificance of AC and BC does not result from improper model specification, but truly reflects the weak intensity of these interactions.
From the perspective of physical mechanism, the coupled electromechanical effects at the interface are essentially a complex process involving multi-factor synergy. According to the variable definitions, the synergistic pathways of AC and BC are inherently constrained by the geometrical parameters and mechanical properties of the materials, limiting their interactive contributions to a small range. Specifically, for the interaction term BC, the variation in crack propagation length is directly related to the stress concentration at the crack tip, while the elastic constant ratio determines the stress-transfer efficiency across the interface. During crack propagation, the regulation of stress transfer by the elastic constant ratio is overshadowed by the stress-concentration effect at the crack tip, producing only weak local effects that cannot be manifested as statistically significant changes in the macroscopic response variables. For the interaction term AC, their synergistic effect theoretically affects the distribution of interfacial electromechanical coupling. However, when the thickness of Material 1 is sufficient to form a stable load-bearing structure, the contribution of the interaction to the response variable becomes weak. Due to its small amplitude, it does not reach the threshold of statistical significance in the macroscopic response model, thus showing statistical insignificance.
In summary, the statistical insignificance of the interaction terms AC and BC is a reasonable reflection of the coupled electromechanical effects at the interface under the constraints of specific geometrical parameters and mechanical properties, rather than a contradiction with physical reality. The key to reconciling them is to clarify that “statistical significance” depends on the detectability of effects within the experimental design and modeling framework, while “physical existence” refers to the inherent synergistic behavior of the factors. These two concepts are complementary rather than mutually exclusive.
Through response-surface analysis of the parameters, A 3D response-surface map (as shown in [Figure 7](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f007)) and a contour map (as shown in [Figure 7](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f007)) were obtained. As shown in [Figure 7](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f007), in the three-dimensional response-surface plot, the horizontal axis represents parameters A and B respectively, and the vertical axis represents the corresponding response values. The shape and gradient changes in the response surface intuitively reflect the variation pattern of the response values under the combined effect of the two factors: the red area corresponds to higher response values, the blue area corresponds to lower response values, and the degree of color gradient clearly shows the increase or decrease range of the response values.
In the contour plot, the same pattern is also presented, that is, the color that is closer to red indicates a higher response value, and the color that is closer to blue indicates a lower response value. Overall, the response-surface plot and the contour plot can intuitively, comprehensively and quantitatively reveal the spatial distribution characteristics, interaction influence rules and extreme distribution areas of the response values with respect to parameters A and B, and can provide reliable visual evidence and theoretical support for the subsequent identification of key parameters, analysis of the influencing mechanism and selection of the optimal parameter range.
In the three-dimensional response-surface plot, the horizontal axis is A and C, and the vertical axis R1 represents the obtained response values (as shown in [Figure 8](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f008)). This surface is used to reflect the continuous pattern of response values changing with the two independent variables. The steeper the surface, the stronger the influence of the corresponding independent variable on the response value in that area. If the surface is flat, the influence of the variable on the response value is relatively weak. For the contour plot, the color represents the gradual change (green-yellow-red), indicating high or low response values. The red area has the highest response value, while the green area has the lowest response value. Points on the same curve have the same response value. The denser the curve, the more intense the change in the response value in that area. Within the 1.5 range, the contours are elliptical, indicating a significant interaction between the two independent variables and the change in response values.
[Figure 9](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f009) shows the three-dimensional response-surface plot and contour plot corresponding to the joint variation in independent variables B and C. From the figure, it can be seen that the response value shows a continuous and smooth trend as variables B and C change. The color in the figure gradually transitions from blue to red, intuitively reflecting the distribution pattern of the response value from low to high. The response surface as a whole presents a morphological characteristic of rapid increase first and then tending to be stable. This indicates that when B and C increase to a certain range, the growth rate of the response value gradually slows down and tends to stabilize.
Further analysis reveals that when the value of variable B is within the range of 1.5 to 2, the response value R1 can reach the optimal range of 1.8 to 2. At the same time, from the contour plot, a significant stretched elliptical distribution feature can be clearly observed. This feature directly indicates that there is a significant interaction between independent variables B and C, and their influence on the response value is not independent of each other but shows a clear coupling effect.
Depicted in [Figure 10](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f010) are the response values of A, B, and C under the maximum stress intensity factor condition. From the figure, it can be seen that the optimal values of A, B and C are 0.4045, 1.6797 and 1.9035 respectively, and the stress intensity factor reaches its maximum value of 1.3375 at these values. Additionally, it is found that the interaction between B and C is the most significant. The combination of B and C has a more significant impact on the response value than the combinations of other factors.
[Figure 10](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f010) shows the parameter combination when the stress-intensity factor reaches its global maximum: the corresponding values of variables A, B, and C are 0.4045, 1.6797, and 1.9035, respectively, at which the peak stress intensity factor is 1.3375. Meanwhile, the response surface analysis reveals that the interaction effect between factors B and C is statistically significant (p \< 0.05), and their coupled influence on the stress intensity factor is significantly stronger than that of other factor combinations (see [Table 2](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-t002) for the significance analysis of interaction terms).
It should be clearly defined that the term “optimal parameters” here refers only to the mathematical extreme solution from the response-surface model, specifically the parameter point that maximizes the stress-intensity factor, rather than the optimal design scheme based on structural safety criteria. As a key parameter characterizing the stress-field intensity at the crack tip, the magnitude of the stress-intensity factor directly determines the crack initiation threshold and propagation rate. A higher stress-intensity factor indicates a more severe stress concentration at the crack tip, a greater probability of brittle fracture or fatigue failure, and lower structural safety. Therefore, the parameter combination shown in [Figure 9](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f009) corresponds to the most unfavorable service condition of the structure, whose core value is to provide a boundary reference for dangerous working conditions in engineering design, rather than a recommended safe parameter range.
In practical structural safety design and fracture control, the high-risk parameter combination shown in [Figure 9](https://www.mdpi.com/1996-1944/19/5/964#materials-19-00964-f009) should be avoided. We have supplemented the distinction between mathematically optimal extreme values and engineering safety-optimized solutions in the revised manuscript to prevent confusion of academic concepts, ensure the consistency between research conclusions and engineering practice, and thus provide a more valuable reference for structural safety design.
## 7\. Conclusions
This paper conducts a systematic study on the dynamic propagation of orthogonal anisotropic interface cracks in piezoelectric bimaterial structures. By introducing Laplace transform and Fourier transform to convert the governing equations and combining the Chebyshev collocation method to discretize the original problem into algebraic equations, the dynamic variation laws of the stress intensity factor at the crack tip and the electric displacement intensity factor are obtained after numerical inverse transformation. Based on this, the influence mechanism of key parameters on the fracture behavior is analyzed by using the response surface method, and the following main conclusions are obtained:
(1)
The elastic parameters have a significant regulatory effect on the dynamic propagation of interface cracks. Within a certain range, the larger the elastic parameters, the higher the dimensionless dynamic stress-intensity factor at the crack tip, indicating that the elastic properties of the material directly affect the degree of stress concentration and the driving ability of the crack.
(2)
The impact-load application time significantly changes the dynamic response characteristics of the crack. With the passage of time, the dimensionless stress-intensity factor shows a pattern of first increasing, reaching a peak, and then oscillating and converging, eventually approaching the corresponding value under static load, reflecting the combined effect of stress wave propagation, reflection, and energy dissipation.
(3)
The crack length is an important geometric factor affecting the fracture risk. The larger the crack length, the more significant the stress concentration at the tip, and the corresponding dimensionless dynamic stress-intensity factor increases. The structural fracture failure risk is higher.
(4)
For non-permeable cracks, the dual effect of forward shear stress promoting crack propagation and negative shear stress first inhibiting and then driving the crack propagation has verified the rationality of the shear-fracture calculation model established in this paper. For the dimensionless dynamic strength factors corresponding to different force–electric coupling coefficients λ, they all show a typical characteristic of rapid increase—reaching a peak—slow decrease—approaching stability as c p t / ( 2 a ) changes. This further proves that the model in this paper can accurately describe the transient fracture-mechanics behavior of piezoelectric materials under coupled loads.
(5)
The constructed prediction model has high accuracy and strong reliability, with a determination coefficient R2 of 0.9886. The real values and predicted values are highly consistent, providing an efficient and accurate data-driven method for the rapid prediction of the stress-intensity factor at the crack tip.
(6)
The key parameter combinations that maximize the stress-intensity factor were obtained through response surface analysis: A = 0.4045, B = 1.6797, C = 1.9035, corresponding to a peak stress intensity factor of 1.3375. Significance analysis indicates that there is a significant interaction between factors B and C (p \< 0.05), and their coupling effect is stronger than other parameter combinations.
(7)
The “optimal parameter combination” obtained in this paper is only a mathematical extremum point. Even the parameter combination that maximizes the stress-intensity factor corresponds to the most dangerous service condition of the structure, rather than an engineering safety optimization scheme. This result can be used to identify high-risk parameter intervals, providing important theoretical references for the fracture prevention, safety design, and boundary determination of dangerous conditions of engineering structures.
Overall, the analytical–numerical combined analysis framework established in this paper not only enriches the dynamic fracture theory of piezoelectric bimaterial interface cracks but also provides feasible methods for the safety assessment, parameter optimization, and failure prevention of intelligent materials and structures. It has certain theoretical value and practical significance for promoting the reliable application of piezoelectric structures in engineering fields. In future research, we will incorporate the thickness ratio (d1/d2) as a key parameter into our analysis and conduct systematic quantitative studies using appropriate numerical simulation software. This extension aims to more comprehensively reveal the influence law of the thickness ratio on the dynamic fracture behavior of piezoelectric bimaterial interface cracks, further improving the completeness and engineering applicability of the proposed theoretical framework.
## Author Contributions
Conceptualization, Y.Z. and J.L.; methodology, J.L.; software, J.L.; validation, Y.Z., J.L. and X.L.; formal analysis, X.L.; investigation, J.M.; resources, Y.Z. and J.M.; data curation, Y.Z.; writing—original draft preparation, Y.Z.; writing—review and editing, X.L.; visualization, J.M.; supervision, J.L.; project administration, Y.Z.; funding acquisition, J.M. and J.L. All authors have read and agreed to the published version of the manuscript.
## Funding
This research was funded by National Natural Science Foundation of China grant number 12301591, and also by the Shanxi Provincial Key Research and Development Project grant number 202102090301027. The APC was funded by the Taiyuan University of Science and Technology doctoral research start-up fund grant number 20252008.
## Data Availability Statement
The original contributions presented in this study are included in the article, further inquiries can be directed to the corresponding author.
## Conflicts of Interest
The authors declare no conflicts of interest.
## Appendix A
ϖ 1 ( k ) \= ( − ε 22 ( k ) c 44 ( k ) ( ( c 55 ( k ) c 44 ( k ) − v 2 c k ¯ 2 ) s 2 \+ 2 v p s i c k ¯ 2 \+ p 2 c k ¯ 2 ) − 2 e 24 ( k ) e 15 ( k ) s 2 c 44 ( k ) c 44 ( k ) − ε 11 ( k ) s 2 c 44 ( k ) ) / ( ( ε 22 ( k ) c 44 ( k ) \+ ( e 24 ( k ) c 44 ( k ) ) 2 ) ,
ϖ 2 ( k ) \= ( ε 11 ( k ) c 44 ( k ) ( ( c 55 ( k ) c 44 ( k ) − v 2 c k ¯ 2 ) s 4 \+ 2 v p s 3 i c k ¯ 2 \+ p 2 s 2 c k ¯ 2 ) \+ e 15 ( k ) e 15 ( k ) s 4 c 44 ( k ) c 44 ( k ) ) / ( ( ε 22 ( k ) c 44 ( k ) \+ ( e 24 ( k ) c 44 ( k ) ) 2 ) ,
ϑ 1 ( k ) \= − 72 ϖ 1 ( k ) ϖ 2 ( k ) − 2 ϖ 1 ( k ) 3 54 \+ ( 72 ϖ 1 ( k ) ϖ 2 ( k ) − 2 ϖ 1 ( k ) 3 ) 2 2916 − ( 12 ϖ 2 ( k ) \+ ϖ 1 ( k ) 2 ) 3 729 3 ,
ϑ 1 ( k ) \= − 72 ϖ 1 ( k ) ϖ 2 ( k ) − 2 ϖ 1 ( k ) 3 54 − ( 72 ϖ 1 ( k ) ϖ 2 ( k ) − 2 ϖ 1 ( k ) 3 54 ) 2 − ( 12 ϖ 2 ( k ) \+ ϖ 1 ( k ) 2 9 ) 3 3 ,
γ 1 ( k ) \= − − 2 ϖ 1 ( k ) / 3 \+ ϑ 1 ( k ) \+ ϑ 2 ( k ) \+ − 4 ϖ 1 ( k ) / 3 − ϑ 1 ( k ) − ϑ 2 ( k ) 2 ,
γ 2 ( k ) \= − − 2 ϖ 1 ( k ) / 3 \+ ϑ 1 ( k ) \+ ϑ 2 ( k ) − − 4 ϖ 1 ( k ) / 3 − ϑ 1 ( k ) − ϑ 2 ( k ) 2 ,
γ 3 ( k ) \= − 2 ϖ 1 ( k ) / 3 \+ ϑ 1 ( k ) \+ ϑ 2 ( k ) \+ − 4 ϖ 1 ( k ) / 3 − ϑ 1 ( k ) − ϑ 2 ( k ) 2 ,
γ 4 ( k ) \= − 2 ϖ 1 ( k ) / 3 \+ ϑ 1 ( k ) \+ ϑ 2 ( k ) − − 4 ϖ 1 ( k ) / 3 − ϑ 1 ( k ) − ϑ 2 ( k ) 2 ,
A \* \= ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) e s γ 1 ( 2 ) d 2 ( − χ 11 ( 1 ) χ 22 ( 2 ) \+ χ 12 ( 2 ) χ 21 ( 1 ) \+ χ 12 ( 1 ) χ 22 ( 2 ) − χ 12 ( 2 ) χ 22 ( 1 ) ) \+ ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) e s γ 2 ( 2 ) d 2 ( − χ 11 ( 2 ) χ 21 ( 1 ) \+ χ 11 ( 1 ) χ 21 ( 2 ) − χ 12 ( 1 ) χ 21 ( 2 ) \+ χ 11 ( 2 ) χ 22 ( 1 ) ) ,
A 1 \= − ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) e s γ 1 ( 2 ) d 2 χ 22 ( 2 ) \+ ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) e s γ 2 ( 2 ) d 2 χ 21 ( 2 ) ,
B 1 \= − ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) e s γ 2 ( 2 ) d 2 χ 11 ( 2 ) \+ ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) e s γ 1 ( 2 ) d 2 χ 12 ( 2 ) ,
A 2 \= − ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) e s γ 2 ( 2 ) d 2 χ 21 ( 2 ) \+ ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) e s γ 1 ( 2 ) d 2 χ 22 ( 2 ) ,
B 2 \= − ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) e s γ 1 ( 2 ) d 2 χ 12 ( 2 ) \+ ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) e s γ 2 ( 2 ) d 2 χ 11 ( 2 ) ,
A 1 \* \= ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) e s γ 2 ( 2 ) d 2 ( χ 21 ( 1 ) − χ 22 ( 1 ) ) ,
B 1 \* \= ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) e s γ 2 ( 2 ) d 2 ( χ 12 ( 1 ) − χ 11 ( 1 ) ) ,
A 2 \* \= ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) e s γ 1 ( 2 ) d 2 ( χ 22 ( 1 ) − χ 21 ( 1 ) ) ,
B 2 \* \= ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) e s γ 1 ( 2 ) d 2 ( χ 11 ( 1 ) − χ 12 ( 1 ) ) ,
C 1 \* \= ( ( c 44 ( 1 ) χ 11 ( 1 ) \+ e 24 ( 1 ) χ 21 ( 1 ) ) γ 1 ( 1 ) − ( c 44 ( 1 ) χ 12 ( 1 ) \+ e 24 ( 1 ) χ 22 ( 1 ) ) γ 2 ( 1 ) ) ( ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) χ 21 ( 2 ) − ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) χ 22 ( 2 ) ) ,
C 2 \* \= ( ( c 44 ( 1 ) χ 11 ( 1 ) \+ e 24 ( 1 ) χ 21 ( 1 ) ) γ 1 ( 1 ) − ( c 44 ( 1 ) χ 12 ( 1 ) \+ e 24 ( 1 ) χ 22 ( 1 ) ) γ 2 ( 1 ) ) ( ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) χ 12 ( 2 ) − ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) χ 11 ( 2 ) ) ,
C 3 \* \= ( ( e 15 ( 1 ) χ 11 ( 1 ) − ε 22 ( 1 ) χ 21 ( 1 ) ) ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) − ( e 15 ( 1 ) χ 12 ( 1 ) − ε 22 ( 1 ) χ 22 ( 1 ) ) ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) ) γ 1 ( 1 ) γ 2 ( 2 ) ( χ 21 ( 1 ) − χ 22 ( 1 ) ) ,
C 4 \* \= ( ( e 15 ( 1 ) χ 11 ( 1 ) − ε 22 ( 1 ) χ 21 ( 1 ) ) ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) − ( e 15 ( 1 ) χ 12 ( 1 ) − ε 22 ( 1 ) χ 22 ( 1 ) ) ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) ) γ 1 ( 1 ) γ 2 ( 2 ) ( χ 12 ( 1 ) − χ 11 ( 1 ) ) ,
B \* \= ( c 44 ( 2 ) χ 11 ( 2 ) \+ e 24 ( 2 ) χ 21 ( 2 ) ) γ 1 ( 2 ) ( − χ 11 ( 1 ) χ 22 ( 2 ) \+ χ 12 ( 2 ) χ 21 ( 1 ) \+ χ 12 ( 1 ) χ 22 ( 2 ) − χ 12 ( 2 ) χ 22 ( 1 ) ) \+ ( c 44 ( 2 ) χ 12 ( 2 ) \+ e 24 ( 2 ) χ 22 ( 2 ) ) γ 2 ( 2 ) ( − χ 11 ( 2 ) χ 21 ( 1 ) \+ χ 11 ( 1 ) χ 21 ( 2 ) − χ 12 ( 1 ) χ 21 ( 2 ) \+ χ 11 ( 2 ) χ 22 ( 1 ) ) .
## References
1. Singh, V.; Jangid, K.; Bui, T.Q. A study of strain and electric field gradient effects on two collinear cracks in an arbitrary polarized piezoelectric layer. Eur. J. Mech.-A/Solids **2025**, 114, 105723. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=A+study+of+strain+and+electric+field+gradient+effects+on+two+collinear+cracks+in+an+arbitrary+polarized+piezoelectric+layer&author=Singh,+V.&author=Jangid,+K.&author=Bui,+T.Q.&publication_year=2025&journal=Eur.+J.+Mech.-A/Solids&volume=114&pages=105723&doi=10.1016/j.euromechsol.2025.105723)\] \[[CrossRef](https://doi.org/10.1016/j.euromechsol.2025.105723)\]
2. Mu, X.; Zhu, Z.; Zhang, L.; Gao, Y. The Green’s functions of two-dimensional piezoelectric quasicrystal semi-infinite crack and finite crack. Appl. Math. Model. **2025**, 137, 115732. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=The+Green%E2%80%99s+functions+of+two-dimensional+piezoelectric+quasicrystal+semi-infinite+crack+and+finite+crack&author=Mu,+X.&author=Zhu,+Z.&author=Zhang,+L.&author=Gao,+Y.&publication_year=2025&journal=Appl.+Math.+Model.&volume=137&pages=115732&doi=10.1016/j.apm.2024.115732)\] \[[CrossRef](https://doi.org/10.1016/j.apm.2024.115732)\]
3. Li, Y.; Yan, S.; Zhao, M.; Ren, J. Fracture Analysis of Planar Cracks in 3D Thermal Piezoelectric Semiconductors. Int. J. Mech. Sci. **2024**, 273, 109212. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Fracture+Analysis+of+Planar+Cracks+in+3D+Thermal+Piezoelectric+Semiconductors&author=Li,+Y.&author=Yan,+S.&author=Zhao,+M.&author=Ren,+J.&publication_year=2024&journal=Int.+J.+Mech.+Sci.&volume=273&pages=109212&doi=10.1016/j.ijmecsci.2024.109212)\] \[[CrossRef](https://doi.org/10.1016/j.ijmecsci.2024.109212)\]
4. Tan, Y.; Rao, W.; Wan, K.; Peng, K.; Zhao, J.; Li, X. Phase-field model for fatigue crack growth in piezoelectrics: Energetically consistent boundary condition. Int. J. Solids Struct. **2025**, 316, 113378. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Phase-field+model+for+fatigue+crack+growth+in+piezoelectrics:+Energetically+consistent+boundary+condition&author=Tan,+Y.&author=Rao,+W.&author=Wan,+K.&author=Peng,+K.&author=Zhao,+J.&author=Li,+X.&publication_year=2025&journal=Int.+J.+Solids+Struct.&volume=316&pages=113378&doi=10.1016/j.ijsolstr.2025.113378)\] \[[CrossRef](https://doi.org/10.1016/j.ijsolstr.2025.113378)\]
5. Nourazar, M.; Yang, W.; Chen, Z. Fracture analysis of a curved crack in a piezoelectric plane under general thermal loading. Eng. Fract. Mech. **2023**, 284, 109208. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Fracture+analysis+of+a+curved+crack+in+a+piezoelectric+plane+under+general+thermal+loading&author=Nourazar,+M.&author=Yang,+W.&author=Chen,+Z.&publication_year=2023&journal=Eng.+Fract.+Mech.&volume=284&pages=109208&doi=10.1016/j.engfracmech.2023.109208)\] \[[CrossRef](https://doi.org/10.1016/j.engfracmech.2023.109208)\]
6. Cao, T.; Zhang, Y.; Qin, T. Research on interaction of two straight cracks embedded in different spaces of piezoelectric biomaterial. Eng. Anal. Bound. Elem. **2025**, 179, 106401. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Research+on+interaction+of+two+straight+cracks+embedded+in+different+spaces+of+piezoelectric+biomaterial&author=Cao,+T.&author=Zhang,+Y.&author=Qin,+T.&publication_year=2025&journal=Eng.+Anal.+Bound.+Elem.&volume=179&pages=106401&doi=10.1016/j.enganabound.2025.106401)\] \[[CrossRef](https://doi.org/10.1016/j.enganabound.2025.106401)\]
7. Yu, H.; Zhu, S.; Ma, H.; Wang, J. Interface crack analysis of piezoelectric laminates considering initial strain. Int. J. Mech. Sci. **2024**, 271, 109104. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Interface+crack+analysis+of+piezoelectric+laminates+considering+initial+strain&author=Yu,+H.&author=Zhu,+S.&author=Ma,+H.&author=Wang,+J.&publication_year=2024&journal=Int.+J.+Mech.+Sci.&volume=271&pages=109104&doi=10.1016/j.ijmecsci.2024.109104)\] \[[CrossRef](https://doi.org/10.1016/j.ijmecsci.2024.109104)\]
8. Zhu, S.; Yu, H.; Guo, L. Analysis of an interfacial crack between two nonhomogeneous piezoelectric materials using a new domain-independent interaction integral. Compos. Struct. **2024**, 331, 117873. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Analysis+of+an+interfacial+crack+between+two+nonhomogeneous+piezoelectric+materials+using+a+new+domain-independent+interaction+integral&author=Zhu,+S.&author=Yu,+H.&author=Guo,+L.&publication_year=2024&journal=Compos.+Struct.&volume=331&pages=117873&doi=10.1016/j.compstruct.2024.117873)\] \[[CrossRef](https://doi.org/10.1016/j.compstruct.2024.117873)\]
9. Zhang, Y.; Li, J.; Xie, X. SH-wave scattering by the interface crack of piezoelectric ceramic polymer composites. J. Eng. Res. **2024**, 12, 259–267. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=SH-wave+scattering+by+the+interface+crack+of+piezoelectric+ceramic+polymer+composites&author=Zhang,+Y.&author=Li,+J.&author=Xie,+X.&publication_year=2024&journal=J.+Eng.+Res.&volume=12&pages=259%E2%80%93267&doi=10.1016/j.jer.2023.08.003)\] \[[CrossRef](https://doi.org/10.1016/j.jer.2023.08.003)\]
10. Zhao, S.; Govorukha, V.; Sheveleva, A.; Loboda, V. On energy release rate for an electrically permeable interface crack between two different 1D hexagonal piezoelectric QCs. Procedia Struct. Integr. **2023**, 51, 69–75. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=On+energy+release+rate+for+an+electrically+permeable+interface+crack+between+two+different+1D+hexagonal+piezoelectric+QCs&author=Zhao,+S.&author=Govorukha,+V.&author=Sheveleva,+A.&author=Loboda,+V.&publication_year=2023&journal=Procedia+Struct.+Integr.&volume=51&pages=69%E2%80%9375&doi=10.1016/j.prostr.2023.10.069)\] \[[CrossRef](https://doi.org/10.1016/j.prostr.2023.10.069)\]
11. Loboda, V.; Sheveleva, A.; Komarov, O.; Chapelle, F.; Lapusta, Y. Arbitrary number of electrically permeable cracks on the interface between two one-dimensional piezoelectric quasicrystals with piezoelectric effect. Eng. Fract. Mech. **2022**, 276, 108878. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Arbitrary+number+of+electrically+permeable+cracks+on+the+interface+between+two+one-dimensional+piezoelectric+quasicrystals+with+piezoelectric+effect&author=Loboda,+V.&author=Sheveleva,+A.&author=Komarov,+O.&author=Chapelle,+F.&author=Lapusta,+Y.&publication_year=2022&journal=Eng.+Fract.+Mech.&volume=276&pages=108878&doi=10.1016/j.engfracmech.2022.108878)\] \[[CrossRef](https://doi.org/10.1016/j.engfracmech.2022.108878)\]
12. Sheveleva, A.; Loboda, V.; Lapusta, Y. A conductive crack and a remote electrode at the interface between two piezoelectric materials. Appl. Math. Model. **2020**, 87, 287–299. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=A+conductive+crack+and+a+remote+electrode+at+the+interface+between+two+piezoelectric+materials&author=Sheveleva,+A.&author=Loboda,+V.&author=Lapusta,+Y.&publication_year=2020&journal=Appl.+Math.+Model.&volume=87&pages=287%E2%80%93299&doi=10.1016/j.apm.2020.06.003)\] \[[CrossRef](https://doi.org/10.1016/j.apm.2020.06.003)\]
13. Sun, W.; Qu, W.; Gu, Y. Arbitrary-order GFDM for crack analysis of homogeneous and bi-material piezoelectric problems. Theor. Appl. Fract. Mech. **2025**, 139, 105107. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Arbitrary-order+GFDM+for+crack+analysis+of+homogeneous+and+bi-material+piezoelectric+problems&author=Sun,+W.&author=Qu,+W.&author=Gu,+Y.&publication_year=2025&journal=Theor.+Appl.+Fract.+Mech.&volume=139&pages=105107&doi=10.1016/j.tafmec.2025.105107)\] \[[CrossRef](https://doi.org/10.1016/j.tafmec.2025.105107)\]
14. Singh, R. Inspection of dynamic fracture behavior of multiple interfacial cracks emanating from circular holes in functionally graded piezoelectric bi-materials. Appl. Math. Model. **2025**, 143, 116041. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Inspection+of+dynamic+fracture+behavior+of+multiple+interfacial+cracks+emanating+from+circular+holes+in+functionally+graded+piezoelectric+bi-materials&author=Singh,+R.&publication_year=2025&journal=Appl.+Math.+Model.&volume=143&pages=116041&doi=10.1016/j.apm.2025.116041)\] \[[CrossRef](https://doi.org/10.1016/j.apm.2025.116041)\]
15. Yan, Z.; Wen, C.; Feng, W.; Zhang, C. Dynamic growth properties of an internal crack in piezoelectric plates based on the improved mechanical energy release rate criterion. Theor. Appl. Fract. Mech. **2025**, 140, 105131. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+growth+properties+of+an+internal+crack+in+piezoelectric+plates+based+on+the+improved+mechanical+energy+release+rate+criterion&author=Yan,+Z.&author=Wen,+C.&author=Feng,+W.&author=Zhang,+C.&publication_year=2025&journal=Theor.+Appl.+Fract.+Mech.&volume=140&pages=105131&doi=10.1016/j.tafmec.2025.105131)\] \[[CrossRef](https://doi.org/10.1016/j.tafmec.2025.105131)\]
16. Zhu, S.; Yu, H.; Wang, Z. Interfacial dynamic impermeable crack analysis in dissimilar piezoelectric materials by a new interaction integral. Compos. Struct. **2025**, 352, 118668. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Interfacial+dynamic+impermeable+crack+analysis+in+dissimilar+piezoelectric+materials+by+a+new+interaction+integral&author=Zhu,+S.&author=Yu,+H.&author=Wang,+Z.&publication_year=2025&journal=Compos.+Struct.&volume=352&pages=118668&doi=10.1016/j.compstruct.2024.118668)\] \[[CrossRef](https://doi.org/10.1016/j.compstruct.2024.118668)\]
17. Liu, H.; Zhu, S. Dynamic analysis of two collinear permeable Mode-I cracks in piezoelectric materials based on non-local piezoelectricity theory. Multidiscip. Model. Mater. Struct. **2019**, 15, 1274–1293. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+analysis+of+two+collinear+permeable+Mode-I+cracks+in+piezoelectric+materials+based+on+non-local+piezoelectricity+theory&author=Liu,+H.&author=Zhu,+S.&publication_year=2019&journal=Multidiscip.+Model.+Mater.+Struct.&volume=15&pages=1274%E2%80%931293&doi=10.1108/MMMS-12-2018-0215)\] \[[CrossRef](https://doi.org/10.1108/MMMS-12-2018-0215)\]
18. Liu, H.-T. Dynamic non-local theory solution to a permeable mode-I crack in a piezoelectric medium. Eng. Fract. Mech. **2017**, 179, 43–58. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+non-local+theory+solution+to+a+permeable+mode-I+crack+in+a+piezoelectric+medium&author=Liu,+H.-T.&publication_year=2017&journal=Eng.+Fract.+Mech.&volume=179&pages=43%E2%80%9358&doi=10.1016/j.engfracmech.2017.04.023)\] \[[CrossRef](https://doi.org/10.1016/j.engfracmech.2017.04.023)\]
19. Liu, H.-T.; Wu, J.-G.; Li, T.-J. Dynamic analytical solution of a limited-permeable mode-I crack in piezoelectric materials based on the non-local theory. Wave Motion **2019**, 90, 82–98. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+analytical+solution+of+a+limited-permeable+mode-I+crack+in+piezoelectric+materials+based+on+the+non-local+theory&author=Liu,+H.-T.&author=Wu,+J.-G.&author=Li,+T.-J.&publication_year=2019&journal=Wave+Motion&volume=90&pages=82%E2%80%9398&doi=10.1016/j.wavemoti.2019.05.003)\] \[[CrossRef](https://doi.org/10.1016/j.wavemoti.2019.05.003)\]
20. Afshar, H.; Bagheri, R. Several embedded cracks in a functionally graded piezoelectric strip under dynamic loading. Comput. Math. Appl. **2018**, 76, 534–550. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Several+embedded+cracks+in+a+functionally+graded+piezoelectric+strip+under+dynamic+loading&author=Afshar,+H.&author=Bagheri,+R.&publication_year=2018&journal=Comput.+Math.+Appl.&volume=76&pages=534%E2%80%93550&doi=10.1016/j.camwa.2018.04.035)\] \[[CrossRef](https://doi.org/10.1016/j.camwa.2018.04.035)\]
21. Wünsche, M.; Sladek, J.; Sladek, V.; Zhang, C.; García-Sánchez, F.; Sáez, A. Dynamic crack analysis in piezoelectric solids under time-harmonic loadings with a symmetric Galerkin boundary element method. Eng. Anal. Bound. Elem. **2017**, 84, 141–153. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+crack+analysis+in+piezoelectric+solids+under+time-harmonic+loadings+with+a+symmetric+Galerkin+boundary+element+method&author=W%C3%BCnsche,+M.&author=Sladek,+J.&author=Sladek,+V.&author=Zhang,+C.&author=Garc%C3%ADa-S%C3%A1nchez,+F.&author=S%C3%A1ez,+A.&publication_year=2017&journal=Eng.+Anal.+Bound.+Elem.&volume=84&pages=141%E2%80%93153&doi=10.1016/j.enganabound.2017.08.013)\] \[[CrossRef](https://doi.org/10.1016/j.enganabound.2017.08.013)\]
22. Yu, T.; Bui, T.Q.; Liu, P.; Zhang, C.; Hirose, S. Interfacial dynamic impermeable cracks analysis in dissimilar piezoelectric materials under coupled electromechanical loading with the extended finite element method. Int. J. Solids Struct. **2015**, 67–68, 205–218. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Interfacial+dynamic+impermeable+cracks+analysis+in+dissimilar+piezoelectric+materials+under+coupled+electromechanical+loading+with+the+extended+finite+element+method&author=Yu,+T.&author=Bui,+T.Q.&author=Liu,+P.&author=Zhang,+C.&author=Hirose,+S.&publication_year=2015&journal=Int.+J.+Solids+Struct.&volume=67%E2%80%9368&pages=205%E2%80%93218&doi=10.1016/j.ijsolstr.2015.03.037)\] \[[CrossRef](https://doi.org/10.1016/j.ijsolstr.2015.03.037)\]
23. Pak, Y.E. Crack extension force in a piezoelectric material. J. Appl. Mech. **1990**, 57, 647–653. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Crack+extension+force+in+a+piezoelectric+material&author=Pak,+Y.E.&publication_year=1990&journal=J.+Appl.+Mech.&volume=57&pages=647%E2%80%93653&doi=10.1115/1.2897071)\] \[[CrossRef](https://doi.org/10.1115/1.2897071)\]
24. Wang, X.-Y.; Yu, S.-W. Transient response of a crack in piezoelectric strip subjected to the mechanical and electrical impacts:mode-III problem. Int. J. Solids Struct. **2000**, 37, 5795–5808. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Transient+response+of+a+crack+in+piezoelectric+strip+subjected+to+the+mechanical+and+electrical+impacts:mode-III+problem&author=Wang,+X.-Y.&author=Yu,+S.-W.&publication_year=2000&journal=Int.+J.+Solids+Struct.&volume=37&pages=5795%E2%80%935808&doi=10.1016/S0020-7683\(99\)00268-1)\] \[[CrossRef](https://doi.org/10.1016/S0020-7683\(99\)00268-1)\]
25. Wang, B.-L.; Han, J.-C.; Du, S.-Y. Dynamic response for non-homogeneous piezoelectric medium with multiple cracks. Eng. Fract. Mech. **1998**, 61, 607–617. \[[Google Scholar](https://scholar.google.com/scholar_lookup?title=Dynamic+response+for+non-homogeneous+piezoelectric+medium+with+multiple+cracks&author=Wang,+B.-L.&author=Han,+J.-C.&author=Du,+S.-Y.&publication_year=1998&journal=Eng.+Fract.+Mech.&volume=61&pages=607%E2%80%93617&doi=10.1016/S0013-7944\(98\)00090-3)\] \[[CrossRef](https://doi.org/10.1016/S0013-7944\(98\)00090-3)\]
**Figure 1.** Structure of the piezoelectric bimaterials model with cracks (The green part in the picture indicates the direction of the crack movement).
**Figure 1.** Structure of the piezoelectric bimaterials model with cracks (The green part in the picture indicates the direction of the crack movement).
**Figure 2.** The variation of K with respect to d 1 / a ((**a**): ς c \= 2; (**b**): ς c \= 1\.5).
**Figure 2.** The variation of K with respect to d 1 / a ((**a**): ς c \= 2; (**b**): ς c \= 1\.5).
**Figure 3.** (**a**) The 2D plot that K varies with c p t / a; (**b**) The three-dimensional graph that K varies with c p t / a.
**Figure 3.** (**a**) The 2D plot that K varies with c p t / a; (**b**) The three-dimensional graph that K varies with c p t / a.
**Figure 4.** The variation in the normalized energy release rate under far-field stress loading.
**Figure 4.** The variation in the normalized energy release rate under far-field stress loading.
**Figure 5.** The variation in the normalized intensity factor under different λ values.
**Figure 5.** The variation in the normalized intensity factor under different λ values.
**Figure 6.** Comparison chart of real values and predicted values.
**Figure 6.** Comparison chart of real values and predicted values.
**Figure 7.** The three-dimensional graph and contour map of the interaction between A and B in terms of the influence on the stress intensity factor.
**Figure 7.** The three-dimensional graph and contour map of the interaction between A and B in terms of the influence on the stress intensity factor.
**Figure 8.** The three-dimensional graph and contour map of the interaction between A and C in terms of the influence on the stress-intensity factor.
**Figure 8.** The three-dimensional graph and contour map of the interaction between A and C in terms of the influence on the stress-intensity factor.
**Figure 9.** The three-dimensional graph and contour map of the interaction between B and C in terms of the influence on the stress intensity factor.
**Figure 9.** The three-dimensional graph and contour map of the interaction between B and C in terms of the influence on the stress intensity factor.
**Figure 10.** The response values of A, B and C when the stress-intensity factor reaches its maximum value (The red and grey dots represent the corresponding data values).
**Figure 10.** The response values of A, B and C when the stress-intensity factor reaches its maximum value (The red and grey dots represent the corresponding data values).
**Table 1.** Variance analysis of predictive models.
**Table 1.** Variance analysis of predictive models.
| Standard Error | Correlation Coefficient | Correction Coefficient | F | p | Significant |
|---|---|---|---|---|---|
| 0\.0285 | 0\.9975 | 0\.9886 | 706\.4 | \<0.0001 | significant |
**Table 2.** Variance analysis of predictive model of response value.
**Table 2.** Variance analysis of predictive model of response value.
| Source | Sum of Square | Df | Mean Square | F-Value | p\-Value | Significant |
|---|---|---|---|---|---|---|
| A-A | 0\.1013 | 1 | 0\.1013 | 125\.03 | \<0.0001 | significant |
| B-B | 3\.42 | 1 | 3\.42 | 4214\.47 | \<0.0001 | significant |
| C-C | 0\.0092 | 1 | 0\.0092 | 11\.32 | 0\.0120 | |
| AB | 0\.0399 | 1 | 0\.0399 | 49\.29 | 0\.0002 | |
| AC | 0\.0025 | 1 | 0\.0025 | 3\.09 | 0\.1224 | |
| BC | 0\.0012 | 1 | 0\.0012 | 1\.47 | 0\.2642 | |
| A2 | 0\.2894 | 1 | 0\.2894 | 357\.08 | \<0.0001 | significant |
| B2 | 0\.8693 | 1 | 0\.8693 | 1072\.74 | \<0.0001 | significant |
| C2 | 0\.2735 | 1 | 0\.2735 | 337\.47 | \<0.0001 | significant |
| Residual | 0\.0057 | 7 | 0\.0008 | | | |
| Lack of Fit | 0\.0034 | 3 | 0\.0011 | 2\.06 | 0\.2480 | not significant |
| Pure Error | 0\.0022 | 4 | 0\.0006 | | | |
| Cor Total | 5\.16 | 16 | | | | |
| | |
|---|---|
| | **Disclaimer/Publisher’s Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2026 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the [Creative Commons Attribution (CC BY) license](https://creativecommons.org/licenses/by/4.0/). |
| Shard | 31 (laksa) |
| Root Hash | 11974975279136771031 |
| Unparsed URL | com,mdpi!www,/1996-1944/19/5/964 s443 |