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

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 elementmethod (CSFEM) andVCCThave been used to simulate the fracturemechanics 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.


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][7][8][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][32][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.

Governing Equations
The constitutive equations for a two-dimensional piezoelectric material in the - axis can be expressed in terms of the strains and the electric field: where , , , and  are the stress tensor, the strain tensor, the electric displacement vector, and the electric field vector, respectively., , and  are the elastic stiffness, piezoelectric, and dielectric constants, respectively.The strain matrix is related to displacements by The strain displacement relation can be expressed using the condensed matrix notation given in [13] where  and  are the displacement in the and directions, respectively.Commas followed by indices represent differentiation with respect to that index (i.e.,  , = /).
The electric field is related to electric potential by The mechanical equilibrium is governed by And the governing electrostatic equilibrium is given by The two-dimensional matrix form of the mechanical and electrical constitutive equations is given by [15] [ [ where   are the elastic compliance constants,   are piezoelectric constants, and    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 where   is the number of nodes of an element; N  , 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 where Using Hamilton's principle, the piezoelectric static equations of an element can be obtained as follows: in which

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 x  in a smoothing element domain Ω  , they are approximated as follows: where Ω  is a constant smoothing function described by where   = ∫ Ω  Ω is the area of the smoothing cell Ω  .The cell-based element approach is illustrated in detail in Figure 1.
Substituting Ω  into (17) and applying the divergence theorem, we obtain a smoothed strain and electric field in the domain where n   and n   are matrices associated with units outward normal to the boundary Ω  , By introducing the finite element approximation of u and , (19) can be transformed into matrix form as follows: where   is the number of subcells (cell-based element approach), When bilinear quadrilateral elements are used for modeling, a linear completed displacement field along the boundary Γ  is guaranteed.Therefore, one Gaussian point is sufficient

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." Figure 2 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).
A high stiffness spring is placed between nodes 1 and 2 to compute the nodal forces at the crack tip by where ( 1 , V 1 ) and  1 are the displacement components and electric potential for node 1 referring to the global coordinate system (, ), while ( 2 , V 2 ) and  2 are those for node 2.   ,   , and   are the spring stiffness corresponding to the , , 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 where ( 3 , V 3 ) and ( 4 , V 4 ) are the displacement components for node 3 and node 4, respectively, referring to the global coordinate system (, ).The crack jump length is the distance between nodes 1 and 5 and is given by where ( 1 ,  1 ) and ( 5 ,  5 ) 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 ( I and  II ) are computed with respect to the local coordinate system (, ) attached to the crack tip as shown in Figure 2. The included angle between  and  is determined by The nodal forces and displacement are projected into the local coordinate system (, ) as  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: where  is the thickness of the body.

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 Figure 3, the example deals with a central crack of 2 = 2 m in a rectangular plate with dimensions ℎ = 8 m and  = 8 m.The crack is subject to mechanical and electrical loads represented by stress,  ∞ = 10 MPa, and electric displacement,  ∞ = 10 −3 C/m 2 .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.
The total energy of the system is given by [26]  = 1 2 ∫ (      +       ) Ω (31) and the error in the energy norm is then given by where  ex  is the exact displacement and  ex  is the exact electric field.The computations were made with an FEM model, whose mesh is structured and has been gradually refined.Figure 4 shows the mesh used for the convergence study in the case of a geometrical enrichment, for  = 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.

A Central Crack in a Rectangular
Plate. Figure 6 shows the polarization direction as .The distance of the central crack along the -direction is 2, side length of the plate is 40 cm,  ∞ is the uniform stress  ∞ = 1 × 10 5 Pa, and  ∞ is the uniform electric displacement  ∞ = 7.5×10 −5 C/m 2 .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: P7: PZT-5H: 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 length 2 = 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 Table 2, 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.664seconds, while using FEM-VCCT the simulation took 31.725seconds.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.
As shown in Table 3, 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.
The length of the crack for PZT-4 was 2 = 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.

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 Figure 7).Uniform tension  ∞ = 1MPa and electric displacement  ∞ = 1 C/m 2 are applied in the -direction and the ratios of crack length to width / = 0.1 (in dimension, ℎ = 1 m,  = 0.1 m).
The mechanical and electrical intensity factors of the piezoelectric material PZT-4 for an inclined crack in a rectangle plate are listed in Table 5.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.

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 ( I ,  II ) and electrical energy release rates (  ) 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

Figure 1 :
Figure 1: A schematic of the smoothing subcells and the values of shape functions at nodes.

Figure 2 :
Figure 2: A schematic of the fracture of piezoelectric element with dummy nodes.

Figure 4 :Figure 5 :
Figure 4: The mesh used in the convergence study, for the case of a geometrical enrichment, for  = 1/12.

Figure 6 :
Figure 6: The piezoelectric model with a central crack.

Figure 7 :
Figure 7: A central and inclined crack in a rectangular plate.

Table 2 :
The energy release rate under different mesh dividing methods.

Table 3 :
Energy release rate under different crack lengths.

Table 4 :
Normalized energy release rate under different smoothing subcells.

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  = 1 C/m 2 .