# Fracture Analysis of Griffith Interface Crack in Fine-Grained Piezoelectric Coating/Substrate under Thermal Loading

Authors: Shuaishuai Hu; Jiansheng Liu; Junlin Li
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

Coating often plays a role in monitoring and protecting substrates in engineering applications. Interface cracks between the coating and the substrate can lead to crack growth under the action of external loading and will cause device failure. In this paper, the behavior of a fine-grained piezoelectric coating/substrate with a Griffith interface crack under steady-state thermal loading is studied. The temperature field, displacement field, and electric field of the coupling of thermal and electromechanical problems are constructed via integral transformation and the principle of superposition. Thus, problems are transformed into a system of singular integral equations, and the expressions of thermal intensity factor, thermal stress intensity factor, and electric displacement intensity factor are obtained. We used a numerical calculation and a system of singular equations to obtain the relationship of strength factor with material parameters, coating thickness, and crack size.


---

## Body

## 1. Introduction

Piezoelectric materials have played an important role in the production of intelligent structures. Due to strength demand, we often bond another piezoelectric material on the surface of the device to protect it or monitor the device. However, piezoelectric composites are prone to cracks at the interface during the process of fabricating defects and load conditions. This can even lead to structural failure. Therefore, problems of interface fracture of the piezoelectric composite under all kinds of loads have attracted wide concern [1–9]. Qin and Mai [10] analyzed the interface crack problem of a piezoelectric bimaterial subjected to combined thermal, mechanical, and electrical loads via Stroh’s formalism and the singular integral equation method. However, the influence of material size on the thermal intensity factor has not been considered in the numerical analysis. Ueda [11, 12] investigated the problem of a parallel crack in a piezoelectric strip under thermoelectric loading and the FGM material strip containing an embedded crack or an edge crack perpendicular to its boundaries. Wang and Noda [13] studied the problems of the piezoelectric material strip with a Griffith interface crack under thermal loading and discussed the effects of polarization direction, crack size, and location on the thermal strength factor. Kang et al. [14] considered the influences of piezoelectric laminated composite shells with interlaminar stresses under electrical, thermal, and mechanical loads. They concluded that the interlaminar shear stress of the laminated shell can be reduced by selecting the appropriate electric field value. Cook and Vel [15] examined a comprehensive multiscale analysis of laminated plates with integrated piezoelectric fiber composite actuators. They used the coupling method of microscale field variables and macroscale field variables to solve the three-dimensional macroscopic equilibrium equation of the piezoelectric laminated plate under arbitrary boundary conditions solved. Li et al. [16] presented a semianalytical technique to solve fracture behaviors of piezoelectric composites under thermal loading. This method could analytically represent the resulting stress and electric displacement distribution along the radial direction obtained. Kuuttia and Virkkunen [17] justified the surface crack behavior under periodic thermal loads. The finite element method and weight function method could simulate the crack behavior over the entire load cycle. Mojahedin et al. [18] considered the mechanical behavior of a solid circular plate which is made of saturated and unsaturated porous material with piezoelectric actuators under thermal loading. They discussed the effect of porosity on the thermal and mechanical stability of the material. Herrmann and Loboda [19] performed fracture analysis on an electrically impermeable interface crack with contact zones in thermopiezoelectrical bimaterials. Using the admissible directions of the heat and the electrical fluxes, the dependencies of the electrical intensity factors on the intensities of the thermal and electrical fluxes were discussed.Many researchers have studied the interface fracture problems of piezoelectric materials under loading [20–24]. However, all of polycrystalline materials studied above are composed of multidomain large grains, and such sizes cannot meet the requirements of the current intelligent devices, but fine-grained piezoelectric materials are gaining increasing interest. However, there are few studies on the mechanical properties of fine-grained piezoelectric materials.As an important fracture parameter that can measure structural safety, the intensity factor plays a crucial role in checking the safety of structures. In this article, a Griffith interface crack in fine-grained piezoelectric coating/substrate under steady-state thermal loading is established. An integral transformation method and superposition principle of the solution of the equation are used to transform the thermal load problem into a singular integral equation to obtain the intensity factor conveniently. Furthermore, the interaction between the intensity factor and material parameters is investigated. The results indicate that the larger elastic modulus and thinner coating thickness improve the safety of the coating/substrate structure.


## 2. Problem Formulation

Figure1 shows that a fine-grained ceramic powder coating and the substrate are bond together by plasma spraying, while the crack with a length of 2l is along the interface between the coating and the substrate. A fine-grained piezoelectric coating/substrate structure was obtained by the polarization treatment of the coating. The coating and substrate are polarized along the Y-axis, and both are transversely isotropic. The thicknesses of the coating and substrate are h1 and h2, respectively. We assumed that the crack faces remain thermally and electrically insulated, and Ta and Tb are the environment temperatures.Figure 1
Mechanical structure of fine-grained piezoelectric coating/substrate with interface crack.The constitutive equations for the elastic field are [12](1)σxx=c11m∂w∂x+c13m∂u∂y+e31m∂ϕ∂y−λ11mT,(2)σyy=c13m∂w∂x+c33m∂u∂y+e33m∂ϕ∂y−λ33mT,(3)σxy=c44m∂w∂y+∂u∂x+e15m∂ϕ∂x,(4)Dx=e15m∂w∂y+∂u∂x−ε11m∂ϕ∂x,(5)Dy=e31m∂w∂x+e33m∂u∂y−ε33m∂ϕ∂y+pymT.In Equations (1), (2), (3), (4) and (5), the w and u are the displacement components; ϕ is the electric potential; T is the temperature change; σxx, σyy, and σxy are stress components; Dx and Dy are electric displacement components;c11m, c13m, and c44m are elastic modulus; e31m, e33m, e15m and ε11m, ε33m are piezoelectric and dielectric constants, respectively. λ11m and λ33m are the stress-temperature coefficients; pym is the pyroelectric constant. The superscript m(m=1,2) stands for the fine-grained piezoelectric coating and piezoelectric substrate, respectively.Assume that the temperature satisfies the Fourier heat conduction equation, as follows:(6)km2∂2T∂x2+∂2T∂y2=0,where km=kxm/kym, kxm, and kym are coefficients of thermal conductivity. m=1,2.The equations of equilibrium are [12](7)c11m∂2w∂x2+c44m∂2w∂y2+c13m+c44m∂2u∂x∂y+e31m+e15m∂2ϕ∂x∂y=λ11m∂T∂x,(8)c44m∂2u∂x2+c33m∂2w∂y2+c13m+c44m∂2w∂x∂y+e15m∂2ϕ∂x2+e33m∂2ϕ∂y2=λ33m∂T∂y,(9)e15m∂2u∂x2+e33m∂2u∂y2+e15m+e31m∂2w∂x∂y−ε11m∂2ϕ∂x2−ε33m∂2ϕ∂y2=−pym∂T∂y.Without loss of generality, the boundary conditions can be written as(10a)σyyx,0+=σyyx,0−=0,x<l,(10b)ux,0+=ux,0−,x≥l,(10c)σxyx,0+=σxyx,0−=0,x<l,(10d)wx,0+=wx,0−,x≥l,(10e)Dyx,0+=Dyx,0−=0,x<l,(10f)ϕx,0+=ϕx,0−,x≥l,(10g)σyyx,0+=σyyx,0−,x≥l,(10h)σxyx,0+=σxyx,0−,x≥l,(10i)Dyx,0+=Dyx,0−,x≥l,(10j)σyyx,h1=σyyx,−h2=0,x≥l,(10k)σxyx,h1=σxyx,−h2=0,x≥l,(10l)Dyx,h1=Dyx,−h2=0,x≥l.The thermal boundary conditions are specified as(11a)Tx,0+=Tx,0−,x≥l,(11b)Tx,h1=Tax,Tx,h2=Tbx,x≥0,(11c)ky1∂Tx,0+∂y=ky2∂Tx,0−∂y=−Q0x,x<l,(11d)ky1∂Tx,0+∂y=ky2∂Tx,0−∂y,x≥l,where Tax and Tbx are environment temperatures, Q0x is the thermal flow, and ky1 and ky2 are coefficients of thermal conductivity.Since the temperature field does not depend on the stress and displacement field, it can therefore be calculated separately, and then, the thermal stress and potential displacement field can be solved.


## 3. Solution to the Problem

### 3.1. Temperature Field

By using the Fourier integral transform, the solution of temperature field can be expressed as(12)Tmx,y=12π∫−∞+∞A1mse−skmy+A2mseskmye−isxds,where km=kxm/kym, A1ms, and A2ms are the unknown functions m=1,2.To solve the problem easily, the density function is defined as follows:(13)Gx=∂∂xTmx,0+−Tmx,0−.According to Equation (11a), the dislocation density function should make the following single-valued conditions valid:(14)∫−llGtdt=0.Substituting formula (12) into Equations (11a), (11b), (11c) and (11d), we obtain(15)A11s+A21s−A12s−A22s=is∫−llGteistdt,(16)−sk1ky1A11s+sk1ky1A21s=−sk1ky2A12s+sk1ky2A21s,(17)A11se−sk1h1+A21sesk1h1=∫−∞+∞Taxeisxdx,(18)A12se−sk2h2+A22sesk2h2=∫−∞+∞Tbxeisxdx.Equations (16), (17) and (18) can be written in the form of a matrix as follows:(19)A⋅B=C,where(20)A=−sk1ky1sk1ky1sk2ky2−sk2ky211−1−1e−sk1h1esk1h10000e−sk2h2esk2h2,B=A11S,A21S,A12S,A22ST,C=0,I0,I1,I2T,I0=is∫−llGteistdt,I1=∫−llTaxeisxdx,I2=∫−llTbxeisxdx.Solving the linear Equation (19), we get(21)A11s=F1I0−F2I1+F3I2,A21=F4I1−F5I0−F6I2,A12s=F7I2+F8I0−F9I1,A22s=F10I1−F11I2−F12I0.The expressions ofFq and Δ, q=1,2,...12, are given in Appendix A. Obviously, the unknown functions A11s, A21s, A12s, and A22s depend on I0. Once I0 is determined, the temperature field could be obtained.By substituting Equation (12) into Equation (11c), we obtain(22)i2π∫−∞+∞sk1ky1F1−F5seit−xds∫−llGtdt=Q0x−d0x,where(23)d0x=12π∫−∞+∞sk1ky1F3I2−F2I1−F4I1−F6I2e−isxds.Furthermore, Equation (23) can be written as(24)1π∫0+∞Mssint−xsds∫−llGtdt=Q0x−d0x,where(25)Ms=sk1ky1F1−F5s.Through changing the order of the integration on the cracks’ surface, we obtain(26)−1π∫−llM0t−x+κ1x,tGtdt=Q0x−d0x,where(27)lims→∞Ms=M0,M0=−k1k2ky1ky2k1ky1+k2ky2,κ1x,t=∫0+∞Ms−M0sint−xsds.Equation (26) is a singular integral equation with Cauchy kernel. The density function Gt in Equation (26) has the square root-type singularity. Thus, we can express Gt as follows:(28)Gt=1l2−t2g0t,where g0t is a continuous function defined in the interval −l,l.The thermal flow intensity factors are defined by [13].(29)k0L=limx→−l−2π−x−lky∂Tmx,0∂y,(30)k0R=limx→l+2πx−lky∂Tmx,0∂y.According to Equation (26), the singular part of ky∂Tmx,0/∂y is(31)limx→−l−ky∂Tmx,0∂y=M0limx→−l−∫−llG−l−x−l,(32)limx→l+ky∂Tmx,0∂y=M0limx→l+∫−llGlx−l.Substituting Equation (28) into Equations (31) and (32), we obtain(33)k0L=M0πlg0−l,k0R=−M0πlg0l.


### 3.2. Thermal Stress and Electric Displacement Fields

Similar to the solution of linear equation, the solutions to Equations (7), (8) and (9) of equilibrium for displacements and electric potential consist of a homogeneous solution and a particular solution.Equations (7), (8) and (9) can be written into homogeneous forms as follows:(34)c11m∂2w∂x2+c44m∂2w∂y2+c13m+c44m∂2u∂x∂y+e31m+e15m∂2ϕ∂x∂y=0,(35)c44m∂2u∂x2+c33m∂2w∂y2+c13m+c44m∂2w∂x∂y+e15m∂2ϕ∂x2+e33m∂2ϕ∂y2=0,(36)e15m∂2u∂x2+e33m∂2u∂y2+e15m+e31m∂2w∂x∂y−ε11m∂2ϕ∂x2−ε33m∂2ϕ∂y2=0.By using the Fourier integral transform, the solutions of Equations (34), (35) and (36) can be expressed as [13](37)upx,y=12π∫−∞+∞∑k=16B1kmsesλkyfkse−isxds,(38)wpx,y=12π∫−∞+∞∑k=16B2kmsesλkyfkse−isxds,(39)ϕpx,y=12π∫−∞+∞∑k=16B3kmsesλkyfkse−isxds,where B1km, B2km, B3km, and fks, m=1,2, are unknown functions to be determined.The special solutions of Equations (7), (8) and (9) can be expressed as(40)uqx,y=12π∫−∞+∞1sB1amA1me−skmy+B1bmA2meskmye−isxds,(41)wqx,y=12π∫−∞+∞1sB2amA1me−skmy+B2bmA2meskmye−isxds,(42)ϕqx,y=12π∫−∞+∞1sB3amA1me−skmy+B3bmA2meskmye−isxds,where B1am, B2am, B3am, B1bm, B2bm, and B3bm satisfy the expression as follows:(43)Λ2km2−Λ1km+Λ0B1amB2amB3am=−iλ11m−is2sλ33mkms2spymkm,Λ2km2+Λ1km+Λ0B1bmB2bmB3bm=−iλ11m−is2sλ33mkms2spymkm,where(44)Λ0=−c11m000−c44m−e15m0−e15mε11m,Λ1=−is2s0c44me31m+e15mc44m00e31m+e15m00,Λ2=c44m000c33me33m0e33m−ε33m.Therefore, the solution of Equations (7), (8) and (9) can be written as(45)ux,y=12π∫−∞+∞∑k=16B1kmsesλkyfkse−isxds+12π∫−∞+∞1sB1amA1me−skmy+B1bmA2meskmye−isxds,(46)wx,y=12π∫−∞+∞∑k=16B2kmsesλkyfkse−isxds+12π∫−∞+∞1sB2amA1me−skmy+B2bmA2meskmye−isxds,(47)ϕx,y=12π∫−∞+∞∑k=16B3kmsesλkyfkse−isxds+12π∫−∞+∞1sB3amA1me−skmy+B3bmA2meskmye−isxds.Substituting Equations (45), (46) and (47) into Equations (2), (3), and (5), the thermal stress and electric displacement fields are obtained:(48)σyy=12π∫−∞+∞B1∗mf0se−isxds+12π∫−∞+∞B10∗mA0e−isxds,(49)σxy=12π∫−∞+∞B2∗mf0se−isxds+12π∫−∞+∞B20∗mA0e−isxds,(50)Dy=12π∫−∞+∞B3∗mf0se−isxds+12π∫−∞+∞B30∗mA0e−isxds,where(51)B1∗m=λ1c33mesλ1yB11m−ic31mesλ1yB21m+λ1e33mesλ1yB31mλ1c33mesλ1yB12m−ic32mesλ1yB22m+λ1e33mesλ1yB32mλ1c33mesλ1yB13m−ic33mesλ1yB23m+λ1e33mesλ1yB33mλ1c33mesλ1yB14m−ic34mesλ1yB24m+λ1e33mesλ1yB34mλ1c33mesλ1yB15m−ic35mesλ1yB25m+λ1e33mesλ1yB35mλ1c33mesλ1yB16m−ic36mesλ1yB26m+λ1e33mesλ1yB36mT,B2∗m=‐ic44mesλ1yB11m+λ1c44mesλ1yB21m−ie15mesλ1yB31m‐ic44mesλ1yB12m+λ1c44mesλ1yB22m−ie15mesλ1yB32m‐ic44mesλ1yB13m+λ1c44mesλ1yB23m−ie15mesλ1yB33m‐ic44mesλ1yB14m+λ1c44mesλ1yB24m−ie15mesλ1yB34m‐ic44mesλ1yB15m+λ1c44mesλ1yB25m−ie15mesλ1yB35m‐ic44mesλ1yB16m+λ1c44mesλ1yB26m−ie15mesλ1yB36mT,B3∗m=λ1e33mesλ1yB11m−ie31mesλ1yB21m−λ1ε33mesλ1yB31mλ1e33mesλ1yB12m−ie31mesλ1yB22m−λ1ε33mesλ1yB32mλ1e33mesλ1yB13m−ie31mesλ1yB23m−λ1ε33mesλ1yB33mλ1e33mesλ1yB14m−ie31mesλ1yB24m−λ1ε33mesλ1yB34mλ1e33mesλ1yB15m−ie31mesλ1yB25m−λ1ε33mesλ1yB35mλ1e33mesλ1yB16m−ie31mesλ1yB26m−λ1ε33mesλ1yB36mT,B10∗m=−ic13mB2am+skmc33msB1am+λ11me−skmyskmc33msB1bm−ic13mB2bm−λ11me−skmyT,B20∗m=−skmc44msB2am+ic44mB1am+ie15mB3ame−skmyskmc44msB2bm−ic44mB1bm−ie15mB3bme−skmyT,B30∗m=−skme33msB1am+skmε33msB3am+Pyme−skmyskme33msB1bm−skmε33msB3bm+Pyme−skmyT,f0=f1f2f3f4f5f6T,A0=A1mA2mT.In order to solve the thermal stress and electric displacement fields, the dislocation density functions are introduced:(52)d1x=ddxux,0+−ux,0−,(53)d2x=ddxwx,0+−wx,0−,(54)d3x=ddxϕx,0+−ϕx,0−.Substituting Equations (52), (53) and (54) into Equations (10a), (10b), (10c), (10d), (10e), (10f), (10g), (10h), (10i), (10j), (10k) and (10l), we obtain(55)∑k=16B1k1fk=H1I1∗−H2I2∗−H3I3∗+a1∗H,(56)∑k=16B1k2fk=H4I1∗−H5I2∗−H6I3∗+a2∗H,(57)∑k=16B2k1fk=H7I1∗−H8I2∗+H9I3∗+a3∗H,(58)∑k=16B2k2fk=H10I1∗−H11I2∗−H12I3∗+a4∗H,(59)∑k=16B3k1fk=H13I1∗−H14I2∗−H15I3∗+a5∗H,(60)∑k=16B3k2fk=H16I1∗−H17I2∗−H18I3∗+a6∗H,where(61)I1∗=is∫−lld1teistdt,(62)I2∗=is∫−lld2teistdt,(63)I3∗=is∫−lld3teistdt.The expressions ofH, ap∗, and Hq, p=1,2,⋯,6, q=1,2,⋯,18, are given in Appendix B. Similar to the process of solving the temperature field, once I1∗, I2∗, and I3∗ are determined, the thermal stress and electric displacement fields could be obtained.Substituting Equations (55), (56), (57), (58), (59) and (60) into Equations (48), (49) and (50), we get(64)σyyx,0−=limy→0−12π∫−∞+∞∑k=16sesλkyH−ic132H10+λkc332H4+λke332H16I1∗+−ic132H11−λkc332H5−λke332H17I2∗+ic132H12−λkc332H6−λke332H18I3∗e−isxds+12π∫−∞+∞−ic132B2aA12+B2bA22−λ112A12+A22+sk2e332s⋅−B3aA12+B3bA22+sk2c332s−B1aA12+B1bA22e−isxds=0,(65)σxyx,0−=limy→0−12π∫−∞+∞∑k=16sesλkyHλkc442H10−iλkc442H4−ie152H16I1∗+ic442H5−λkc442H11+ie152H17I2∗+ic442H6−λkc442H12+ie152H18I3∗e−isxds+12π∫−∞+∞−ic442B1aA12+B1bA22+sk2c442s−B2aA12+B2bA22−ie152−B3aA12+B3bA22e−isxds=0,(66)Dyx,0−=limy→0−12π∫−∞+∞∑k=16sesλkyHλke332H4−ie312H10−λkε332H16I1∗+ie312H11−λke332H5+λkε332H17I2∗+ie312H12−λke332H6+λkε332H18I3∗e−isxds+12π∫−∞+∞−ie312B2aA12+B2bA22−py2A12+A22−sk2ε332s⋅−B3aA12+B3bA22+sk2e332s−B1aA12+B1bA22e−isxds=0.Equations (64), (65) and (66) are divided into two parts, as follows:(67)σyyx,0−=σyy1x,0−+σyy2x,0−,σxyx,0−=σxy1x,0−+σxy2x,0−,Dyx,0−=Dy1x,0−+Dy2x,0−,where(68)σyy1x,0−=limy→0−12π∫−∞+∞∑k=16sesλkyH−ic132H10+λkc332H4+λke332H16I1∗+−ic132H11−λkc332H5−λke332H17I2∗+ic132H12−λkc332H6−λke332H18I3∗e−isxds,(69)σxy1x,0−=limy→0−12π∫−∞+∞∑k=16sesλkyHλkc442H10−iλkc442H4−ie152H16I1∗+ic442H5−λkc442H11+ie152H17I2∗+ic442H6−λkc442H12+ie152H18I3∗e−isxds,(70)Dy1x,0−=limy→0−12π∫−∞+∞∑k=16sesλkyHλke332H4−ie312H10−λkε332H16I1∗+ie312H11−λke332H5+λkε332H17I2∗+ie312H12−λke332H6+λkε332H18I3∗e−isxds,(71)σyy2x,0−=12π∫−∞+∞−ic132B2aA12+B2bA22−λ112A12+A22+sk2e332s⋅−B3aA12+B3bA22+sk2c332s−B1aA12+B1bA22e−isxds,(72)σxy2x,0−=12π∫−∞+∞−ic442B1aA12+B1bA22+sk2c442s−B2aA12+B2bA22−ie152−B3aA12+B3bA22e−isxds,(73)Dy2x,0−=12π∫−∞+∞−ie312B2aA12+B2bA22−py2A12+A22−sk2ε332s⋅−B3aA12+B3bA22+sk2e332s−B1aA12+B1bA22e−isxds.Obviously, it is easy to solve the solution of Equations (71), (72) and (73) according to Equations (39), (40), (41) and (42). Next, we need to transform Equations (68), (69) and (70) into singular integral equations of the first kind with Cauchy kernel.In another form, integral Equations (68), (69) and (70) will be(74)σyy1x,0−=limy→0−12π∫−∞+∞sesλyHQ1I1∗+Q2I2∗+Q3I3∗e−isxds,(75)σxy1x,0−=limy→0−12π∫−∞+∞sesλyHQ4I1∗+Q5I2∗+Q6I3∗e−isxds,(76)Dy1x,0−=limy→0−12π∫−∞+∞sesλyHQ7I1∗+Q8I2∗+Q9I3∗e−isxds.The expressions ofesλy and Qp, p=1,2,⋯,9, are given in Appendix C.Substituting Equations (61), (62) and (63) into Equations (74), (75) and (76), we get(77)12π∫−ll∫−∞+∞iesλyHQ1eist−xdsd1tdt+12π∫−ll∫−∞+∞iesλyHQ2eist−xdsd2tdt+12π∫−ll∫−∞+∞iesλyHQ3eist−xdsd3tdt=σ1,12π∫−ll∫−∞+∞iesλyHQ4eist−xdsd4tdt+12π∫−ll∫−∞+∞iesλyHQ5eist−xdsd5tdt+12π∫−ll∫−∞+∞iesλyHQ6eist−xdsd6tdt=σ2,12π∫−ll∫−∞+∞iesλyHQ7eist−xdsd7tdt+12π∫−ll∫−∞+∞iesλyHQ8eist−xdsd8tdt+12π∫−ll∫−∞+∞iesλyHQ9eist−xdsd9tdt=σ3,where(78)σ1=−12π∫−∞+∞−ic132B2aA12e−sk2y+B2bA22esk2y−λ112A12+A22+sk2e332s⋅−B3aA12e−sk2y+B3bA22esk2y+sk2c332s−B1aA12e−sk2y+B1bA22esk2ye−isxds,σ2=−12π∫−∞+∞−ic442B1aA12+B1bA22+sk2c442s−B2aA12+B2bA22−ie152−B3aA12+B3bA22e−isxds,σ3=−12π∫−∞+∞−ie312B2aA12e−sk2y+B2bA22esk2y−py2A12+A22−sk2ε332s⋅−B3aA12e−sk2y+B3bA22esk2y+sk2e332s−B1aA12e−sk2y+B1bA22esk2ye−isxds.On the cracks’ surface,(79)∫−llQ1∗H∗1t−x+κ2x,td1tdt+∫−llQ2∗H∗1t−x+κ3x,td2tdt+∫−llQ3∗H∗1t−x+κ4x,td3tdt=πσ1,(80)∫−llQ4∗H∗1t−x+κ5x,td1tdt+∫−llQ5∗H∗1t−x+κ6x,td2tdt+∫−llQ6∗H∗1t−x+κ7x,td3tdt=πσ2,(81)∫−llQ7∗H∗1t−x+κ8x,td1tdt+∫−llQ8∗H∗1t−x+κ9x,td2tdt+∫−llQ9∗H∗1t−x+κ10x,td3tdt=πσ3,where(82)κ2x,t=∫0∞Q1sH−Q1∗H∗sinst−xds,κ3x,t=∫0∞Q2sH−Q2∗H∗sinst−xds,κ4x,t=∫0∞Q3sH−Q3∗H∗sinst−xds,κ5x,t=∫0∞Q4sH−Q4∗H∗sinst−xds,κ6x,t=∫0∞Q5sH−Q5∗H∗sinst−xds,κ7x,t=∫0∞Q6sH−Q6∗H∗sinst−xds,κ8x,t=∫0∞Q7sH−Q7∗H∗sinst−xds,κ9x,t=∫0∞Q8sH−Q8∗H∗sinst−xds,κ10x,t=∫0∞Q9sH−Q9∗H∗sinst−xds.The expressions ofH∗ and Qp∗, p=1,2,⋯,10, are given in Appendix C. Ultimately, Equations (68), (69) and (70) are reduced to coupling the first kind Cauchy singular integral equations.In Equations (61), (62) and (63), the unknown function dkt, k=1,2,3 with the following conditions:(83)∫−lldktdt=0,k=1,2,3.Considering that the solution of the singular integral equation will use Chebychev polynomials which are defined in the interval [−1, 1], we introduce the following normalized quantities:(84)r=xl,s=tl,(85)dkx=dk~x,Lir,s=l⋅κix,t,k=1,2,3,i=1,⋯,10.Using Equations (84) and (85), integral Equations (79), (80) and (81) may be expressed as(86)∫−11Q1∗H∗1s−r+L2s,rd~1sds+∫−11Q2∗H∗1s−r+L3s,rd~2sds+∫−11Q3∗H∗1s−r+L4s,rd~3sds=πσ1,∫−11Q4∗H∗1s−r+L5s,rd~1sds+∫−11Q5∗H∗1s−r+L6s,rd~2sds+∫−11Q6∗H∗1s−r+L7s,rd~3sds=πσ2,∫−11Q7∗H∗1s−r+L8s,rd~1sds+∫−11Q8∗H∗1s−r+L9s,rd~2sds+∫−11Q9∗H∗1s−r+L10s,rd~3sds=πσ3.The crack-tip behavior can be characterized by square-root singularity weighed down by the oscillation. However, because the zone of oscillation is very small [25], the unknown function dk~s, k=1,2,3, can be expressed as follows:(87)dk~s=11−s2gks,k=1,2,3,where gks, k=1,2,3, are continuous functions defined in the interval −1,1.Using the Lobatto-Chebyshev collocation method, one can transform Equations (79), (80), (81) and (83) into a system of algebraic equations:(88)1N∑R=0NχRQ1∗H∗1sR−rq+L2sR,rqd~1sR+Q2∗H∗1sR−rq+L3sR,rqd~2sR+Q3∗H∗1sR−rq+L4sR,rqd~3sR=πσ1,1N∑R=0NχRQ4∗H∗1sR−rq+L5sR,rqd~1sR+Q5∗H∗1sR−rq+L6sR,rqd~2sR+Q6∗H∗1sR−rq+L7sR,rqd~3sR=πσ2,1N∑R=0NχRQ7∗H∗1sR−rq+L8sR,rqd~1sR+Q8∗H∗1sR−rq+L9sR,rqd~2sR+Q9∗H∗1sR−rq+L10sR,rqd~3sR=πσ3,∑R=0NχRgksR=0,k=1,2,3,q=1,2,⋯,N,where(89)χ0=χN=12,χ1=⋯=χN−1=1,rq=cos2q−1π2N,sR=cosRπN.N is the node number of quadrature formula; sR and rq are the zero points of the first and second kinds of Chebyshev polynomials. According to Theocaris and Ioakimidis [26] as well as Lu [27], the convergence rate of the Lobatto-Chebyshev collocation method is first order (i.e., oh, where h stands for step size).The intensity factors are defined by [13](90)KIL=limx→−l−2π−l−xσyyx,0,(91)KIR=limx→l+2πx−lσyyx,0,(92)KIIL=limx→−l−2π−l−xσxyx,0,(93)KIIR=limx→l+2πx−lσxyx,0,(94)KDL=limx→−l−2π−l−xDyx,0,(95)KDR=limx→l+2πx−lDyx,0.Substituting Equation (87) into Equations (90), (91), (92), (93), (94) and (95), we obtain(96)KIL=πlQ1∗H∗g1−1+Q2∗H∗g2−1+Q3∗H∗g3−1,KIR=−πlQ1∗H∗g11+Q2∗H∗g21+Q3∗H∗g31,KIIL=πlQ4∗H∗g1−1+Q5∗H∗g2−1+Q6∗H∗g3−1,KIIR=−πlQ4∗H∗g11+Q5∗H∗g21+Q6∗H∗g31,KDL=πlQ7∗H∗g1−1+Q8∗H∗g2−1+Q9∗H∗g3−1,KDR=−πlQ7∗H∗g11+Q8∗H∗g21+Q9∗H∗g31.According to the numerical process adopted in [28], the thermal stress and electric displacement fields can be calculated finally.


## 3.1. Temperature Field

By using the Fourier integral transform, the solution of temperature field can be expressed as(12)Tmx,y=12π∫−∞+∞A1mse−skmy+A2mseskmye−isxds,where km=kxm/kym, A1ms, and A2ms are the unknown functions m=1,2.To solve the problem easily, the density function is defined as follows:(13)Gx=∂∂xTmx,0+−Tmx,0−.According to Equation (11a), the dislocation density function should make the following single-valued conditions valid:(14)∫−llGtdt=0.Substituting formula (12) into Equations (11a), (11b), (11c) and (11d), we obtain(15)A11s+A21s−A12s−A22s=is∫−llGteistdt,(16)−sk1ky1A11s+sk1ky1A21s=−sk1ky2A12s+sk1ky2A21s,(17)A11se−sk1h1+A21sesk1h1=∫−∞+∞Taxeisxdx,(18)A12se−sk2h2+A22sesk2h2=∫−∞+∞Tbxeisxdx.Equations (16), (17) and (18) can be written in the form of a matrix as follows:(19)A⋅B=C,where(20)A=−sk1ky1sk1ky1sk2ky2−sk2ky211−1−1e−sk1h1esk1h10000e−sk2h2esk2h2,B=A11S,A21S,A12S,A22ST,C=0,I0,I1,I2T,I0=is∫−llGteistdt,I1=∫−llTaxeisxdx,I2=∫−llTbxeisxdx.Solving the linear Equation (19), we get(21)A11s=F1I0−F2I1+F3I2,A21=F4I1−F5I0−F6I2,A12s=F7I2+F8I0−F9I1,A22s=F10I1−F11I2−F12I0.The expressions ofFq and Δ, q=1,2,...12, are given in Appendix A. Obviously, the unknown functions A11s, A21s, A12s, and A22s depend on I0. Once I0 is determined, the temperature field could be obtained.By substituting Equation (12) into Equation (11c), we obtain(22)i2π∫−∞+∞sk1ky1F1−F5seit−xds∫−llGtdt=Q0x−d0x,where(23)d0x=12π∫−∞+∞sk1ky1F3I2−F2I1−F4I1−F6I2e−isxds.Furthermore, Equation (23) can be written as(24)1π∫0+∞Mssint−xsds∫−llGtdt=Q0x−d0x,where(25)Ms=sk1ky1F1−F5s.Through changing the order of the integration on the cracks’ surface, we obtain(26)−1π∫−llM0t−x+κ1x,tGtdt=Q0x−d0x,where(27)lims→∞Ms=M0,M0=−k1k2ky1ky2k1ky1+k2ky2,κ1x,t=∫0+∞Ms−M0sint−xsds.Equation (26) is a singular integral equation with Cauchy kernel. The density function Gt in Equation (26) has the square root-type singularity. Thus, we can express Gt as follows:(28)Gt=1l2−t2g0t,where g0t is a continuous function defined in the interval −l,l.The thermal flow intensity factors are defined by [13].(29)k0L=limx→−l−2π−x−lky∂Tmx,0∂y,(30)k0R=limx→l+2πx−lky∂Tmx,0∂y.According to Equation (26), the singular part of ky∂Tmx,0/∂y is(31)limx→−l−ky∂Tmx,0∂y=M0limx→−l−∫−llG−l−x−l,(32)limx→l+ky∂Tmx,0∂y=M0limx→l+∫−llGlx−l.Substituting Equation (28) into Equations (31) and (32), we obtain(33)k0L=M0πlg0−l,k0R=−M0πlg0l.


## 3.2. Thermal Stress and Electric Displacement Fields

Similar to the solution of linear equation, the solutions to Equations (7), (8) and (9) of equilibrium for displacements and electric potential consist of a homogeneous solution and a particular solution.Equations (7), (8) and (9) can be written into homogeneous forms as follows:(34)c11m∂2w∂x2+c44m∂2w∂y2+c13m+c44m∂2u∂x∂y+e31m+e15m∂2ϕ∂x∂y=0,(35)c44m∂2u∂x2+c33m∂2w∂y2+c13m+c44m∂2w∂x∂y+e15m∂2ϕ∂x2+e33m∂2ϕ∂y2=0,(36)e15m∂2u∂x2+e33m∂2u∂y2+e15m+e31m∂2w∂x∂y−ε11m∂2ϕ∂x2−ε33m∂2ϕ∂y2=0.By using the Fourier integral transform, the solutions of Equations (34), (35) and (36) can be expressed as [13](37)upx,y=12π∫−∞+∞∑k=16B1kmsesλkyfkse−isxds,(38)wpx,y=12π∫−∞+∞∑k=16B2kmsesλkyfkse−isxds,(39)ϕpx,y=12π∫−∞+∞∑k=16B3kmsesλkyfkse−isxds,where B1km, B2km, B3km, and fks, m=1,2, are unknown functions to be determined.The special solutions of Equations (7), (8) and (9) can be expressed as(40)uqx,y=12π∫−∞+∞1sB1amA1me−skmy+B1bmA2meskmye−isxds,(41)wqx,y=12π∫−∞+∞1sB2amA1me−skmy+B2bmA2meskmye−isxds,(42)ϕqx,y=12π∫−∞+∞1sB3amA1me−skmy+B3bmA2meskmye−isxds,where B1am, B2am, B3am, B1bm, B2bm, and B3bm satisfy the expression as follows:(43)Λ2km2−Λ1km+Λ0B1amB2amB3am=−iλ11m−is2sλ33mkms2spymkm,Λ2km2+Λ1km+Λ0B1bmB2bmB3bm=−iλ11m−is2sλ33mkms2spymkm,where(44)Λ0=−c11m000−c44m−e15m0−e15mε11m,Λ1=−is2s0c44me31m+e15mc44m00e31m+e15m00,Λ2=c44m000c33me33m0e33m−ε33m.Therefore, the solution of Equations (7), (8) and (9) can be written as(45)ux,y=12π∫−∞+∞∑k=16B1kmsesλkyfkse−isxds+12π∫−∞+∞1sB1amA1me−skmy+B1bmA2meskmye−isxds,(46)wx,y=12π∫−∞+∞∑k=16B2kmsesλkyfkse−isxds+12π∫−∞+∞1sB2amA1me−skmy+B2bmA2meskmye−isxds,(47)ϕx,y=12π∫−∞+∞∑k=16B3kmsesλkyfkse−isxds+12π∫−∞+∞1sB3amA1me−skmy+B3bmA2meskmye−isxds.Substituting Equations (45), (46) and (47) into Equations (2), (3), and (5), the thermal stress and electric displacement fields are obtained:(48)σyy=12π∫−∞+∞B1∗mf0se−isxds+12π∫−∞+∞B10∗mA0e−isxds,(49)σxy=12π∫−∞+∞B2∗mf0se−isxds+12π∫−∞+∞B20∗mA0e−isxds,(50)Dy=12π∫−∞+∞B3∗mf0se−isxds+12π∫−∞+∞B30∗mA0e−isxds,where(51)B1∗m=λ1c33mesλ1yB11m−ic31mesλ1yB21m+λ1e33mesλ1yB31mλ1c33mesλ1yB12m−ic32mesλ1yB22m+λ1e33mesλ1yB32mλ1c33mesλ1yB13m−ic33mesλ1yB23m+λ1e33mesλ1yB33mλ1c33mesλ1yB14m−ic34mesλ1yB24m+λ1e33mesλ1yB34mλ1c33mesλ1yB15m−ic35mesλ1yB25m+λ1e33mesλ1yB35mλ1c33mesλ1yB16m−ic36mesλ1yB26m+λ1e33mesλ1yB36mT,B2∗m=‐ic44mesλ1yB11m+λ1c44mesλ1yB21m−ie15mesλ1yB31m‐ic44mesλ1yB12m+λ1c44mesλ1yB22m−ie15mesλ1yB32m‐ic44mesλ1yB13m+λ1c44mesλ1yB23m−ie15mesλ1yB33m‐ic44mesλ1yB14m+λ1c44mesλ1yB24m−ie15mesλ1yB34m‐ic44mesλ1yB15m+λ1c44mesλ1yB25m−ie15mesλ1yB35m‐ic44mesλ1yB16m+λ1c44mesλ1yB26m−ie15mesλ1yB36mT,B3∗m=λ1e33mesλ1yB11m−ie31mesλ1yB21m−λ1ε33mesλ1yB31mλ1e33mesλ1yB12m−ie31mesλ1yB22m−λ1ε33mesλ1yB32mλ1e33mesλ1yB13m−ie31mesλ1yB23m−λ1ε33mesλ1yB33mλ1e33mesλ1yB14m−ie31mesλ1yB24m−λ1ε33mesλ1yB34mλ1e33mesλ1yB15m−ie31mesλ1yB25m−λ1ε33mesλ1yB35mλ1e33mesλ1yB16m−ie31mesλ1yB26m−λ1ε33mesλ1yB36mT,B10∗m=−ic13mB2am+skmc33msB1am+λ11me−skmyskmc33msB1bm−ic13mB2bm−λ11me−skmyT,B20∗m=−skmc44msB2am+ic44mB1am+ie15mB3ame−skmyskmc44msB2bm−ic44mB1bm−ie15mB3bme−skmyT,B30∗m=−skme33msB1am+skmε33msB3am+Pyme−skmyskme33msB1bm−skmε33msB3bm+Pyme−skmyT,f0=f1f2f3f4f5f6T,A0=A1mA2mT.In order to solve the thermal stress and electric displacement fields, the dislocation density functions are introduced:(52)d1x=ddxux,0+−ux,0−,(53)d2x=ddxwx,0+−wx,0−,(54)d3x=ddxϕx,0+−ϕx,0−.Substituting Equations (52), (53) and (54) into Equations (10a), (10b), (10c), (10d), (10e), (10f), (10g), (10h), (10i), (10j), (10k) and (10l), we obtain(55)∑k=16B1k1fk=H1I1∗−H2I2∗−H3I3∗+a1∗H,(56)∑k=16B1k2fk=H4I1∗−H5I2∗−H6I3∗+a2∗H,(57)∑k=16B2k1fk=H7I1∗−H8I2∗+H9I3∗+a3∗H,(58)∑k=16B2k2fk=H10I1∗−H11I2∗−H12I3∗+a4∗H,(59)∑k=16B3k1fk=H13I1∗−H14I2∗−H15I3∗+a5∗H,(60)∑k=16B3k2fk=H16I1∗−H17I2∗−H18I3∗+a6∗H,where(61)I1∗=is∫−lld1teistdt,(62)I2∗=is∫−lld2teistdt,(63)I3∗=is∫−lld3teistdt.The expressions ofH, ap∗, and Hq, p=1,2,⋯,6, q=1,2,⋯,18, are given in Appendix B. Similar to the process of solving the temperature field, once I1∗, I2∗, and I3∗ are determined, the thermal stress and electric displacement fields could be obtained.Substituting Equations (55), (56), (57), (58), (59) and (60) into Equations (48), (49) and (50), we get(64)σyyx,0−=limy→0−12π∫−∞+∞∑k=16sesλkyH−ic132H10+λkc332H4+λke332H16I1∗+−ic132H11−λkc332H5−λke332H17I2∗+ic132H12−λkc332H6−λke332H18I3∗e−isxds+12π∫−∞+∞−ic132B2aA12+B2bA22−λ112A12+A22+sk2e332s⋅−B3aA12+B3bA22+sk2c332s−B1aA12+B1bA22e−isxds=0,(65)σxyx,0−=limy→0−12π∫−∞+∞∑k=16sesλkyHλkc442H10−iλkc442H4−ie152H16I1∗+ic442H5−λkc442H11+ie152H17I2∗+ic442H6−λkc442H12+ie152H18I3∗e−isxds+12π∫−∞+∞−ic442B1aA12+B1bA22+sk2c442s−B2aA12+B2bA22−ie152−B3aA12+B3bA22e−isxds=0,(66)Dyx,0−=limy→0−12π∫−∞+∞∑k=16sesλkyHλke332H4−ie312H10−λkε332H16I1∗+ie312H11−λke332H5+λkε332H17I2∗+ie312H12−λke332H6+λkε332H18I3∗e−isxds+12π∫−∞+∞−ie312B2aA12+B2bA22−py2A12+A22−sk2ε332s⋅−B3aA12+B3bA22+sk2e332s−B1aA12+B1bA22e−isxds=0.Equations (64), (65) and (66) are divided into two parts, as follows:(67)σyyx,0−=σyy1x,0−+σyy2x,0−,σxyx,0−=σxy1x,0−+σxy2x,0−,Dyx,0−=Dy1x,0−+Dy2x,0−,where(68)σyy1x,0−=limy→0−12π∫−∞+∞∑k=16sesλkyH−ic132H10+λkc332H4+λke332H16I1∗+−ic132H11−λkc332H5−λke332H17I2∗+ic132H12−λkc332H6−λke332H18I3∗e−isxds,(69)σxy1x,0−=limy→0−12π∫−∞+∞∑k=16sesλkyHλkc442H10−iλkc442H4−ie152H16I1∗+ic442H5−λkc442H11+ie152H17I2∗+ic442H6−λkc442H12+ie152H18I3∗e−isxds,(70)Dy1x,0−=limy→0−12π∫−∞+∞∑k=16sesλkyHλke332H4−ie312H10−λkε332H16I1∗+ie312H11−λke332H5+λkε332H17I2∗+ie312H12−λke332H6+λkε332H18I3∗e−isxds,(71)σyy2x,0−=12π∫−∞+∞−ic132B2aA12+B2bA22−λ112A12+A22+sk2e332s⋅−B3aA12+B3bA22+sk2c332s−B1aA12+B1bA22e−isxds,(72)σxy2x,0−=12π∫−∞+∞−ic442B1aA12+B1bA22+sk2c442s−B2aA12+B2bA22−ie152−B3aA12+B3bA22e−isxds,(73)Dy2x,0−=12π∫−∞+∞−ie312B2aA12+B2bA22−py2A12+A22−sk2ε332s⋅−B3aA12+B3bA22+sk2e332s−B1aA12+B1bA22e−isxds.Obviously, it is easy to solve the solution of Equations (71), (72) and (73) according to Equations (39), (40), (41) and (42). Next, we need to transform Equations (68), (69) and (70) into singular integral equations of the first kind with Cauchy kernel.In another form, integral Equations (68), (69) and (70) will be(74)σyy1x,0−=limy→0−12π∫−∞+∞sesλyHQ1I1∗+Q2I2∗+Q3I3∗e−isxds,(75)σxy1x,0−=limy→0−12π∫−∞+∞sesλyHQ4I1∗+Q5I2∗+Q6I3∗e−isxds,(76)Dy1x,0−=limy→0−12π∫−∞+∞sesλyHQ7I1∗+Q8I2∗+Q9I3∗e−isxds.The expressions ofesλy and Qp, p=1,2,⋯,9, are given in Appendix C.Substituting Equations (61), (62) and (63) into Equations (74), (75) and (76), we get(77)12π∫−ll∫−∞+∞iesλyHQ1eist−xdsd1tdt+12π∫−ll∫−∞+∞iesλyHQ2eist−xdsd2tdt+12π∫−ll∫−∞+∞iesλyHQ3eist−xdsd3tdt=σ1,12π∫−ll∫−∞+∞iesλyHQ4eist−xdsd4tdt+12π∫−ll∫−∞+∞iesλyHQ5eist−xdsd5tdt+12π∫−ll∫−∞+∞iesλyHQ6eist−xdsd6tdt=σ2,12π∫−ll∫−∞+∞iesλyHQ7eist−xdsd7tdt+12π∫−ll∫−∞+∞iesλyHQ8eist−xdsd8tdt+12π∫−ll∫−∞+∞iesλyHQ9eist−xdsd9tdt=σ3,where(78)σ1=−12π∫−∞+∞−ic132B2aA12e−sk2y+B2bA22esk2y−λ112A12+A22+sk2e332s⋅−B3aA12e−sk2y+B3bA22esk2y+sk2c332s−B1aA12e−sk2y+B1bA22esk2ye−isxds,σ2=−12π∫−∞+∞−ic442B1aA12+B1bA22+sk2c442s−B2aA12+B2bA22−ie152−B3aA12+B3bA22e−isxds,σ3=−12π∫−∞+∞−ie312B2aA12e−sk2y+B2bA22esk2y−py2A12+A22−sk2ε332s⋅−B3aA12e−sk2y+B3bA22esk2y+sk2e332s−B1aA12e−sk2y+B1bA22esk2ye−isxds.On the cracks’ surface,(79)∫−llQ1∗H∗1t−x+κ2x,td1tdt+∫−llQ2∗H∗1t−x+κ3x,td2tdt+∫−llQ3∗H∗1t−x+κ4x,td3tdt=πσ1,(80)∫−llQ4∗H∗1t−x+κ5x,td1tdt+∫−llQ5∗H∗1t−x+κ6x,td2tdt+∫−llQ6∗H∗1t−x+κ7x,td3tdt=πσ2,(81)∫−llQ7∗H∗1t−x+κ8x,td1tdt+∫−llQ8∗H∗1t−x+κ9x,td2tdt+∫−llQ9∗H∗1t−x+κ10x,td3tdt=πσ3,where(82)κ2x,t=∫0∞Q1sH−Q1∗H∗sinst−xds,κ3x,t=∫0∞Q2sH−Q2∗H∗sinst−xds,κ4x,t=∫0∞Q3sH−Q3∗H∗sinst−xds,κ5x,t=∫0∞Q4sH−Q4∗H∗sinst−xds,κ6x,t=∫0∞Q5sH−Q5∗H∗sinst−xds,κ7x,t=∫0∞Q6sH−Q6∗H∗sinst−xds,κ8x,t=∫0∞Q7sH−Q7∗H∗sinst−xds,κ9x,t=∫0∞Q8sH−Q8∗H∗sinst−xds,κ10x,t=∫0∞Q9sH−Q9∗H∗sinst−xds.The expressions ofH∗ and Qp∗, p=1,2,⋯,10, are given in Appendix C. Ultimately, Equations (68), (69) and (70) are reduced to coupling the first kind Cauchy singular integral equations.In Equations (61), (62) and (63), the unknown function dkt, k=1,2,3 with the following conditions:(83)∫−lldktdt=0,k=1,2,3.Considering that the solution of the singular integral equation will use Chebychev polynomials which are defined in the interval [−1, 1], we introduce the following normalized quantities:(84)r=xl,s=tl,(85)dkx=dk~x,Lir,s=l⋅κix,t,k=1,2,3,i=1,⋯,10.Using Equations (84) and (85), integral Equations (79), (80) and (81) may be expressed as(86)∫−11Q1∗H∗1s−r+L2s,rd~1sds+∫−11Q2∗H∗1s−r+L3s,rd~2sds+∫−11Q3∗H∗1s−r+L4s,rd~3sds=πσ1,∫−11Q4∗H∗1s−r+L5s,rd~1sds+∫−11Q5∗H∗1s−r+L6s,rd~2sds+∫−11Q6∗H∗1s−r+L7s,rd~3sds=πσ2,∫−11Q7∗H∗1s−r+L8s,rd~1sds+∫−11Q8∗H∗1s−r+L9s,rd~2sds+∫−11Q9∗H∗1s−r+L10s,rd~3sds=πσ3.The crack-tip behavior can be characterized by square-root singularity weighed down by the oscillation. However, because the zone of oscillation is very small [25], the unknown function dk~s, k=1,2,3, can be expressed as follows:(87)dk~s=11−s2gks,k=1,2,3,where gks, k=1,2,3, are continuous functions defined in the interval −1,1.Using the Lobatto-Chebyshev collocation method, one can transform Equations (79), (80), (81) and (83) into a system of algebraic equations:(88)1N∑R=0NχRQ1∗H∗1sR−rq+L2sR,rqd~1sR+Q2∗H∗1sR−rq+L3sR,rqd~2sR+Q3∗H∗1sR−rq+L4sR,rqd~3sR=πσ1,1N∑R=0NχRQ4∗H∗1sR−rq+L5sR,rqd~1sR+Q5∗H∗1sR−rq+L6sR,rqd~2sR+Q6∗H∗1sR−rq+L7sR,rqd~3sR=πσ2,1N∑R=0NχRQ7∗H∗1sR−rq+L8sR,rqd~1sR+Q8∗H∗1sR−rq+L9sR,rqd~2sR+Q9∗H∗1sR−rq+L10sR,rqd~3sR=πσ3,∑R=0NχRgksR=0,k=1,2,3,q=1,2,⋯,N,where(89)χ0=χN=12,χ1=⋯=χN−1=1,rq=cos2q−1π2N,sR=cosRπN.N is the node number of quadrature formula; sR and rq are the zero points of the first and second kinds of Chebyshev polynomials. According to Theocaris and Ioakimidis [26] as well as Lu [27], the convergence rate of the Lobatto-Chebyshev collocation method is first order (i.e., oh, where h stands for step size).The intensity factors are defined by [13](90)KIL=limx→−l−2π−l−xσyyx,0,(91)KIR=limx→l+2πx−lσyyx,0,(92)KIIL=limx→−l−2π−l−xσxyx,0,(93)KIIR=limx→l+2πx−lσxyx,0,(94)KDL=limx→−l−2π−l−xDyx,0,(95)KDR=limx→l+2πx−lDyx,0.Substituting Equation (87) into Equations (90), (91), (92), (93), (94) and (95), we obtain(96)KIL=πlQ1∗H∗g1−1+Q2∗H∗g2−1+Q3∗H∗g3−1,KIR=−πlQ1∗H∗g11+Q2∗H∗g21+Q3∗H∗g31,KIIL=πlQ4∗H∗g1−1+Q5∗H∗g2−1+Q6∗H∗g3−1,KIIR=−πlQ4∗H∗g11+Q5∗H∗g21+Q6∗H∗g31,KDL=πlQ7∗H∗g1−1+Q8∗H∗g2−1+Q9∗H∗g3−1,KDR=−πlQ7∗H∗g11+Q8∗H∗g21+Q9∗H∗g31.According to the numerical process adopted in [28], the thermal stress and electric displacement fields can be calculated finally.


## 4. Numerical Process and Discussion

In what follows, we assume that the piezoelectric substrate is cadmium selenide piezoelectric ceramic. Its parameters are given as [24](97)c112=7.41×1010N/m2,e152=−0.138C/m2,c132=3.93×1010N/m2,ε112=0.825×10−10C/Vm,c442=1.32×1010N/m2,ε332=0.903×10−10C/Vm,e312=−0.16C/m2,λ112=0.621×106N/m2K,e332=0.347C/m2,λ332=0.551×106N/m2K,where superscript 2 represents the substrate material. Assume that the coefficients of thermal conductivity are km=kxm/kym=1, m=1,2.


### 4.1. Verification

We chose a special case to verify the foregoing derivation. We first assumed that the fine-grained piezoelectric coating and substrate are the same material. The problem in Figure1 then becomes that of a Griffith crack in a single material plate with thermal loading.Therefore, in Equation (26), M0=−k1k2ky1ky2/k1ky1+k2ky2 becomes(98)k1=k2=k=kxky,ky1=ky2=ky,M0→M∗=−kxky2.Thus, the thermal flow intensity factors (29) and (30) can be written as follows:(99)k0L=−kxky2πlg0−l,k0R=kxky2πlg0l.Because the crack is symmetrical about theY-axis, we get(100)k0=k0L=k0L.In Figure2, the effects of h1+h2/l on thermal flow intensity factor for different strip height are displayed (where K0∗=kxky/2πl). It is shown that with an increase in the thickness of the material plate, the size of the thermal strength factor decreases gradually and tends to stabilize. The results show that the thickness of the sheet has a great influence on the thermal strength factor, which is consistent with the conclusion of Wang and Noda [13].Figure 2
Effects ofh1+h2/l on thermal flow intensity factor for different strip height.


### 4.2. Thermal Intensity Factors of Fine-Grained Piezoelectric Coating/Substrate Structure

#### 4.2.1. Studies of the Effect of Elastic Modulus on Thermal Intensity Factor

Here, we select material 1, which is a fine-grained coating made of cadmium selenide piezoelectric ceramic; material 2 is normal cadmium selenide piezoelectric ceramic (where 1 and 2 represent the fine-grained coating and substrate material, respectively). This means that material 1 and material 2 have approximately the same piezoelectric and dielectric constants under suitable conditions, and the elastic modulus of material 1 is larger than that of material 2.Figures3–5 show the effects of the elastic modulus on the thermal intensity factor for different l/h1 (where h1/h2=1/100, h1=5mm, KI0=liml/h1→0KI, KII0=liml/h1→0KII, and KD0=liml/h1→0KD). These data indicate that with an increase in the crack length, the thermal stress intensity factor and the electric displacement intensity factor increase, and with the increase of the ratio of the elastic modulus between the fine-grained piezoelectric coating and substrate. The peaks of the thermal stress intensity factor and the electric displacement intensity factor increase correspondingly with the increase in the ratio of the elastic modulus between the fine-grained piezoelectric coating and substrate. Concurrently, when l/h1 increases from 2, the growth rates of KI and KII are greater than that of KD, and the interface fracture will eventually occur with a continuous increase in l/h1.Figure 3
Effects ofl/h1 on mode I thermal stress intensity factor for different c441/c442.Figure 4
Effects ofl/h1 on mode II thermal stress intensity factor for different c441/c442.Figure 5
Effects ofl/h1 on electric displacement intensity factor for different c441/c442.


#### 4.2.2. Studies of the Effect of the Thickness of Coating on Thermal Intensity Factor

Figures6–8 illustrate the effects of the thickness of the coating on thermal intensity factor for different l/h2 (where c441/c442=2.5, h2=100mm, KI0=liml/h2→0KI, KII0=liml/h2→0KII, and KD0=liml/h1→0KD). This shows that with an increase in crack length, the thermal stress intensity factor and the electric displacement intensity factor increase monotonically. The thermal stress intensity factor and the electric displacement intensity factor have a higher peak value when the coating thickness is small. However, when h1/h2=0.01, the influence of l/h2 on KD is not significant enough. The influence of l/h2 on KD is obviously significant with an increase of h1/h2.Figure 6
Effects ofl/h2 on mode I thermal stress intensity factor for different h1/h2.Figure 7
Effects ofl/h2 on mode II thermal stress intensity factor for different h1/h2.Figure 8
Effects ofl/h2 on electric displacement intensity factor for different h1/h2.


## 4.1. Verification

We chose a special case to verify the foregoing derivation. We first assumed that the fine-grained piezoelectric coating and substrate are the same material. The problem in Figure1 then becomes that of a Griffith crack in a single material plate with thermal loading.Therefore, in Equation (26), M0=−k1k2ky1ky2/k1ky1+k2ky2 becomes(98)k1=k2=k=kxky,ky1=ky2=ky,M0→M∗=−kxky2.Thus, the thermal flow intensity factors (29) and (30) can be written as follows:(99)k0L=−kxky2πlg0−l,k0R=kxky2πlg0l.Because the crack is symmetrical about theY-axis, we get(100)k0=k0L=k0L.In Figure2, the effects of h1+h2/l on thermal flow intensity factor for different strip height are displayed (where K0∗=kxky/2πl). It is shown that with an increase in the thickness of the material plate, the size of the thermal strength factor decreases gradually and tends to stabilize. The results show that the thickness of the sheet has a great influence on the thermal strength factor, which is consistent with the conclusion of Wang and Noda [13].Figure 2
Effects ofh1+h2/l on thermal flow intensity factor for different strip height.


## 4.2. Thermal Intensity Factors of Fine-Grained Piezoelectric Coating/Substrate Structure

### 4.2.1. Studies of the Effect of Elastic Modulus on Thermal Intensity Factor

Here, we select material 1, which is a fine-grained coating made of cadmium selenide piezoelectric ceramic; material 2 is normal cadmium selenide piezoelectric ceramic (where 1 and 2 represent the fine-grained coating and substrate material, respectively). This means that material 1 and material 2 have approximately the same piezoelectric and dielectric constants under suitable conditions, and the elastic modulus of material 1 is larger than that of material 2.Figures3–5 show the effects of the elastic modulus on the thermal intensity factor for different l/h1 (where h1/h2=1/100, h1=5mm, KI0=liml/h1→0KI, KII0=liml/h1→0KII, and KD0=liml/h1→0KD). These data indicate that with an increase in the crack length, the thermal stress intensity factor and the electric displacement intensity factor increase, and with the increase of the ratio of the elastic modulus between the fine-grained piezoelectric coating and substrate. The peaks of the thermal stress intensity factor and the electric displacement intensity factor increase correspondingly with the increase in the ratio of the elastic modulus between the fine-grained piezoelectric coating and substrate. Concurrently, when l/h1 increases from 2, the growth rates of KI and KII are greater than that of KD, and the interface fracture will eventually occur with a continuous increase in l/h1.Figure 3
Effects ofl/h1 on mode I thermal stress intensity factor for different c441/c442.Figure 4
Effects ofl/h1 on mode II thermal stress intensity factor for different c441/c442.Figure 5
Effects ofl/h1 on electric displacement intensity factor for different c441/c442.


### 4.2.2. Studies of the Effect of the Thickness of Coating on Thermal Intensity Factor

Figures6–8 illustrate the effects of the thickness of the coating on thermal intensity factor for different l/h2 (where c441/c442=2.5, h2=100mm, KI0=liml/h2→0KI, KII0=liml/h2→0KII, and KD0=liml/h1→0KD). This shows that with an increase in crack length, the thermal stress intensity factor and the electric displacement intensity factor increase monotonically. The thermal stress intensity factor and the electric displacement intensity factor have a higher peak value when the coating thickness is small. However, when h1/h2=0.01, the influence of l/h2 on KD is not significant enough. The influence of l/h2 on KD is obviously significant with an increase of h1/h2.Figure 6
Effects ofl/h2 on mode I thermal stress intensity factor for different h1/h2.Figure 7
Effects ofl/h2 on mode II thermal stress intensity factor for different h1/h2.Figure 8
Effects ofl/h2 on electric displacement intensity factor for different h1/h2.


## 4.2.1. Studies of the Effect of Elastic Modulus on Thermal Intensity Factor

Here, we select material 1, which is a fine-grained coating made of cadmium selenide piezoelectric ceramic; material 2 is normal cadmium selenide piezoelectric ceramic (where 1 and 2 represent the fine-grained coating and substrate material, respectively). This means that material 1 and material 2 have approximately the same piezoelectric and dielectric constants under suitable conditions, and the elastic modulus of material 1 is larger than that of material 2.Figures3–5 show the effects of the elastic modulus on the thermal intensity factor for different l/h1 (where h1/h2=1/100, h1=5mm, KI0=liml/h1→0KI, KII0=liml/h1→0KII, and KD0=liml/h1→0KD). These data indicate that with an increase in the crack length, the thermal stress intensity factor and the electric displacement intensity factor increase, and with the increase of the ratio of the elastic modulus between the fine-grained piezoelectric coating and substrate. The peaks of the thermal stress intensity factor and the electric displacement intensity factor increase correspondingly with the increase in the ratio of the elastic modulus between the fine-grained piezoelectric coating and substrate. Concurrently, when l/h1 increases from 2, the growth rates of KI and KII are greater than that of KD, and the interface fracture will eventually occur with a continuous increase in l/h1.Figure 3
Effects ofl/h1 on mode I thermal stress intensity factor for different c441/c442.Figure 4
Effects ofl/h1 on mode II thermal stress intensity factor for different c441/c442.Figure 5
Effects ofl/h1 on electric displacement intensity factor for different c441/c442.


## 4.2.2. Studies of the Effect of the Thickness of Coating on Thermal Intensity Factor

Figures6–8 illustrate the effects of the thickness of the coating on thermal intensity factor for different l/h2 (where c441/c442=2.5, h2=100mm, KI0=liml/h2→0KI, KII0=liml/h2→0KII, and KD0=liml/h1→0KD). This shows that with an increase in crack length, the thermal stress intensity factor and the electric displacement intensity factor increase monotonically. The thermal stress intensity factor and the electric displacement intensity factor have a higher peak value when the coating thickness is small. However, when h1/h2=0.01, the influence of l/h2 on KD is not significant enough. The influence of l/h2 on KD is obviously significant with an increase of h1/h2.Figure 6
Effects ofl/h2 on mode I thermal stress intensity factor for different h1/h2.Figure 7
Effects ofl/h2 on mode II thermal stress intensity factor for different h1/h2.Figure 8
Effects ofl/h2 on electric displacement intensity factor for different h1/h2.


## 5. Conclusions

The Griffith crack problem of fine-grained piezoelectric coating/substrates under steady-state thermal loading is studied. Based on the principle of integral transformation and superposition of equation solutions, the thermomechanical problem is transformed into a solution of thermal stress and electric displacement intensity factors. We then analyzed the effects of crack geometry, coating thickness, and material modulus of elasticity on the fracture behavior of the crack under thermal load. The results show that a larger elastic modulus and thinner coating thickness improve the safety of the coating/substrate structure. This reflects the advantages of fine-grained structure materials.



---

*Source: 4201591-2020-09-15.xml*