# Cell-Based Smoothed Finite Element Method-Virtual Crack Closure Technique for a Piezoelectric Material of Crack

Authors: Li Ming Zhou; Guang Wei Meng; Feng Li; Hui Wang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

In order to improve the accuracy and efficiency of solving fracture parameters of piezoelectric materials, a piezoelectric element, tailored for the virtual crack closure technique (VCCT), was used to study piezoelectric materials containing a crack. Recently, the cell-based smoothed finite element method (CSFEM) and VCCT have been used to simulate the fracture mechanics of piezoelectric materials. A center cracked piezoelectric materials with different material properties, crack length, mesh, and smoothing subcells at various strain energy release rates are discussed and compared with finite element method-virtual crack closure technique (FEM-VCCT). Numerical examples show that CSFEM-VCCT gives an improved simulation compared to FEM-VCCT, which generally simulates materials as too stiff with lower accuracy and efficiency. Due to its simplicity, the VCCT piezoelectric element demonstrated in this study could be a potential tool for engineers to practice piezoelectric fracture analysis. CSFEM-VCCT is an efficient numerical method for fracture analysis of piezoelectric materials.


---

## Body

## 1. Introduction

Piezoelectric materials have been widely used in high technology fields due to their attractive electromechanical coupling characteristics. Piezoelectric materials are typically brittle materials. Therefore, pores and cracks often arise in their manufacture or application process due to the electromechanical joint effect. The main cause of cracks is material failure. Solving the fracture parameters of piezoelectric materials accurately will have significant impact on their applications and may lead to device performance improvements.Pak [1], Sosa [2], Suo et al. [3], Wang [4], and Zhang and Hack [5] began research on the fracture mechanics of piezoelectric materials in the early 1990s and have since become the focus of attention in this field [6–9]. In the 20 years to date, wider research has been conducted by both domestic and foreign researchers, with a remarkable progress as a result. The theoretical framework of the fracture mechanics of piezoelectric materials has been established. However, the theoretical model applies only to simple questions and in order to solve more complex problems, one still has to resort to numerical methods. The first significant numerical attempt using finite element implementation for piezoelectric phenomenon was a piezoelectric vibration analysis proposed by Allik and Hughes [10].Until now, displacement finite element method (FEM) models have been used mostly for engineering problems. However, it is well known that FEM produces overestimations of the stiffness matrix [11, 12]. As a consequence, the solution is always smaller than the real result. Additionally, since mapping and coordinate transforms are involved in the FEM, elements are not allowed to be of arbitrary shape. In the effort of overcoming these problems, Liu et al. proposed for the first time a cell-based smoothed finite element method (CSFEM) by combining the existing FEM technology with the strain smoothing technique of mesh-free methods [13]. No derivative of the shape functions is involved in computing the field gradients to form the stiffness matrix. Correspondingly, the element shape in CSFEM can be of arbitrary shape. In CSFEM, the strain in an element is modified by smoothing the compatible strains over quadrilateral smoothing domains, which gives important softening effects. CSFEM can improve the accuracy and convergence rate of the FEM-Q4 model using the same mesh. The SFEM was extended to various problems such as shells [14], piezoelectric material [15], fracture mechanics [16], heat transfer [17], and structural acoustics [18] among others. CS-FEM has been combined with the extended FEM to address problems involving discontinuities.Rabczuk et al. [19] presented an extension of the phantom-node method by allowing crack tips to be placed within a finite element. Thereby, the crack growth in the phantom-node method became almost independent of the finite element mesh. Wu et al. [20] applied the NMM to investigate the cracking behavior of a sedimentary rock under dynamic loading. By incorporating the NMM with the cracking processes, crack initiation, propagation, and coalescence were successfully modeled. The element free Galerkin method (EFGM) [21, 22], developed by Belytschko et al., has a unique feature in solving the problems of crack growth. The notable feature of this method is that there is no mesh required in establishing a discrete equation. Moreover, it only needs to arrange discrete points in the global domain. Thus, the complicated process of mesh formation is avoided and influences from mesh distortion are reduced. A new method for treating crack growth by particle methods has been proposed by Rabczuk and Belytschko [23]. The crack is treated as a collection of cracked particles. At each cracked particle, a discontinuity along a line in 2D or a plane in 3D is introduced, where the normal depends on the complete constitutive model of the material. Shi et al. [24] presented an extended meshless method based on the partition of units used for concurrent multiple crack situations and multiple crack simulations. This method describes the discontinuous displacement field and crack tip singularity field caused by embedding discontinuous items and the crack tip singularity field function into a conventional meshless approximation function. Nanthakumar et al. [25] developed an algorithm to detect and quantify defects in piezoelectric plates. The inverse problem is solved iteratively where XFEM is used for solving the forward problem in each iteration. Béchet et al. [26] applied XFEM to the fracture of piezoelectric materials. Nguyen-Vinh et al. [27] present an extended finite element formulation for dynamic fracture of piezoelectric materials.VCCT was put forward in 1977 by Rybicki and Kanninen [28]. Xie and Biggers [29, 30] had done a lot of research work for VCCT. Compared with the extrapolation method and local or entire equivalent domain integrals, VCCT has an obvious advantage in solving fracture parameters [31–33]. It only uses the nodal force and displacement to calculate the strain energy release rate and only requires a single step in the numerical analysis, thereby simplifying the problem and giving the additional advantages such as high precision and efficiency, no need for special processing of the crack tip unit and small grid size requirements [29, 34, 35]. To date, there are no reports on the virtual crack closure of the electromechanical coupling field.In this paper, a piezoelectric element tailored for VCCT was used to study the crack of piezoelectric materials. CSFEM and VCCT were introduced into fracture mechanics of piezoelectric materials and CSFEM-VCCT for piezoelectric material with cracks was put forward. The energy release rates of different piezoelectric materials with cracks are discussed and compared with FEM-VCCT.


## 2. Governing Equations

The constitutive equations for a two-dimensional piezoelectric material in thex-z axis can be expressed in terms of the strains and the electric field:(1)σp=cpqεq-ekpEk,Di=eiqεq+ξikEk,where σ, ε, D, and E are the stress tensor, the strain tensor, the electric displacement vector, and the electric field vector, respectively. c, e, and ξ are the elastic stiffness, piezoelectric, and dielectric constants, respectively.The strain matrix is related to displacements by(2)εij=(ui,j+uj,i)2.The strain displacement relation can be expressed using the condensed matrix notation given in [13](3)εx=εxx=u,x,εz=εzz=w,z,γxz=2εxz=u,z+w,x,where u and w are the displacement in the x- and z-directions, respectively. Commas followed by indices represent differentiation with respect to that index (i.e., u,x=∂u/∂x).The electric field is related to electric potential by(4)Ei=-ϕ,i.The mechanical equilibrium is governed by(5)σij,j=0.And the governing electrostatic equilibrium is given by (6)Di,i=0.The two-dimensional matrix form of the mechanical and electrical constitutive equations is given by [15](7)εxεzγxz=s11s130s31s33000s55σxσzτxz-0d310d33d150DxDz,ExEz=00d15d31d330σxσzτxz+ξ11σ00ξ33σDxDz,where sij are the elastic compliance constants, dij are piezoelectric constants, and ξiiσ are the dielectric constants. The superscript σ represents quantities measured at constant stress.The finite element solution for 2D piezoelectric problems using the standard linear element can be expressed as(8)u=∑i=1npNuiqi=Nuq,ϕ=∑i=1npNϕiφi=Nϕφ,where np is the number of nodes of an element; Nu, Nϕ are shape function matrices; and q and φ are the nodal displacement and nodal electric potential vectors, respectively.The corresponding approximations of the linear strainε and electric field E are(9)ε=∇us=Buq,E=-∇ϕ=-Bϕφ,where(10)Bui=Ni,x00Ni,zNi,zNi,x,Bϕi=Ni,xNi,z.Using Hamilton’s principle, the piezoelectric static equations of an element can be obtained as follows:(11)KuueKuϕeKuϕeTKϕϕeqφ=FQin which (12)Kuue=∫ΩBuTCEBudΩ,(13)Kuϕe=∫ΩBuTeTBϕdΩ,(14)Kϕϕe=-∫ΩBϕTξTBϕdΩ,(15)F=∫ΩNTfdΩ+∫ΓNTt-dΓ,(16)Q=∫ΓqNϕTq-dΓ.


## 3. Cell-Based Smoothed Finite Element Method

In the stabilized conforming nodal integration technique, the strainε- and the electric field E- used to evaluate the stiffness matrix are computed by a weighted average of the standard strain and electric field of the finite element method. In particular, at an arbitrary point xk in a smoothing element domain Ωk, they are approximated as follows:(17)ε-xk=∫ΩkεxΦkx-xkdΩ,E-xk=∫ΩkExΦkx-xkdΩ,where Ωk is a constant smoothing function described by(18)Φkx-xk=1Ak,x∈Ωk,0,x∉Ωk,where Ak=∫ΩkdΩ is the area of the smoothing cell Ωk. The cell-based element approach is illustrated in detail in Figure 1.Figure 1
A schematic of the smoothing subcells and the values of shape functions at nodes.SubstitutingΩk into (17) and applying the divergence theorem, we obtain a smoothed strain and electric field in the domain Ωk(19)ε-xk=1Ak∫ΓknukudΓ,E-xk=1Ak∫ΓknϕkϕdΓ,where nuk and nϕk are matrices associated with units outward normal to the boundary Ωk,(20)nuk=nxk00nzknzknxk,nϕk=nxknzk.By introducing the finite element approximation ofu and ϕ, (19) can be transformed into matrix form as follows:(21)ε-xk=∑i=1ncB-uixkqi,E-xk=∑i=1ncB-ϕixkϕi,where nc is the number of subcells (cell-based element approach),(22)B-uixk=1Ak∫ΓkNinxk00NinzkNinzkNinxkdΓ,B-ϕixk=1Ak∫ΓkNinxkNinzkdΓ.When bilinear quadrilateral elements are used for modeling, a linear completed displacement field along the boundaryΓk is guaranteed. Therefore, one Gaussian point is sufficient for accurate boundary integration along each line segment Γik of the contour Γk of the domain Ωk. This allows (22) to be evaluated as(23)B-uixk=1Ak∑b=1nbNixbGnxk00NixbGnzkNixbGnzkNixbGnxklbk,B-ϕixk=1Ak∑b=1nbNixbGnxkNixbGnzklbk,where nb is the total number of the line segments of the contour Γk and xbG are the midpoint (Gauss point) of each line segment Γbk, whose length and outward unit normal are lbk and nk, respectively. Finally, the element stiffness matrices in (12)–(14) can be rewritten as follows:(24)Kuuk=∑i=1ncB-uiTCEB-uiAk,Kuϕk=∑i=1ncB-uiTeTB-uiAk,Kϕϕk=-∑k=1ncB-ϕiTξTB-ϕiAk.


## 4. Electromechanical Virtual Crack Closure Technique

The electromechanical VCCT is put forward based on incremental expansion which required equivalent work of the crack closure in potential energy. The VCCT conducts lateral expansions, based on the assumption of the included potential and displacement functions. These are processed using the virtual crack extension and correspond to the piezoelectric element, where the potential was considered a component of “displacement.”Figure2 shows the definition and node numbering of a typical VCCT piezoelectric element for 2D fracture problems. Each element has five nodes. When such an element is applied, it is placed in such a way that nodes 1 and 2 are located at the crack tip, with nodes 3 and 4 behind and node 5 ahead of the crack tip. The element contains two sets of node groups: the top set (nodes 1, 3, and 5) and the bottom set (nodes 2 and 4).Figure 2
A schematic of the fracture of piezoelectric element with dummy nodes.A high stiffness spring is placed between nodes 1 and 2 to compute the nodal forces at the crack tip by(25)Fx=Kxu1-u2,Fy=Kyv1-v2,Q=Kϕϕ1-ϕ2,where (u1, v1) and ϕ1 are the displacement components and electric potential for node 1 referring to the global coordinate system (X, Y), while (u2, v2) and ϕ2 are those for node 2. Kx, Ky, and Kϕ are the spring stiffness corresponding to the X, Y, and ϕ, respectively. Initially, these parameters are set to be large numbers [36]; then once the crack is predicted to grow, they are set to zero.Dummy nodes 3, 4, and 5 do not have contributions to the stiffness matrix and they are introduced only to extract information for displacement opening behind the crack tip and the crack jump length ahead of the crack tip. Since nodes 3 and 4 are behind the crack tip, the displacement openings are(26)Δu=u3-u4,Δv=v3-v4,where (u3, v3) and (u4, v4) are the displacement components for node 3 and node 4, respectively, referring to the global coordinate system (X, Y). The crack jump length is the distance between nodes 1 and 5 and is given by(27)Δa=x5-x12+y5-y12,where (x1,y1) and (x5,y5) are the global coordinates for node 1 and node 5, respectively.When updated at each step, the crack orientation is also updated. This is of particular interest when large deformations cannot be neglected.In order to separate the fracture modes I and II, the strain energy release rates (GI and GII) are computed with respect to the local coordinate system (X-,Y-) attached to the crack tip as shown in Figure 2. The included angle between X and Y is determined by(28)cos⁡θ=x5-x1Δa,sinθ=y5-y1Δa.The nodal forces and displacement are projected into the local coordinate system (X-,Y-) as(29)F-x=Fxcos⁡θ+Fysinθ,F-y=-Fxsinθ+Fycos⁡θ,Q-=Q,Δu~=Δucos⁡θ+Δvsinθ,Δv~=-Δusinθ+Δvcos⁡θ,ϕ~=Δϕ.Based on 2D-VCCT, the energy release rates are approximated as the product of the nodal forces at the crack tip and the nodal displacement openings behind the crack tip by the relations:(30)GI=F~yΔv~2BΔa,GII=F~xΔu~2BΔa,GD=Q~Δϕ~2BΔa,where B is the thickness of the body.


## 5. Numerical Results

### 5.1. Convergence Study in the Energy Norm

The Griffith-Irwin crack in an infinite plate is the first example of a convergence study in this field. As shown in Figure3, the example deals with a central crack of 2a=2 m in a rectangular plate with dimensions h=8 m and w=8 m. The crack is subject to mechanical and electrical loads represented by stress, σ∞=10 MPa, and electric displacement, D∞=10-3 C/m2. The calculations are performed using the piezoelectric material, PZT-4 with a poling direction perpendicular to the crack faces. The material constants of PZT-4 are given in Table 1.Table 1
Material constants.
Material
Elastic constants/(1010 N/m2)
Piezoelectric constant/(C/m2)
Dielectric constant/(1010 C/Vm)
c
11
c
12
c
13
c
44
c
55
e
31
e
33
e
15
d
11
d
33
PZT-4
13.9
7.78
7.43
2.56
11.3
-
6.98
13.84
13.44
60.0
54.7
P-7
13.0
8.3
8.3
2.5
11.9
-
10.3
14.7
13.5
171.0
186.0
PZT-5H
12.6
7.95
8.41
2.3
11.7
-
6.5
23.3
17.44
150.3
130.0Figure 3
Griffith crack under electromechanical loading.The total energy of the system is given by [26](31)W=12∫cijklεijεkl+κijEiEjdΩand the error in the energy norm is then given by (32)err⁡w=12∫cijklεij-εijexεkl-εklex+κijEi-EiexEj-EjexdΩ121/2,where εijex is the exact displacement and Eiex is the exact electric field.The computations were made with an FEM model, whose mesh is structured and has been gradually refined. Figure4 shows the mesh used for the convergence study in the case of a geometrical enrichment, for s=1/12. Figure 5 shows the relationship between the mesh density and the error in the energy norm on a log-log scale. The comparison of the errors among these three different methods shows that the order of numerical accuracy from best to worst is CSFEM-VCCT, XFEM, and FEM-VCCT. The results indicate that the use of CSFEM-VCCT for solving the fracture problem in piezoelectric structures is correct and effective.Figure 4
The mesh used in the convergence study, for the case of a geometrical enrichment, fors=1/12.Figure 5
Convergence for three methods, with coupling loadingσ∞=10 MPa and D∞=10-3 C/m2.


### 5.2. A Central Crack in a Rectangular Plate

Figure6 shows the polarization direction as P. The distance of the central crack along the x-direction is 2a, side length of the plate is 40 cm, σ∞ is the uniform stress σ∞=1×105Pa, and D∞ is the uniform electric displacement D∞=7.5×10-5C/m2. The piezoelectric materials PZT-4, P7, and PZT-5H were adopted for numerical simulation. The material parameters are shown in Table 1. The theoretical solution of the energy release rate for this problem is as follows:PZT-4:(33)G=0.3629σ∞2+0.373σ∞E∞-138.3E∞2×10-10a,P7:(34)G=0.4248σ∞2-0.6952σ∞E∞-398.44E∞2×10-10a,PZT-5H:(35)G=0.4068σ∞2-0.446σ∞E∞-428.5E∞2×10-10a.Figure 6
The piezoelectric model with a central crack.As a consequence of the symmetry of the problem, only a quarter of the plate needs to be modeled. In order to verify the reliability of CSFEM-VCCT for the crack length2a=2 cm, a discrete model of three grids (I: 30 × 30 elements, II: 60 × 60 elements, III: 90 × 90 elements) is simulated to the point of fracture. Every 4-node calculation grid adopts four smooth elements and is compared with the results of the FEM and the theoretical solution.From Table2, we can observe that the result of the release rate of the three materials in each of the three models has a high precision. The precision of the result calculated by CSFEM-VCCT is higher than that of FEM-VCCT. Also model III of PZT-4 was used under the CSFEM-VCCT and FEM-VCCT using four-gauss integral calculation efficiency. The simulation using SFEM-VCCT took 30.664 seconds, while using FEM-VCCT the simulation took 31.725 seconds. The CPU used was an Intel(R) Core(TM) i5-3470 3.20 GHz, RAM: 8 G. The efficiency of CSFEM-VCCT simulation showed some improvement but was not obvious.Table 2
The energy release rate under different mesh dividing methods.
Element
Material
G
I/(×10−3 N/m)
G
D/(×10−3 N/m)
CSFEM-VCCT
Error/(%)
FEM-VCCT
Error/(%)
Analytical solu.
CSFEM-VCCT
Error/(%)
FEM-VCCT
Error/(%)
Analytical solu.
I
PZT-4
3.6614
1.63
3.6603
1.66
3.7221
-
3.2960
1.61
-
3.4072
1.71
-
3.35
P-7
3.9268
2.62
3.9171
2.86
4.0325
-
1.0992
1.75
-
1.1392
1.83
-
1.1188
PZT-5H
4.1140
2.27
4.1052
2.48
4.2096
-
0.4253
2.55
-
0.4480
2.65
-
0.4365
II
PZT-4
3.6685
1.44
3.6632
1.58
3.7221
-
3.2990
1.52
-
3.2964
1.6
-
3.35
P-7
3.9353
2.41
3.9292
2.56
4.0325
-
1.1011
1.58
-
1.1003
1.65
-
1.1188
PZT-5H
4.1174
2.19
4.106
2.46
4.2096
-
0.4260
2.39
-
0.4255
2.51
-
0.4365
III
PZT-4
3.6781
1.18
3.6702
1.39
3.7221
-
3.3071
1.28
-
3.3031
1.4
0
-
3.35
P-7
3.9425
2.23
3.9381
2.34
4.0325
-
1.1043
1.29
-
1.1031
1.39
-
1.1188
PZT-5H
4.1275
1.95
4.1139
2.27
4.2096
-
0.4303
1.41
-
0.4299
1.51
-
0.4365As shown in Table3, for the different length of crack, every four-node grids adopting four smoothed elements, the results show that the CSFEM-VCCT has a high precision.Table 3
Energy release rate under different crack lengths.
Material
Methods
Energy release rate/(×10−3 N/m)
a
0.5 cm
1.0 cm
1.5 cm
2.0 cm
2.5 cm
PZT-4
CSFEM-VCCT
G
I
1.8474
3.6781
5.507
7.4288
9.1384
G
D
-
1.6534
-
3.3071
-
4.9261
-
6.6558
-
8.1681
FEM-VCCT
G
I
1.8353
3.6702
5.5023
7.3404
9.1275
G
D
-
1.6512
-
3.3031
-
4.9054
-
6.6062
-
8.1507
Analytical solu.
G
I
1.861
3.7221
5.5833
7.4440
9.3049
G
D
-
1.6724
-
3.3500
-
5.0174
-
6.6899
-
8.3624
P-7
CSFEM-VCCT
G
I
1.9817
3.9425
5.9412
7.9503
9.8892
G
D
-
0.5499
-
1.1043
-
1.6569
-
2.216
-
2.7444
FEM-VCCT
G
I
1.9689
3.9381
5.9137
7.8850
9.8562
G
D
-
0.5455
-
1.1031
-
1.6524
-
2.2086
-
2.7361
Analytical solu.
G
I
2.0163
4.0325
6.0188
8.0650
10.0813
G
D
-
0.5593
-
1.1188
-
1.6782
-
2.2376
-
2.7969
PZT-5H
CSFEM-VCCT
G
I
2.0536
4.1275
6.1715
8.2270
10.3746
G
D
-
0.2107
-
0.4303
-
0.6461
-
0.8463
-
1.0581
FEM-VCCT
G
I
2.0504
4.1139
6.1676
8.2178
10.2847
G
D
-
0.2082
-
0.4299
-
0.6401
-
0.8359
-
1.0274
Analytical solu.
G
I
2.1048
4.2096
6.3144
8.4192
10.524
0
G
D
-
0.2183
-
0.4365
-
0.6514
-
0.8730
-
1.0912The length of the crack for PZT-4 was2a=2 cm. The result of the normalized energy release rate when the number of the smoothed elements is 1, 2, 4, 8, and 16 is depicted in Table 4. The precision of the results calculated by CSFEM-VCCT is higher than that of the FEM-VCCT when the number of the smoothed elements is changed to 2. As the number of the smoothed elements increase, the accuracy of the whole problem is improved. When the number of the smoothed elements is 4, the precision of the results calculated by CSFEM-VCCT is very high. Therefore, we have shown that the validity and reliability of the simulation using CSFEM-VCCT method are proved.Table 4
Normalized energy release rate under different smoothing subcells.
Normalized energy release rate
FEM-VCCT
CSFEM-VCCT
1
2
4
8
16
PZT-4
G
I
/
G
I exact
0.9861
0.9688
0.9799
0.9882
1.0095
1.011
G
D
/
G
D
exact
0.9860
0.9741
0.9789
0.9872
0.9912
0.9916
P-7
G
I
/
G
I exact
0.9766
0.9630
0.9699
0.9777
0.9777
0.9802
G
D
/
G
D
exact
0.9861
0.9744
0.9784
0.9871
0.9873
0.9901
PZT-5H
G
I
/
G
I exact
0.9773
0.9646
0.9715
0.9805
0.9803
0.9841
G
D
/
G
D
exact
0.9849
0.9758
0.9815
0.9859
0.9891
1.0142


### 5.3. Central Inclined Crack in a Rectangular Plate

As the final example, a central crack inclined 45° to the horizontal direction in a rectangular PZT-4 plate is considered (see Figure7). Uniform tension σ∞=1 MPa and electric displacement D∞=1 C/m2 are applied in the y-direction and the ratios of crack length to width a/w=0.1 (in dimension, h=1 m, a=0.1 m).Figure 7
A central and inclined crack in a rectangular plate.The mechanical and electrical intensity factors of the piezoelectric material PZT-4 for an inclined crack in a rectangle plate are listed in Table5. From the comparison, it is obvious that the CSFEM-VCCT can produce more accurate results than FEM-VCCT using the same meshes (3248 elements). These results show that the smoothing technique adopted in this work improves the calculation of fracture parameters in piezoelectric materials.Table 5
Error comparison of two methods of the mechanical and electrical intensity factors for an inclined crack in a rectangular plate under coupledσ=1 MPa and D=1 C/m2.
Methods
Analytical solu.
CSFEM-VCCT
FEM-VCCT
K
I/MNm−3/2
0.28025
0.28713
0.29317
Error (%)
—
2.455
4.610
K
II/MNm−3/2
0.28025
0.28682
0.29126
Error (%)
—
2.344
3.928
K
D/Cm−3/2
0.39633
0.40036
0.41234
Error (%)
—
1.016
3.030


## 5.1. Convergence Study in the Energy Norm

The Griffith-Irwin crack in an infinite plate is the first example of a convergence study in this field. As shown in Figure3, the example deals with a central crack of 2a=2 m in a rectangular plate with dimensions h=8 m and w=8 m. The crack is subject to mechanical and electrical loads represented by stress, σ∞=10 MPa, and electric displacement, D∞=10-3 C/m2. The calculations are performed using the piezoelectric material, PZT-4 with a poling direction perpendicular to the crack faces. The material constants of PZT-4 are given in Table 1.Table 1
Material constants.
Material
Elastic constants/(1010 N/m2)
Piezoelectric constant/(C/m2)
Dielectric constant/(1010 C/Vm)
c
11
c
12
c
13
c
44
c
55
e
31
e
33
e
15
d
11
d
33
PZT-4
13.9
7.78
7.43
2.56
11.3
-
6.98
13.84
13.44
60.0
54.7
P-7
13.0
8.3
8.3
2.5
11.9
-
10.3
14.7
13.5
171.0
186.0
PZT-5H
12.6
7.95
8.41
2.3
11.7
-
6.5
23.3
17.44
150.3
130.0Figure 3
Griffith crack under electromechanical loading.The total energy of the system is given by [26](31)W=12∫cijklεijεkl+κijEiEjdΩand the error in the energy norm is then given by (32)err⁡w=12∫cijklεij-εijexεkl-εklex+κijEi-EiexEj-EjexdΩ121/2,where εijex is the exact displacement and Eiex is the exact electric field.The computations were made with an FEM model, whose mesh is structured and has been gradually refined. Figure4 shows the mesh used for the convergence study in the case of a geometrical enrichment, for s=1/12. Figure 5 shows the relationship between the mesh density and the error in the energy norm on a log-log scale. The comparison of the errors among these three different methods shows that the order of numerical accuracy from best to worst is CSFEM-VCCT, XFEM, and FEM-VCCT. The results indicate that the use of CSFEM-VCCT for solving the fracture problem in piezoelectric structures is correct and effective.Figure 4
The mesh used in the convergence study, for the case of a geometrical enrichment, fors=1/12.Figure 5
Convergence for three methods, with coupling loadingσ∞=10 MPa and D∞=10-3 C/m2.


## 5.2. A Central Crack in a Rectangular Plate

Figure6 shows the polarization direction as P. The distance of the central crack along the x-direction is 2a, side length of the plate is 40 cm, σ∞ is the uniform stress σ∞=1×105Pa, and D∞ is the uniform electric displacement D∞=7.5×10-5C/m2. The piezoelectric materials PZT-4, P7, and PZT-5H were adopted for numerical simulation. The material parameters are shown in Table 1. The theoretical solution of the energy release rate for this problem is as follows:PZT-4:(33)G=0.3629σ∞2+0.373σ∞E∞-138.3E∞2×10-10a,P7:(34)G=0.4248σ∞2-0.6952σ∞E∞-398.44E∞2×10-10a,PZT-5H:(35)G=0.4068σ∞2-0.446σ∞E∞-428.5E∞2×10-10a.Figure 6
The piezoelectric model with a central crack.As a consequence of the symmetry of the problem, only a quarter of the plate needs to be modeled. In order to verify the reliability of CSFEM-VCCT for the crack length2a=2 cm, a discrete model of three grids (I: 30 × 30 elements, II: 60 × 60 elements, III: 90 × 90 elements) is simulated to the point of fracture. Every 4-node calculation grid adopts four smooth elements and is compared with the results of the FEM and the theoretical solution.From Table2, we can observe that the result of the release rate of the three materials in each of the three models has a high precision. The precision of the result calculated by CSFEM-VCCT is higher than that of FEM-VCCT. Also model III of PZT-4 was used under the CSFEM-VCCT and FEM-VCCT using four-gauss integral calculation efficiency. The simulation using SFEM-VCCT took 30.664 seconds, while using FEM-VCCT the simulation took 31.725 seconds. The CPU used was an Intel(R) Core(TM) i5-3470 3.20 GHz, RAM: 8 G. The efficiency of CSFEM-VCCT simulation showed some improvement but was not obvious.Table 2
The energy release rate under different mesh dividing methods.
Element
Material
G
I/(×10−3 N/m)
G
D/(×10−3 N/m)
CSFEM-VCCT
Error/(%)
FEM-VCCT
Error/(%)
Analytical solu.
CSFEM-VCCT
Error/(%)
FEM-VCCT
Error/(%)
Analytical solu.
I
PZT-4
3.6614
1.63
3.6603
1.66
3.7221
-
3.2960
1.61
-
3.4072
1.71
-
3.35
P-7
3.9268
2.62
3.9171
2.86
4.0325
-
1.0992
1.75
-
1.1392
1.83
-
1.1188
PZT-5H
4.1140
2.27
4.1052
2.48
4.2096
-
0.4253
2.55
-
0.4480
2.65
-
0.4365
II
PZT-4
3.6685
1.44
3.6632
1.58
3.7221
-
3.2990
1.52
-
3.2964
1.6
-
3.35
P-7
3.9353
2.41
3.9292
2.56
4.0325
-
1.1011
1.58
-
1.1003
1.65
-
1.1188
PZT-5H
4.1174
2.19
4.106
2.46
4.2096
-
0.4260
2.39
-
0.4255
2.51
-
0.4365
III
PZT-4
3.6781
1.18
3.6702
1.39
3.7221
-
3.3071
1.28
-
3.3031
1.4
0
-
3.35
P-7
3.9425
2.23
3.9381
2.34
4.0325
-
1.1043
1.29
-
1.1031
1.39
-
1.1188
PZT-5H
4.1275
1.95
4.1139
2.27
4.2096
-
0.4303
1.41
-
0.4299
1.51
-
0.4365As shown in Table3, for the different length of crack, every four-node grids adopting four smoothed elements, the results show that the CSFEM-VCCT has a high precision.Table 3
Energy release rate under different crack lengths.
Material
Methods
Energy release rate/(×10−3 N/m)
a
0.5 cm
1.0 cm
1.5 cm
2.0 cm
2.5 cm
PZT-4
CSFEM-VCCT
G
I
1.8474
3.6781
5.507
7.4288
9.1384
G
D
-
1.6534
-
3.3071
-
4.9261
-
6.6558
-
8.1681
FEM-VCCT
G
I
1.8353
3.6702
5.5023
7.3404
9.1275
G
D
-
1.6512
-
3.3031
-
4.9054
-
6.6062
-
8.1507
Analytical solu.
G
I
1.861
3.7221
5.5833
7.4440
9.3049
G
D
-
1.6724
-
3.3500
-
5.0174
-
6.6899
-
8.3624
P-7
CSFEM-VCCT
G
I
1.9817
3.9425
5.9412
7.9503
9.8892
G
D
-
0.5499
-
1.1043
-
1.6569
-
2.216
-
2.7444
FEM-VCCT
G
I
1.9689
3.9381
5.9137
7.8850
9.8562
G
D
-
0.5455
-
1.1031
-
1.6524
-
2.2086
-
2.7361
Analytical solu.
G
I
2.0163
4.0325
6.0188
8.0650
10.0813
G
D
-
0.5593
-
1.1188
-
1.6782
-
2.2376
-
2.7969
PZT-5H
CSFEM-VCCT
G
I
2.0536
4.1275
6.1715
8.2270
10.3746
G
D
-
0.2107
-
0.4303
-
0.6461
-
0.8463
-
1.0581
FEM-VCCT
G
I
2.0504
4.1139
6.1676
8.2178
10.2847
G
D
-
0.2082
-
0.4299
-
0.6401
-
0.8359
-
1.0274
Analytical solu.
G
I
2.1048
4.2096
6.3144
8.4192
10.524
0
G
D
-
0.2183
-
0.4365
-
0.6514
-
0.8730
-
1.0912The length of the crack for PZT-4 was2a=2 cm. The result of the normalized energy release rate when the number of the smoothed elements is 1, 2, 4, 8, and 16 is depicted in Table 4. The precision of the results calculated by CSFEM-VCCT is higher than that of the FEM-VCCT when the number of the smoothed elements is changed to 2. As the number of the smoothed elements increase, the accuracy of the whole problem is improved. When the number of the smoothed elements is 4, the precision of the results calculated by CSFEM-VCCT is very high. Therefore, we have shown that the validity and reliability of the simulation using CSFEM-VCCT method are proved.Table 4
Normalized energy release rate under different smoothing subcells.
Normalized energy release rate
FEM-VCCT
CSFEM-VCCT
1
2
4
8
16
PZT-4
G
I
/
G
I exact
0.9861
0.9688
0.9799
0.9882
1.0095
1.011
G
D
/
G
D
exact
0.9860
0.9741
0.9789
0.9872
0.9912
0.9916
P-7
G
I
/
G
I exact
0.9766
0.9630
0.9699
0.9777
0.9777
0.9802
G
D
/
G
D
exact
0.9861
0.9744
0.9784
0.9871
0.9873
0.9901
PZT-5H
G
I
/
G
I exact
0.9773
0.9646
0.9715
0.9805
0.9803
0.9841
G
D
/
G
D
exact
0.9849
0.9758
0.9815
0.9859
0.9891
1.0142


## 5.3. Central Inclined Crack in a Rectangular Plate

As the final example, a central crack inclined 45° to the horizontal direction in a rectangular PZT-4 plate is considered (see Figure7). Uniform tension σ∞=1 MPa and electric displacement D∞=1 C/m2 are applied in the y-direction and the ratios of crack length to width a/w=0.1 (in dimension, h=1 m, a=0.1 m).Figure 7
A central and inclined crack in a rectangular plate.The mechanical and electrical intensity factors of the piezoelectric material PZT-4 for an inclined crack in a rectangle plate are listed in Table5. From the comparison, it is obvious that the CSFEM-VCCT can produce more accurate results than FEM-VCCT using the same meshes (3248 elements). These results show that the smoothing technique adopted in this work improves the calculation of fracture parameters in piezoelectric materials.Table 5
Error comparison of two methods of the mechanical and electrical intensity factors for an inclined crack in a rectangular plate under coupledσ=1 MPa and D=1 C/m2.
Methods
Analytical solu.
CSFEM-VCCT
FEM-VCCT
K
I/MNm−3/2
0.28025
0.28713
0.29317
Error (%)
—
2.455
4.610
K
II/MNm−3/2
0.28025
0.28682
0.29126
Error (%)
—
2.344
3.928
K
D/Cm−3/2
0.39633
0.40036
0.41234
Error (%)
—
1.016
3.030


## 6. Conclusions

In this paper, a piezoelectric element tailored for VCCT was used to study the fracture parameters in a piezoelectric material. CSFEM and VCCT were introduced into fracture mechanics of piezoelectric materials and CSFEM-VCCT for piezoelectric material with cracks was put forward. The energy release rates of different piezoelectric materials with cracks are discussed and compared with FEM-VCCT. Numerical examples show that the CSFEM-VCCT shows improvement over the FEM-VCCT, whose stiffness is generally too stiff. Moreover CSFEM-VCCT has a higher precision and increased efficiency. The simple method of CSFEM-VCCT can get the mechanical energy release rates (GI,GII) and electrical energy release rates (Ge) by one step. Due to its simplicity, the VCCT piezoelectric element could be a potential tool for engineers to practice piezoelectric fracture analysis. Therefore, CSFEM-VCCT is an efficient numerical method for simulation fracture problems in piezoelectric materials.



---

*Source: 371083-2015-07-14.xml*