Small-Scale Domain Switching Near Sharp Piezoelectric Bi-Material Notches

RESEARCH Small-scale domain switching near sharp piezoelectric bimaterial notches Miroslav Hrstka .Michal Kotoul .Toma ´s ˇProfant .Marta Kianicova ´ Received: 29 July 2024 / Accepted: 4 November 2024 / Published online: 3 March 2025 ÓThe Author(s) 2025 Abstract Assuming a scenario of small-scale domain switching, the dimensions and configuration of the domain switching region preceding a clearly defined primarily monoclinic piezoelectric bi-material notch are determined by embracing the energetic switching principle and micromechanical domain switching framework proposed by Hwang et al. (Acta Metall Mater 43(5):2073–2084, 1995. https://doi.org/10.1016/09567151(94)00379-V) for a given set of materials, structure, and polarization alignment. The piezoelectric bi-material under consideration comprises piezoelectric ceramics PZT-5H and BaTiO 3 . The analysis of the asymptotic in-plane field around a bi-material sharp notch is conducted utilizing the extended Lekhnitskii–Eshelby– Stroh formalism (Ting in Anisotropic elasticity, Oxford University Press. 1996. https://doi.org/10.1093/oso/ 9780195074475.001.0001). Subsequently, the boundary value problem with theprescribed spontaneous strain and polarization within the switching domain is solved and their influence on the in-plane intensity of singularity at the tip of interface crack is computed. The effects of the initial poling direction on the resulting variation of the energy release rates are discussed. Keywords Small-scale domain switching Bimaterial piezoelectric sharp notch Expanded Lekhnitskii–Eshelby–Stroh formalism Two-state Hintegral List of symbols CE ij Material compliance matrix component [Pa] ^ C0E ij Reduced material compliance matrix component [Pa] DiElectric flux density [Cm -2 ] EcCoercive electric field [Vm -1 ] EiElectric intensity [Vm -1 ] GiEnergy release rate [Nm -1 ] Gtip i  Tip energy release rate [Nm -1 ] HiGeneralized stress intensity factor [Pa m1di] Htip i  Generalized tip stress intensity factor [Pa m1di] KI;KII Stress intensity factor (mechanical) [Pa ffiffiffiffi m p] Ktip I;Ktip II  Tip stress intensity factor (mechanical) [Pa ffiffiffiffi m p] KIV Electric stress intensity factor [C ffiffiffiffi m p] Ktip IV  Electric tip stress intensity factor [C ffiffiffiffi m p] M. Hrstka (&)M. Kotoul T. Profant Institute of Solid Mechanics, Mechatronics and Biomechanics, Brno University of Technology, Technicka ´ 2896/2, 616 69 Brno, Czech Republic e-mail: [email protected] M. Kotoul M. Kianicova ´ Faculty of Special Technology, Alexander Dubc ˇek University, Studentska ´2, 911 50 Trenc ˇı ´n, Slovak Republic 123 Int J Fract (2025) 249:1–30 https://doi.org/10.1007/s10704-024-00823-1(0123456789().,-volV)(0123456789().,-volV) SD ij Material stiffness matrix component [m 2 N -1 ] ^ S0D ij Reduced material stiffness matrix component [m 2 N -1 ] eij Piezoelectric constant [Cm -2 ] ^ e0 ij Reduced piezoelectric constant [Cm -2 ] gij Piezoelectric constant [C -1 m 2 ] ^ g0 ij Reduced piezoelectric constant [C -1 m 2 ] niUnit outward normal [–] PSRemanent polarization [Cm -2 ] r;hPolar coordinates [m], [rad] tiTractions [Nm -2 ], [Cm -2 ] TDGeneralized stress function related to electric displacement [Cm -1 ] T1;T2Generalized stress function related to stress [Nm -1 ] uiDisplacement [m], x1;x2Cartesian coordinates [m] x;yCartesian coordinates [m] a1;a2Poling direction [rad] br ij Dielectric non-permittivity [VmC -1 ] ^ b0r ij Reduced dielectric non-permittivity [VmC -1 ] diSingularity exponent [–] Deij Spontaneous strain change due to the domain switching [–] eij Strain tensor component [–] ccSpontaneous strain [–] liMaterial eigenvalue [–] rij Stress tensor component [Pa] rappl 2  Applied stress (loading) [Pa] /Electric potential [V] xe ij Dielectric permittivity [C(Vm) -1 ] ^ x0e ij Reduced dielectric permittivity [C(Vm) -1 ] x1;x2Notch face angles [rad] FEM Finite element method 1 Introduction Smart structures, like sensors and actuators, are composed of ferroelectric ceramics that possess inherent brittleness and susceptibility to cracking, thereby highlighting the crucial matter of their reliability. These devices typically showcase a multilayered structure frequently incorporating dissimilar material-bonded joints. The presence of a stress singularity at the interface’s edge arises due to a mismatch in material properties across interfaces. Numerous experimental studies have provided evidence indicating that significant singular stresses at joint vertices may trigger fractures and failures in the joints. For instance, actuators commonly experience cracking around the perimeters of internal electrodes. The singularity behaviour of dissimilar piezoelectric material bonded joints has been investigated in various scholarly works, such as those by Xu and Rajapakse (2000), Chue and Chen (2002), Ou and Wu (2003), Chen (2006), Ou and Chen (2004), Hwu (2014), Hirai et al. (2012), Abe et al. (2017), Hrstka et al. (2019), Luangarpa and Koguchi (2018) along with references therein. In order to comprehend the singularity behaviour of bonded joints, it is essential to explore two key parameters: the order of stress singularity and the intensity of singularity. The conservative integral method, grounded in Betti’s reciprocal principle as shown by Sinclair et al. (1984), Hwu (2010), Belov and Kirchner (1996), Stern et al. (1976), Banks-Sills and Sherer (2002), Profant et al. (2013) has been demonstrated as a highly effective tool for determining singularity intensities at a joint’s interface vertex in piezoelectric bi-material bonded joints (Hwu 2010,2014; Hirai et al. 2012; Abe et al. 2017;Hrstkaetal.2019; Luangarpa and Koguchi 2018; Sinclair et al. 1984; Belov and Kirchner 1996). However, this analytical approach is inherently linear and does not account for any nonlinear phenomena that may dominate in the vicinity of the vertex. It is widely recognized that ferroelectric ceramics exhibit pronounced nonlinear and hysteresis behaviours at high field strengths, with domain switching serving as the origin of the hysteresis loop. A number of researchers have analyzed domain switching as a potential source of toughening under electrical or mechanical loading in ferroelectrics, including studies byPark and Sun (1995), Hwang et al. (1995), Zhu and Yang (1997), Yang and Zhu (1998), Ricoeur and Kuna (2003), Zhang et al. (2006), Sheng and Landis (2007), Kuna (2010), Yang et al. (2001), Fulton and Gao (2001a,2001b), Tan et al. (2014), Rajapakse and Zeng (2001), Zeng and Rajapakse (2001), Zeng et al. (2003). It has been observed that domain switching significantly influences toughness variations. As a consequence of concentrated stress and electric fields near the crack tip, domain 123 2 Research reorientation occurs, inducing incompatible strain due to the constraint of the surrounding unswitched material, leading to stress redistribution near the crack tip. Thus, the apparent toughness of ferroelectrics fluctuates due to domain switching. Assuming that these non-elastic processes are confined to regions around the crack tip, which are small relative to relevant crack lengths, the principles of linear elastic fracture mechanics can be applied. Consequently, several basic domain switching models have been put forth in literature to predict crack tip process zones in ferroelectrics accurately, thereby anticipating potential toughening or softening effects (Hwang et al. 1995; Ricoeur and Kuna 2003; Rajapakse and Zeng 2001; Wang and Zhang 2007). The I-integral method combined with the phase field model successfully predicts crack-tip intensity factor variations due to domain switching, highlighting a toughening effect under displacement-controlled mechanical loading and a promotion of fracture under stress-controlled loading (Yu et al. 2016). Research indicates that under cyclic electric fields, a shielding effect occurs at the crack tip due to domain switching, making the crack more vulnerable during unloading due to residual fields (El Khatib et al. 2020). Additionally, flexoelectricity plays a crucial role in affecting the domain configuration near crack tips, with increasing flexoelectric coefficients leading to asymmetric domain configurations (Zhao and Soh 2018). Furthermore, micromechanical domain switching processes at crack tips are studied through numerical simulations, showing that the favoured domain orientation distribution depends on the position within the specimen and correlates with stresses and strains, impacting fracture toughness (Kozinov and Kuna 2019). These findings collectively emphasize the intricate interplay between small scale domain switching and crack behaviour in piezoelectric materials. It is important to highlight that a comprehensive investigation was conducted on large-scale domain switching utilizing the phase field model, which is not limited to small-scale domain switching. The impact on the fracture behaviour of ferroelectric materials was examined in several studies (Schrade et al. 2007; Wang and Kamlah 2009; Xu et al. 2010; Sluka et al. 2012; Abdollahi and Arias 2012; Su et al. 2015; Van Lich et al. 2015; Yu et al. 2016; Zhang and Su 2018). Nevertheless, it is crucial to emphasize that these analyses were specifically focused on crack-like defects in homogeneous materials. To the best of our knowledge there is no literature where the effect of the domain switching on the change in singularity intensity at the interface vertex in piezoelectric bi-material bonded joints or at the interface crack tip has been studied. This particular configuration is common in sensors and actuators within intelligent structures, and it will be the subject of discussion in the present study. 2 Problem formulation and solution methods The subsequent paragraph complies with the typical geometric characteristics of the bi-material notch under investigation as depicted in Fig. 1. Each wedge occupies the region 0 \h\x 1 or x 2 \h\0. The first part of the theory closely refers to the paper (Hrstka et al. 2019) but in order to keep the work comprehensive for the reader, the governing equations are stated in the following text. The reader is also encouraged to see Hrstka (2019) or Hwu and Ikeda (2008) for more detailed derivations. With respect to the material symmetry, piezoelectric ceramics with tetragonal symmetry are assumed to form bi-material bonded joints. However, because the crystallographic reference frame can be arbitrary rotated about x 3 axis with respect to global coordinate system in Fig. 1, the elasticity and piezoelectricity matrices possess the structure of monoclinic materials in the global coordinate system, hence the problem is formulated for the monoclinic materials with the symmetry plane at x 3 = 0. The constitutive laws for a linear elastic Fig. 1 Geometry of a bi-material notch characterized by two regions I and II. Notch faces are defined by angles x 1 and x 2 . Material interface is always considered at h= 0. Angles a 1 and a 2 denote poling direction of the materials I and II, respectively 123 Small-scale domain switching 3 piezoelectric material can be written in the matrix form as r D  ¼CEeT exe  e E  ; e E  ¼SDgT gbr  r D  ;ð1Þ where ris the stress tensor, eis the strain tensor (both written in a Voigt notation in the vector form), Eand D are the vectors of electric intensity and electric flux density,CEis the elastic stiffness tensor at constant electric field, CE ðÞis the stress/charge piezoelectric tensor, and xeis the dielectric permittivity tensor at constant strain, SDis the elastic compliance tensor at constant electric induction, gis the strain/voltage piezoelectric tensor, and bris the dielectric nonpermittivity tensor at constant stress (Hwu and Ikeda 2008; Hwu and Kuo 2009,2010). The material matrices are related through an inversion as CEeT exe  SDgT gbr  ¼I:ð2Þ The elasticity and piezoelectricity matrices for a monoclinic material poled in the x 1 x 2 -plane have the following structure: CE¼ CE 11 CE 12 CE 13 00CE 16 CE 12 CE 22 CE 23 00CE 26 CE 13 CE 23 CE 33 00CE 36 000CE 44 CE 45 0 000CE 45 CE 55 0 CE 16 CE 26 CE 36 00CE 66 2 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 5 ; r¼ r1 r2 r3 r4 r5 r6 8 > > > > > > > < > > > > > > > : 9 > > > > > > > = > > > > > > > ; ¼ r11 r22 r33 r23 r13 r12 8 > > > > > > > < > > > > > > > : 9 > > > > > > > = > > > > > > > ; ;e¼ e1 e2 e3 e4 e5 e6 8 > > > > > > > < > > > > > > > : 9 > > > > > > > = > > > > > > > ; ¼ e11 e22 e33 2e23 2e13 2e12 8 > > > > > > > < > > > > > > > : 9 > > > > > > > = > > > > > > > ; ; ð3Þ e¼ e11 e12 e13 00e16 e21 e22 e23 00e26 000e34 e35 0 2 6 43 7 5; xe¼ xe 11 xe 12 0 xe 12 xe 22 0 00xe 33 2 6 43 7 5;E¼ E1 E2 E3 8 > < > :9 > = > ;;D¼ D1 D2 D3 8 > < > :9 > = > ;: The structure of the matrices SD,gand bris similar to the structure of their inverse counterparts and is not stated here for the sake of brevity. The directional properties of the matrices depend on the poling axis, which can attain two limit configurations, either it coincides with x 1 -axis or with x 2 -axis. Between these states their structure corresponds to the above mentioned monoclinic one. The material constants in the general poling direction are obtained by using transformation relations SD¼K1  TS DK1;g¼XgK1;br¼Xb rX1; ð4Þ where S D,gand b rare material coefficient matrices with the poling direction aligned with certain principal direction (e.g. in which the properties were measured). The form of the transformation matrices Kand Xare given in Appendix A. External loads are assumed to be parallel to the plane defined by x 3 = 0 and plane strain deformations characterized by linear piezoelectricity are considered. These assumptions allow to decouple the inplane and anti-plane relations and to solve the in-plane and anti-plane problem separately. In the following only the in-plane problem is considered. The g-type constitutive relations for the in-plane problem can be expressed as (Hrstka 2019; Hrstka et al. 2019; Hwu and Ikeda 2008) e E  ¼bS0 Dbg0T bg0bb0 r "# r D  ;ð5Þ where r¼ r1 r2 r6 8 < :9 = ;¼ r11 r22 r12 8 < :9 = ;;e¼ e1 e2 e6 8 < :9 = ;¼ e11 e22 2e12 8 < :9 = ;; E¼E1 E2  ¼o/ ox1 o/ ox2 8 > < > :9 > = > ;;D¼D1 D2  ð6Þ are the stress vector, the strain vector, the electric field vector, the electric potential, and the electric displacement vector, respectively, 123 4 Research ^ S0 D¼ ^ S0D 11 ^ S0D 12 ^ S0D 16 ^ S0D 12 ^ S0D 22 ^ S0D 26 ^ S0D 16 ^ S0D 26 ^ S0D 66 2 6 43 7 5; 2 6 4 ^ g0¼ ^ g0 11 ^ g0 12 ^ g0 16 ^ g0 21 ^ g0 22 ^ g0 26  ; ^ b0 r¼ ^ b0r 11 ^ b0r 12 ^ b0r 12 ^ b0r 22 "## ð7Þ where ^ S0 D  ,^ g0 ðÞ, ^ b0 r  are the compliance matrix at constant induction, the piezoelectric strain/voltage matrix, and the dielectric impermeability matrix at constant stress, respectively, evaluated under the assumption of the generalized plane strain a short circuit e3= 0 and E 3 =0as ^ S0D ij ¼ ^ SD ij þ ^ g3i ^ g3j ^ br 33 ¼ ^ S0D ji ;^ g0 ij ¼^ gij  ^ br 3i ^ g3j ^ br 33 ; " ^ b0r ij ¼ ^ br ij  ^ br 3i ^ br 3j ^ br 33 ¼ ^ b0r ji ;#ð8Þ ^ SD ij ¼SD ij SD 3iSD 3j SD 33 ¼ ^ SD ji ;^ gij ¼gij gi3SD 3j SD 33;; ^ br ij ¼br ij gi3gj3 SD 33 ¼ ^ br ji;ð9Þ for i;j6¼ 3, where SD¼C1 EC1 EeTx1 reC1 E;xr¼ eC1 EeTþxe;g¼x1 reC1 E. The inverse form of the constitutive relations in Eq. (5) reads r D b C0 Ebe0T be0b x0 e  e E  ð10Þ where the in-plane elasticity and piezoelectricity matrices under the assumption of the generalized plane strain and short circuit have the following structure: ^ C0 E¼ ^ C0E 11 ^ C0E 12 ^ C0E 16 ^ C0E 12 ^ C0E 22 ^ C0E 26 ^ C0E 16 ^ C0E 26 ^ C0E 66 2 6 43 7 5; 2 6 4 ^ e0¼ ^ e0 11 ^ e0 12 ^ e0 16 ^ e0 21 ^ e0 22 ^ e0 26  ;^ x0 e¼ ^ x0e 11 ^ x0e 12 ^ x0e 21 ^ x0e 22  ð11Þ with e0 ij ¼eij xe 3ie3j xe 33 ;^ eij ¼eij e3iCE 3j CE 33;;^ e0 ij ¼^ eij  ^ xe 3i ^ e3j ^ xe 33 ; x0e ij ¼xe ij xe 3ixe 3j xe 33 ;^ xe ij ¼xe ij þe3ie3j CE 33;; ^ x0e ij ¼^ xe ij  ^ xe 3i ^ xe 3j ^ xe 33 ¼^ x0e ji ^ C0E ij ¼C0E ij ¼CE ij ;^ e0 ij ¼e0 ij ¼eij: ð12Þ Note that as the generalized plane strain and short circuit is considered and the poling direction lies in the x 1 x 2 -plane, zero elements with the index 3 in the permittivity matrix cause the equalities in Eq. (12) below. The in-plane problem of the stress singularity at the sharp notch composed of the two monoclinic piezoelectric materials is characterized by two generally complex exponents d 1 –1, d 2 –1 and one real exponent d 3 –1. The resulting displacements and stresses are obtained as the superposition of these particular singular contributions weighted by the generally complex amplitudes. The amplitudes are introduced in the similar manner as for a crack and referred to as the generalized stress intensity factors (GSIFs). The following vectors are introduced u¼ u1 u2 / 8 < :9 = ;;T¼ T1 T2 TD 8 < :9 = ;;ð13Þ where u 1 ,u 2 are the displacement components, /is the electric potential, T 1 ,T 2 and T D are the components of the generalized stress function vector of the mechanical and electrical quantities. The generalized stress function components T 1 ,T 2 and T D are related to the stresses and electric displacements by ri1 ¼Ti;2;ri2 ¼Ti;1;i¼1;2;D1¼TD;2;D2 ¼TD;1: ð14Þ The complex potential for expression of the stress singularity has the form Zd¼diag zd 1;zd 2;zd 3  ;zi¼x1þlix2;i¼1;2;3; ð15Þ where liare the material eigenvalues and they are obtained by solving the characteristic equation Eq. (A5) (see Appendix A). The coupled electromechanical field near the notch vertex is sought in the 123 Small-scale domain switching 5 form of the ansatz for unknown singularity exponents and corresponding eigenvectors as ur;hðÞ¼AZdhðÞvþAZdhðÞw; Tr;h ðÞ ¼LZdh ðÞ vþLZdh ðÞ w;ð16Þ where the complex function Zdand matrices Aand L are defined in Appendix A, the bar above the symbols denotes complex conjugate quantities. viand wiare eigenvectors pertinent to the singularity exponent eigenvalue problem of the considered sharp notch in Fig. 1. The unknown singularity exponents and eigenvectors are determined through the satisfaction of the boundary conditions at the notch tip. The material eigenvectors matrices A,Lin Eq. (16) are non-degenerate when the poling direction is perpendicular to the x 3 axis and semi-degenerate or degenerate otherwise. Only non-degenerate cases are considered thereinafter. Considering traction and charge free notch faces the following boundary conditions are imposed by: TIx1 ðÞ¼0; TII x2 ðÞ¼0:ð17Þ The displacement and traction continuity conditions are prescribed along the interface h¼0as uI0ðÞ¼uII 0ðÞ; TI0ðÞ¼TII 0ðÞ:ð18Þ By substituting vector functions kihðÞand gihðÞ from Eq. (16) into Eqs. (17) and (18), one gets the singularity exponent eigenvalue problem formed by twelve homogeneous algebraic equations, which can be written in the matrix form as XI 1XI 100 00XII 2XII 2 BI 0BI 0BII 0BII 0 IIII 2 6 6 6 43 7 7 7 5: LIvI LIwI LIIvII LIIwII 8 > < > :9 > = > ;¼0;ð19Þ where Xj¼LZd jxj  LðÞ 1;Xj¼LZd jxj L  1;j¼1;2ðÞ B0¼iAL1;B0¼iAL1 ð20Þ and 0denotes 3 93 zero matrix on the left-hand side and 12 91 zero vector on the right-hand side of Eq. (19). After some algebraic manipulations one receives the characteristic equation for the singularity exponent values in the form det KdðÞ½¼0;ð21Þ where K¼BI 0þBI 0YI 1BII 0þBII 0YII 2  IYII 2  1IYI 1  and Yj¼X1 jXj;j¼1;2ðÞ: ð22Þ For details see Hrstka et al. (2019) and Hrstka (2019). In the case of an in-plane problem and monoclinic piezoelectric materials, there are three generally complex roots of the Eq. (21). In a case of interfacial crack studied in this paper, there are two complex conjugate roots d 1, d 2 and one real root d 3 . In order to simplify the numerical algorithm and corresponding relations for the two-state integral, let us introduce the following angular functions: gihðÞ¼AZdihðÞviþAZdihðÞwi; kihðÞ¼LZdihðÞviþLZdihðÞwi;ð23Þ where i= 1, 2, 3 denotes the individual singularity exponent. Then, the displacement and stress function vectors are defined as their combinations as: ur;hðÞ¼H1rd1g1hðÞþH2rd2g2hðÞþH3rd3g3hðÞ; Tr;hðÞ¼H1rd1k1hðÞþH2rd2k2hðÞþH3rd3k3hðÞ; ð24Þ where H i are generalized stress intensity factors, rand hare polar coordinates, see Fig. 1. Let us introduce following functions: d ~ kihðÞ dx1¼LdiZdi1hðÞviþLdiZdi1hðÞwi;i¼1;2;3 d ~ kihðÞ dx2¼LdiZdi1hðÞlviþLdiZdi1hðÞlwi;i¼1;2;3; 2 6 6 6 43 7 7 7 5 ð25Þ where the parameter ~ kihðÞdenotes the function kihðÞ in Eq. (23) with function Zexpressed in the Cartesian coordinates x 1 and x 2 . Using Eq. (25), the asymptotic stresses and electrical displacements can be written as 123 6 Research r1r;hðÞ¼H1rd11d ~ k1hðÞ dx2H2rd21d ~ k2hðÞ dx2H3rd31d ~ k3hðÞ dx2 ; r2r;hðÞ¼H1rd11d ~ k1hðÞ dx1þH2rd21d ~ k2hðÞ dx1þH3rd31d ~ k3hðÞ dx1 ; 2 6 6 6 43 7 7 7 5 ð26Þ where r1¼ r11 r12 D1 8 > < > :9 > = > ;;r2¼ r21 r22 D2 8 > < > :9 > = > ;; l¼ l100 0l20 00l3 2 6 43 7 5; l¼ l100 0l20 00 l3 2 6 43 7 5: To determine generalized stress intensity factors in Eq. (24), a conservative line integral for the bimaterial notch developed from Betti’s reciprocal principle mentioned in the Introduction is applied. Contrary to the J-integral, the path-independence of the line integral developed from Betti’s reciprocal principle is also preserved for multi-material stress concentrators if the body forces and charges are neglected. This line integral referred to as the Hintegral can be written for a bimaterial notch characterized by angles x 1 and x 2, see Fig. 1,as Hu;^ ui ðÞ¼ Zx1 x2 uT^ ti^ uT it  rdh:ð27Þ The vectors ^ uT iand ^ tiare the auxiliary solutions to the displacements, tractions, electric potential and the charge and correspond to the exponent ^ di¼di. The auxiliary solutions are defined as ^ uir;hðÞ¼rdi^ gihðÞ; ^ tir;hðÞ¼ 1 r o ^ Tir;hðÞ oh¼rdi1^ k0 ihðÞ;i¼1;2;3ðÞ 2 6 43 7 5 ð28Þ where ^ gihðÞ¼AZdihðÞ ^ viþAZdihðÞ ^ wi; ^ ki0hðÞ¼LZ dihðÞ  0^ viþL ZdihðÞ  0^ wi; 2 43 5 where (.)0denotes the differentiation with respect to h. Observe that the vectors uand tin Eq. (27) represent either the regular asymptotic or a full field solution obtained numerically e.g. by FEM. In the first case, the vector uis given by Eq. (16) 1 and the vector tis given by the derivative of Eq. (16) 2 with respect to h similarly as ^ tiin Eq. (28) 2 . The full-field solution reduces to the asymptotic solution if the integration contour shrinks to the notch tip. Since the regular and corresponding auxiliary solutions are orthogonal with respect to the ‘‘scalar product’’ defined by the integral (27), i.e. Hr djgjhðÞ;rdi^ gihðÞ  ¼const 6¼ 0 for i¼j; 0 for i6¼ j; ð29Þ an important result for the GSIFs evaluation follows as Hi¼ HuFEM;rdi c ^ gihðÞ  Hr digihðÞ;rdi ^ gihðÞ  i¼1;2;3ðÞ; 2 43 5ð30Þ where HuFEM;rdi c ^ gihðÞ  ¼Zx1 x2 uFEM  Trdi1 c ^ ki0hðÞ  2 4þrdi c ^ gi ThðÞtFEMrcdhð31Þ with uFEM  and tFEM  standing for the FEM approximation to the vectors uand t, respectively, and with r c denoting the radius of the circular path remote from the notch singularity. Elements of the vector tFEM  along the integrating contour have to be computed from the stresses using the Cauchy formula t i =r ij n j , in the matrix form written as tFEM ¼rFEMn  , where rFEM  is the twodimensional generalized stress tensor and nis the outer normal to the domain enclosed by the circular integrating path of the radius r c defined as 123 Small-scale domain switching 7 rFEM ¼ rFEM 11 rFEM 12 rFEM 21 rFEM 22 DFEM 1DFEM 2 2 6 43 7 5;n¼cos hðÞ sin hðÞ  : 2 6 43 7 5 ð32Þ In order to predict the domain switching zone, the energy-based criterion proposed in Hwang et al. (1995) is applied rijDeij þEiDPi2PsEc;ð33Þ where Deij and DPiare the changes in the spontaneous strain and the spontaneous polarization during switching, respectively. Psis the magnitude of the spontaneous polarization (remanent polarization); and Ecthe coercive electric field. As a first approximation it is assumed that stresses rij and electrical fields Eiremain unchanged during switching. That means that the linear asymptotic field in Eq. (16) dominates at the notch tip and the loading is thus controlled by the GSIFs H i . This case is referred to as small scale switching when the size of the switching zone is much smaller than other specimen dimensions. Note that the left-hand side of Eq. (33) represents the specific work dissipated during switching. The threshold value on the right-hand side of Eq. (33) represents approximately half of the area of a polarization hysteresis. Due to electric and/or stress loading the spontaneous polarization of a domain near the notch tip can rotate by 180°,?90°or -90°. Considering that a ferroelectric domain forms an angle awith x 1 axis the changes in spontaneous strain Deij and polarization DPifor 90°domain switchings can be expressed in matrix form as (Hwang et al. 1995; Zeng and Rajapakse 2001): De¼ De1 De2 De6 8 < :9 = ;¼ De11 De22 2De12 8 < :9 = ;¼cscos 2aðÞ cos 2aðÞ 2 sin 2aðÞ 8 < :9 = ;; 2 43 5 ð34Þ DP¼ffiffiffi2 pPs cos a3p 4  sin a3p 4  2 6 6 43 7 7 5;ð35Þ where csdenotes the spontaneous strain associated with 90°domain switching and 3p 4and 3p 4in Eq. (35) correspond to clockwise and counterclockwise 90° switching, respectively. Note that for 180°Deij ¼0 and DP¼2Pscosasina½ T. A variation of the local stress/electric fields near the vertex induced by domain switching leads to the additional GSIFs DH i which are evaluated again using Betti’s reciprocal principle. To start with, consider initial generalized stress tensor rrur ðÞ;Drur ðÞ fg T, and initial generalized strain tensor erur ðÞ;Erur ðÞ fg Tin the body. Generalized stress rrur ðÞ;Drur ðÞ fg Tdeveloped due to changes in spontaneous strain Deand polarization DPsatisfies equilibrium orr Dr  ¼0inX;ð36Þ where the matrix differential operator is defined as o¼ o ox1 0o ox2 00 0o ox2 o ox1 00 00 0o ox1 o ox2 2 6 6 6 6 6 4 3 7 7 7 7 7 5 and Xdenotes the whole bimaterial body with the body forces and free charges absent. Further, rrur ðÞ;Drur ðÞ fg Tsatisfies the constitutive law in X rrDrDP fg ¼^ C0 E ^ e0T^ e0^ x0 e  erDe Er  ; ð37Þ or rewritten as rrDr fg ¼^ C0 E ^ e0T^ e0^ x0 e  er Er   ^ C0 EDe ^ e0DeDP  ;ð38Þ with ^ C0 E;^ e0;^ x0 e  having the same structure as given in Eq. (11) and De=DP= 0 outside the switching zone XSZ :The boundary and continuity conditions are as follows: tr¼n½ rr Dr  ¼ n10n200 0n2n100 000n1n2 2 43 5rr Dr  ¼0; on oX; ð39Þ where n 1 ,n 2 are components of the unit outer normal vector. Across the boundary of the switching zone 123 8 Research oXSZ the tractions and the displacements are continuous tr ½½¼ n½ rr DrDP  ¼0; ur ½½¼ u1r u2r /r 8 < :9 = ; 2 43 5 2 43 5¼0onoXSZ ;ð40Þ where the double brackets :½½denote a jump of the quantity inside. Finally, the tractions and displacements are required to be continuous along the interface uI r¼ uI 1r uI 2r /I r 8 > < > :9 > = > ;¼ uII 1r uII 2r /II r 8 > < > :9 > = > ;¼uII r;tI r¼tII rð41Þ The generalized strain tensor erur ðÞ;Erur ðÞ fg T satisfies er Er  ¼oTur:ð42Þ For arbitrary virtual displacement duthe weak form of the problem given by the Eqs. (36)–(41) reads ZX ^ C0 Eer^ e0TEr  dedXþZX ^ e0erþ^ x0 eEr  dEdX ¼Z XSZ ^ C0 EDedeþ^ e0DeþDPðÞdE  dX: 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 ð43Þ Equation (43) forms a basis of mixed finite element formulation. Considering the initial generalized stress and strain tensor, the reciprocal theorem for two arbitrary admissible fields uand ^ u(the auxiliary field in our case) reads Z D¼DIþDII ruðÞ;DuðÞ fg rrur ðÞ;Drur ðÞ fg ðÞ e^ uðÞ E^ uðÞ  2 6 4 r^ uðÞ;D^ uðÞ fg euðÞ EuðÞ  erur ðÞ Erur ðÞ  dS¼0;# ð44Þ or Z D¼DIþDII ruðÞ;DuðÞ fg oT^ ur^ uðÞ;D^ uðÞ fg oTurrur ðÞ; f  2 6 4 Drur ðÞg e^ uðÞ E^ uðÞ  þr^ uðÞ;D^ uðÞ fg erur ðÞ ð45ÞErur ðÞ  dS¼0; ð45Þ where the domain Dis any subset of the original domain obtained by excluding the crack tip. The boundary of Dis made up of arbitrary contours C 1 and C 2 circumventing the crack tip and connecting the traction free crack faces. C 2 is a remote path whereas C 1 is a path very close to the crack tip. Clearly, the integration over the subdomains DI  and DII  is performed using pertinent generalized stress and strain fields in the respective subdomains. Observe that it holds r^ uðÞ;D^ uðÞ fg ¼e^ uðÞ E^ uðÞ  T ^ C0 E ^ e0T^ e0^ x0 e  T "# ð46Þ and the last term in Eq. (45) can be written using Eq. (38)as r^ uðÞ;D^ uðÞ fg erur ðÞ Erur ðÞ  ¼erur ðÞ;Erur ðÞ fg ^ C0 E ^ e0T ^ e0^ x0 e "# e^ uðÞ E^ uðÞ  ¼rrur ðÞ;Drur ðÞ fg e^ uðÞ E^ uðÞ  þ^ C0 EDe^ e0DeDP  e^ uðÞ E^ uðÞ  : 2 6 6 6 6 4 3 7 7 7 7 5 ð47Þ Applying the divergence theorem to the first two terms in Eq. (45) one obtains Z C¼C1[C2 ruðÞ;DuðÞ fg n½ T^ ur^ uðÞ;D^ uðÞ fg n½ Tu  ds ¼ Z D¼DIþDII ^ C0 EDe^ e0DeDP  e^ uðÞ E^ uðÞ  dS: 2 6 6 6 6 6 4 3 7 7 7 7 7 5 ð48Þ Reversing the flow direction of the contour C 1 and using FE approximation to the field ualong the contour C 2 we obtain 123 Small-scale domain switching 9 Fig. 6 Domain switching zones for an interface crack loaded by rappl 2¼5 MPa for various orientations of the initial poling direction with respect to the interface; on the left-hand side the detailed view is shown 123 16 Research Fig. 7 Domain switching zones for a bi-material notch with the angle x1¼150loaded by rappl 2¼5 MPa for various orientations of the initial poling direction with respect to the interface; on the left-hand side the detailed view is shown 123 Small-scale domain switching 17 Fig. 8 The displacements, stress components, electric displacement components and electric potential along the circular path r=1 mm for a PZT-5H/BaTiO 3 interface crack, a¼0;loading rappl 2¼5 MPa 123 18 Research Fig. 9 The displacements, stress components, electric displacement components and electric potential along the circular path r=1 mm for a PZT-5H/BaTiO 3 interface crack, a¼70;loading rappl 2¼5 MPa 123 Small-scale domain switching 19 123 20 Research accompanied by a large increase of the electric displacement component D 2 ahead of the crack tip. Another point is also interesting—whilst before switching the stress r 11 is discontinuous across the interface, after switching it gets closer for a290 130and the jump is not very significant (the stress is almost continuous). Figure 13 shows the conventional stress intensity factors calculated using Eq. (55) as functions of the initial poling direction before switching K II ,K I ,K IV , and after switching Ktip II ;Ktip I;Ktip IV . It is seen that while before switching the stress intensity factors depend only weakly on the initial poling direction; after switching they strongly vary with the initial poling direction. Mechanical SIFs Ktip II ;Ktip Iare significantly reduced for the initial poling direction values a240;120140;while the electrical SIF Ktip IV is significantly amplified within this interval of a, which is in accord with the asymptotic fields shown in Figs.8,9,10,11 and 12. The ratios Ktip II =KII;Ktip I=KI;Ktip IV =KIV are plotted in Fig. 14, which better represent the influence of switching domain independently on the applied load under assumption of small-scale switching. It is known, see e.g. Hwu and Ikeda (2008), that in case of interface crack the mechanical load alone can induce nonzero value of the electrical intensity factor KIV for in-plane poling. The presented results show that this effect can be distinctly intensified due to poling switching. As already mentioned in the previous section, it is more convenient to compare the energy release rate values to assess the overall influence of the domain switching zone in case of interface crack between two dissimilar piezoelectric materials. The energy release rate before switching, G, and after switching, Gtip, was calculated using Eq. (59) as function of the orientation of the initial poling direction. The total energy release rate, the mechanical energy release rate (labelled in Fig. 15 as ‘‘without K2 IV ’’) and the pure electrical energy release rate were plotted both before and after switching. It is well known that linear piezoelectricity always predicts that the electrical contribution to the energy release rate is negative. Here, the pure electrical energy release rate before switching is negligible in comparison to mechanical one, as it can be seen in Fig. 15, because only mechanical load of the interface crack is considered and the induced electrical intensity factor KIV is low, see red line in Fig. 13c. Nevertheless, the switching model predicts that the electrical contribution to the energy release rate after switching becomes dominant and makes the total energy release rate negative for the initial poling direction values a240;150, see Fig. 15. Figure 16 then shows the ratios Gtip=Gof the total energy release rate and the mechanical energy release rate only plotted as functions of the initial poling direction. The strong negative electrical contribution to the energy release rate due to switching can be explained using the crack closure integral expression G¼lim Da!0 1 2aZDa 0 Du1DasðÞr12 sðÞþDu2DasðÞ½ r22 sðÞþD2sðÞD/DasðÞds ð70Þ With Dui;D/standing for the jumps of displacements and of the electric potential over the crack faces, and results shown in Figs. 9,10 and 11. In these figures one can observe that poling switching induces a remarkable negative increase of potential difference over the crack faces accompanied by a large increase of the electric displacement component D 2 ahead of the crack tip. As a results, the last term in in the integrand of the crack closure integral obtains large negative values for the initial poling direction values a240;150. 4 Conclusions For the first time a theoretical model was proposed to analyse small-scale domain switching near an interface crack in PZT-5H/BaTiO 3 bi-material and its impact upon the energy release rate under pure mechanical loading. The simple energy-based criterion proposed in Hwang et al. (1995) was applied for prediction of the domain switching zone ahead of the interface crack and the bi-material notch with the bFig. 10 The displacements, stress components, electric displacement components and electric potential along the circular path r= 1 mm for a PZT-5H/BaTiO 3 interface crack, a¼90; loading rappl 2¼5 MPa 123 Small-scale domain switching 21 Fig. 11 The displacements, stress components, electric displacement components and electric potential along the circular path r=1 mm for a PZT-5H/BaTiO 3 interface crack, a¼120;loading rappl 2¼5 MPa 123 22 Research Fig. 12 The displacements, stress components, electric displacement components and electric potential along the circular path r=1 mm for a PZT-5H/BaTiO 3 interface crack, a¼180;loading rappl 2¼5 MPa 123 Small-scale domain switching 23 angle x1¼150:As boundary conditions only traction and charge free faces were considered. It was shown that the orientation of the initial poling has a marking influence on the size and the shape of the switching zone. The effects of the switching zone upon the local electromechanical field was calculated only for the special case of an interface crack under pure mechanical loading. Consequently, the energy release rates before switching, G, and after switching, Gtipthe initial poling direction values a240;150due to dominant electrical contribution to the energy release rate after switching. To conclude, a robust computational procedure based upon the expanded Lekhnitskii–Eshelby–Stroh formalism was developed to evaluate the influence of domain switching zone on the in-plane asymptotic field of the bi-material sharp notch and the interface crack in a generally monoclinic piezoelectric bimaterial. Betti’s reciprocal principle was applied to calculate the general stress intensity factors and also to capture the influence the domain switching zone, thus avoiding a derivation of weight functions. The supplementary data files can be downloaded from the Zenodo online repository, see below. As a next step, formulation of fracture criteria and experimental verification would be desirable. Also application of the semi-permeable crack boundary conditions and socalled energetically consistent boundary conditions is required to analyse the role of electric crack face boundary conditions on domain switching in case of interface cracks and bi-material notches. 5 Supplementary materials Regarding the computational procedures see Hrstka, M. (2024). Data for ‘‘Small-Scale Domain Switching Near Sharp Piezoelectric Bi-Material Notches’’ Fig. 13 Stress intensity factors as functions of the orientation of the initial poling direction calculated from Eq. (55) 123 24 Research (1.0.0). Zenodo. https://doi.org/https://doi.org/10. 5281/zenodo.12570340 Acknowledgements The authors acknowledge the supports by Scientific Grant Agency of the Ministry of Education, Science, Research and Sport of the Slovak Republic and Slovak Academy of Sciences via the project registered under the number VEGA-1/0327/21 and Horizon Europe via the project HORIZON-WIDERA-2021-ACCESS-03, SEP-210806308 and by the project ‘‘Mechanical Engineering of Biological and Bioinspired Systems’’, funded as project No. CZ.02.01.01/00/ 22_008/0004634 by Programme Johannes Amos Commenius, call Excellent Research and by the project GACR 22-14387J Design and manufacturing of 4D metamaterials based on printed structures with embedded elements of smart materials. Computational resources were provided by the e-INFRA CZ project (ID:90254), supported by the Ministry of Education, Youth and Sports of the Czech Republic. Author contributions Miroslav Hrstka performed the numerical simulations, contributed with the numerical methods to the manuscript and discussion. Michal Kotoul wrote the theoretical part of the manuscript and discussion. Toma ´s ˇProfant performed the numerical simulations. Marta Kianicova ´participated at the manuscript writing and organized financial support and administration. All authors reviewed the manuscript. Data availability Regarding the computational procedures see Hrstka, M. (2024). Data for ‘‘Small-Scale Domain Switching Near Sharp Piezoelectric Bi-Material Notches’’ (1.0.0). Zenodo. https://doi.org/https://doi.org/10.5281/zenodo. 12570340. Declarations Competing interests The authors declare no competing interests. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your Fig. 14 Ratios Ktip II =KII ;Ktip I=KI;Ktip IV =KIV plotted as functions of the initial poling direction 123 Small-scale domain switching 25