A Meshless Radial Basis Function Based on Partition of Unity Method for Piezoelectric Structures

Ameshless radial basis function based on partition of unitymethod (RBF-PUM) is proposed to analyse static problems of piezoelectric structures. The methods of radial basis functions (RBFs) possess some merits: the shape functions have high order continuity; h-adaptivity is simpler than mesh-based methods; the shape functions are easily implemented in high dimensional space. The partition of unity method (PUM) easily constructs local approximation. The character of local approximate space can be varied and regarded as p-adaptivity. Considering the good properties of the two methods, the RBFs are used for local approximation and the local supported weight functions are used in the partition of unity method. The system equations of the RBF-PUM are derived using the variational principle. The field variables are approximated using the RBF-PUM shape functions which inherit all the advantages of the RBF shape functions such as the delta function property. The boundary conditions can be implemented easily. Numerical examples of piezoelectric structures are investigated to illustrate the efficiency of the proposedmethod and the obtained results are compared with analytical solutions and available numerical solutions. The behaviors of some parameters that probably influenced numerical results are also studied in detail.


Introduction
Piezoelectric material is a new advanced material which can be used for developing smart structures with self-controlling and self-monitoring capabilities.The essential feature of piezoelectric material is the capability of energy transformation between electrical energy and mechanical energy.Hence, piezoelectric materials have been widely used as intelligent control devices such as sensors, actuators, ultrasonic transducers, and active damping devices.
Due to the complicated properties of piezoelectric materials, various approaches including analytical and numerical methods have been proposed to deal with the problems involving piezoelectric materials.However, the analytical solutions are rather restricted in general and difficult to find unless some geometry and boundary are relatively simple.Accurate, efficient, and robust numerical methods are required to simulate piezoelectric problems.Some numerical methods have been proposed to analyse piezoelectric problems such as the finite element method (FEM) [1][2][3][4], the boundary element method (BEM) [5][6][7], and the smoothed FEM [8].Among them, the FEM is considered to be an effective method for mechanical and electromechanical coupling problems.Since the FEM is a mesh-based method, the FEM may lead to large errors when mesh distortion appears or low quality meshes are adopted.Mesh distortion may occur in the problems of large deformation, crack propagation, and so forth.An alternative method is called meshless method which can eliminate mesh-based problems.In the meshless method, only a set of scattered nodes is required to represent problem domain and boundary and therefore the adaptivity is simpler than mesh-based methods.
The initial ideal of meshless method dates back to the smooth particle hydrodynamics (SPH) method proposed by Lucy [9] and Gingold and Monaghan [10] for modeling astrophysical phenomena.Then several meshless methods including the diffuse element method (DEM) [11], the element-free Galerkin (EFG) method [12], the reproducing 2 Mathematical Problems in Engineering kernel particle method (RKPM) [13], the ℎ- clouds method [14], the point interpolation method (PIM) [15][16][17], the radial basis function point interpolation method (RPIM) [18,19], the moving Kriging (MK) meshless method [20], and the meshless local Petrov-Galerkin (MLPG) method [21][22][23][24][25] have been proposed to solve partial differential equations and mechanical problems.Most of these methods are carried out based on weak form equations and a set of background cells or local subdomains is required to implement numerical integration.In the MLPG method, the moving least square (MLS) approach [26] is used to construct approximate functions and the major drawback is that the MLS shape functions lack interpolation property.The boundary conditions cannot be imposed directly and some methods should be proposed to overcome this problem, such as the Lagrange multiplier method [12] and the penalty function method [18].In the RPIM, radial basis functions (RBFs) are used to construct the shape functions.The variable is only related to the distance between interpolation point and node and therefore the method is easily implemented in high dimension space.The behaviors of radial basis functions [27][28][29][30][31][32] have been widely investigated.Radial basis functions are also used in the meshless collocation method [33][34][35] which is based on the strong form equations.These methods based on the RBFs have obtained high accurate solutions and have been shown excellent interpolation property.Some methods coupling with the partition of unity method (PUM) [36,37] have also been investigated in mechanical problems in recent years.The main advantages of the PUM are the flexibility in choosing the local approximation which can be used the approaches of the MLS, the PIM, and the RPIM.Rajendran and Zhang [38] developed a new "FE-Meshfree" method using the concepts of partition of unity method.In this method, the FEM shape functions are used for weight functions and the meshless approach is used to construct local approximation.Some of important works including Q4-LSPIM [38,39], Q4-RPIM [40], and T3-RPIM [41] have been presented for static and vibration analysis of mechanical problems and obtained very accurate results.
In this paper, a meshless radial basis function based on partition of unity method (RBF-PUM) is proposed to analyse 2D piezoelectric structures.The RBFs are used for local approximation and Shepard's method [44] is used for the construction of weight functions in the partition of unity method.The RBF-PUM method has been studied by Wendland [45] and Fasshauer [29] and successfully applied in the convection diffusion equations [46] using the collocation method.The field variables are approximated using the RBF-PUM shape functions which inherits all good advantages of the PUM and RBFs.The performance of the RBF-PUM shape functions will be discussed in detail.
The rest of this paper is organized as follows.The formulation of the RBF-PUM method is presented in Section 2. The governing equations and variational equations for 2D piezoelectric problems are given in Section 3. Numerical examples are considered to verify the effectiveness of the proposed method in Section 4. Finally, we end this paper with some conclusions in Section 5.

The Formulation of the RBF-PUM Method
2.1.The Partition of Unity Method.This section recalls briefly the partition of unity method for 2D problem.Let Ω ⊂ R 2 be a problem domain and {Ω  }  =1 be an open cover of Ω satisfying some mild overlapping conditions such that where (x) is an index set of the patches in which the interpolation point x locates and  is the total number of all patches.
For every patch Ω  , we construct a local approximation   (x).According to the concept of the partition of unity method, the global approximation function  ℎ (x) can be expressed as where {  (x)}  =1 are weight functions.The weight functions are compactly supported, nonnegative, and continuous and satisfy the partition of unity: Equation (3) shows that   (x) = 0 for  ∉ (x).Therefore, (2) can be rewritten as (4)

Construction of Weight Functions with Partition of Unity
Property.For the partition of unity method, the patch can be of any shape, such as square, circular, spherical, or hyperspherical.The basic requirement for all patch domains is that they cover problem domain and boundary.In this paper, the circular patch is used for the analysis in the PUM.Let {  }  =1 be a set of central points of patches {Ω  }

𝑖=1
and {  }  =1 be the radiuses as shown in Figure 1.Nonnegative and compactly supported weight functions {  }  =1 can be constructed using Shepard's method [44] as follows: where function   (x) is given by In the present study, Wendland function [47] is used for function (): in which where R(, ) and p(, ) represent radial basis and polynomial functions, respectively.  is the total number of nodes in patch Ω  and  represents the number of polynomial terms.a and b are unknown vectors to be determined.
In the radial basis function   (, ), the variable is only related to the space distance between interpolation point x and node x  and therefore the RBF method is easily implemented in high dimensional problems.Some radial basis functions such as multiquadric radial basis function (MQ-RBF), thin plate spline radial basis function (TPS-RBF), and Gaussian radial basis function (EXP-RBF) have been widely investigated in [18,[27][28][29][30][31][32].For the present method, the MQ-RBF is adopted and given by where   = √ ( −   ) 2 + ( −   ) 2 and  and  are two shape parameters.
Enforcing (9) to pass through all node values in patch Ω  , the following equations are obtained: where ] , Adding polynomial terms is an extra requirement.In order to guarantee unique approximation [48], the following conditions are usually imposed: which can be rewritten in matrix form Vectors a and b can be obtained by (11) and (14), and then substitute into (9), the local approximation function   (, ) is eventually expressed as where The existence of inverse matrix R −1 0 for arbitrary scattered nodes has been given by Wendland [30].Finally, the shape functions associated to all nodes in the patch Ω  can be expressed as

RBF-PUM Shape Functions and Mathematical Properties.
The radial basis functions are selected as local approximation and the global approximation  ℎ (x) in ( 4) can be expressed as where Ũ is nodal displacement vector composed of all the nodes in the union of the patches Ω = ⋃ ∈(x) Ω  and Φ(x) is the vector of the RBF-PUM shape functions given by where n is the total number of nodes in domain Ω.It is noted that the dimension of the RBF shape functions Φ  (x) may be less than n, but they can be expanded to vector with dimension n.
The derivative functions of the RBF-PUM shape functions can be obtained using Leibniz's rule: where ,  ∈ N 2 0 are multi-indices and The RBF-PUM method inherits the advantages of the RBF meshless method and the following are the mathematical properties of the RBF-PUM shape functions: (1) The RBF-PUM shape functions possess the delta function property: (2) The RBF-PUM shape functions possess the property of partition of unity: (3) If the RBFs include at least linear polynomial terms, the RBF-PUM shape functions can ensure an exact reproduction of linear polynomials; that is, (4) Local compactly supported property.
(5) Especially for card((x)) = 1, the RBF-PUM shape functions are the RBF shape functions as well as their derivative functions.
The RBF-PUM shape functions and their derivative functions are shown in Figures 2 and 3 for 1D and 2D space, respectively.In Figure 2, six regularly distributed nodes in interval [0, 1] are used for the construction of the shape functions and plotted in different colours.In Figure 3

Variational Equations for Piezoelectric Structures
The variational equations for linear piezoelectricity are presented in this section.Considering a two-dimensional piezoelectric problem in domain Ω with boundary Γ, the governing equations are given by where , D are stress and electric displacement vector, respectively.f is body force density.The operators ∇  and ∇ in - plane are given by The boundary conditions are given by where u and  denote prescribed displacements and electric potential, respectively.n is the unit outward normal vector.t and  are the surface traction and charge, respectively.The strain vector related to displacement is and the electric field vector is where ] The constitutive relations for piezoelectric structures are given by  = c − e  E D = e + E, (30) where c, e, and  are the elastic stiffness matrix, piezoelectric matrix, and dielectric matrix, respectively.
The constitutive relations can also be expressed explicitly in matrix form: The mechanical and electrical constitutive relations can also be given by where s is the elastic compliance matrix, d is the piezoelectric strain matrix, and  is the dielectric matrix.The relationships between these constant matrixes can be transforming each other.
The generalized energy functional Π is determined by a summation of strain energy, electric energy, and energy of external work as follows: Using the generalized variational principle, the variational form of piezoelectric structures can be derived: Substituting ( 27)-( 30) into (34), the linear piezoelectric equations can be obtained: where

Numerical Examples
Numerical examples for piezoelectric structures are used to study the RBF-PUM.The choice of the shape parameters  and  is an open issue and has been discussed in the RBFbased methods [18,31].When shape parameter  is chosen as an integer, the moment matrix is ill-conditioned.Therefore, shape parameter  cannot be an integer and, generally, the optimal parameters may change with the practical problems.The effects of the shape parameters on numerical results are also discussed in this paper.In some problems, unless otherwise stated, the two shape parameters  = 1.95 and  = 2.0 investigated the performance of the present method.Gauss integration scheme is used for all numerical examples: the boundary and each background cell adopt 4 and 4 × 4 Gauss integration, respectively.
In the RBF-PUM, the patch fill distance is used as a measure of the density of central points and defined as which indicates how well the central points {  }  =1 fill in problem domain.For regularly distributed central points, the patch fill distance is proportional to the distance of central points.
The support radius or the patch size is defined as where   is a nondimensional coefficient.The patch size usually controls the number of nodes in each patch and therefore parameter   can be adjusted according to node density and the patch fill distance.
For the purpose of error estimation and convergence studies, the relative error is defined as where "ex" and "nu" represent exact and numerical solutions, respectively.The strip is subjected a uniform load  0 on the top and bottom boundaries and an applied voltage  0 on the left and right boundaries as shown in Figure 4. Since the applied electric field is perpendicular to polarization direction, it induces a shear strain.The overall deformation of the piezoelectric strip is a superposition of the deformation due to the shear strain and the applied load.The analytical solutions for this shear problem are given by [49] The problem domain is represented by both 11 × 11 regular nodes and 121 irregular nodes in Figure 5.The PU cover consists of  =  2 (5 × 5) regular circular patches and the coefficient   = 1.2 is taken here.10 × 10 square background cells are used for numerical integration.
The obtained results of the present method are plotted in Figure 6, in which symbols R and IR denote regular and irregular nodes, respectively.Figure 6(a) shows the horizontal displacement  along  = 0 and the deflection  along  = 0 is plotted in Figure 6(b).The variation of the electric potential  is depicted in Figure 6(c).The results of the RBF method are also presented in these figures and the size of influence domain of the RBF method is taken as 4.0 times average nodal distance.The results of the RBF-PUM method and the RBF method are very accurate comparing with the analytical solutions.Moreover, the present results are slightly better than those of the RBF method.
The computed mechanical deformation results of the present method for irregular nodes are plotted in Figure 7, which seems quite reasonable for the shear deformation.The convergence curves and convergence rates of displacements  and  and electric potential  using the linear regression are shown in Figure 8 for regular nodes.It can be seen that the present results are very accurate even for coarse node distribution.The relative errors will decrease with the increase of nodes.It is due to the fact that the number of nodes using in the RBF-PUM shape functions increases with the node density and the numerical results become better.
In the RBF-PUM, the common requirement for all patches is that they cover the problem domain and boundary.Numerical results are usually influenced by the patch size, since the number of nodes using in the construction of the RBF-PUM shape function are influenced by the patch size.It is great of interest to study the effects of the patch size on numerical results.The results of displacements and electric potential at point  ( = 3/4,  = 0) are listed in Table 1.The present method gives stable and accurate results for different patch sizes and nodal patterns.The solutions will converge to the analytical solutions for a larger patch size.It is due to the fact that increasing the patch size will lead to an increase of the number of nodes in each patch and the RBF-PUM shape function for a given interpolation point contains more data information for a larger patch size.However, it should be noted that the optimal parameter   may change with the patch fill distance and the coefficient   can be adjusted according to the patch fill distance.

Bending Deformation of Piezoelectric Strip.
The same material and geometry of the piezoelectric strip in bending are studied.The top and bottom surfaces are subjected to an applied voltage and the right side is applied to a linear stress as depicted in Figure 9.The analytical solutions of displacements and electric potential are given by [   and the boundary conditions for this problem are given by To demonstrate the efficiency of the proposed method for irregular central points, 25 circular patches with irregular central points are used for the bending analysis of the piezoelectric strip.The central points and nodes are plotted in Figure 10, in which the patch fill distance is about  = 0.1938 and the coefficient   = 1.2 is also taken here.
The solutions of the RBF-PUM, the RBF method, and the analytical method are shown in Figure 11 an applied voltage , the tip deflection  can be approximated as [51] where  is the length of the beam,  is the total thickness, and  31 is the piezoelectric strain constant.The size of the beam is taken as  = 5 mm and  = 0.4 mm and the tip deflection computed by ( 44) is  = 1.0206 × 10 −8 m.
Considering the thickness-length ratio of the beam, a total of 21 × 4 regular patches and 51 × 5 regular nodes are used to represent the problem domain.Due to its geometric symmetry, only one-half of the patches of the top beam are shown in Figure 13, where the central points are at the upper and down quarter through the thickness direction and the coefficient   = 1.2 is taken.
The tip deflections for different node densities are listed in Table 2 and the present results are very close to the results given by the ABAQUS [19].The effect of the patch size on Mathematical Problems in Engineering x (mm) the tip deflection is investigated in Table 3.As the patch size increases, the number of nodes in each patch will increase as well as the nodes using in the construction of the RBF-PUM shape functions.The obtained solutions will converge for a larger patch size.The shape parameters  and  are also studied in Tables 4 and 5, respectively.It can be observed that increasing the parameter  leads to an increase of the tip deflection, but it changes very slowly and the tip deflection keeps a constant for different values of parameter .All the results computed by present method are very stable and slightly larger than the analytical solution.It may be due to using the equations derived by the Euler-Bernoulli model and the numerical method results in a less stiff model [19].

Piezoelectric Cook's Membrane.
The piezoelectric Cook's membrane problem is solved using the present method.The geometry and applied unity tip load  are shown in Figure 14.The lower surface of the membrane is prescribed by zero voltage.Cook's membrane is made of PZT-4 material provided as follows: 11 : 139 × 10 The analytical solutions of this problem are unavailable and the best known numerical results are given by Long et al. [42] using the finite element method with a fine mesh.The results of the displacement, the electric potential, the first principal stress, and the electric flux density at nodes , ,  are The problem is represented by 7 × 4 circular patches and 25 × 7 regular nodes as shown in Figure 15.The patch fill distance is about  = 5.341 and the coefficient is   = 1.5.The results and relative errors percentage of the present  method for different nodal densities are listed in Table 6 together with the results of the smoothed finite element method [8].In the smoothed FEM, a fine mesh with 24 × 24 elements is used for this problem.It can be seen that the present results of the displacement   and the first principal stress  1 are in very good accord with the results given by the FEM.Although the errors of the electric potential   and the electric flux density   are slightly lager, the present method still achieves reasonable prediction.

A Piezoelectric Bimorph Optical Microscanner.
The purpose of this problem is to simulate the linear title angle of the reflected light through a mirror of an optical microscanner.The device is composed of two parallel bimorphs connected by a mirror at the tip center as depicted in Figure 16.
Considering the geometry of the optical microscanner device, node arrangements and patches are very similar to the  The tip displacements calculated for several applied voltages with different nodal densities are listed in Table 7, in which the first two rows represent background cells and nodal densities, respectively.The numerical results obtained by the meshless point collocation method (PCM) [43] are also given in the table.In the PCM, the reproducing kernel particle approximation is used for the displacement field and electric potential.The governing equations and boundary conditions are satisfied only at specific nodes.In the present method, the system equations are derived by the variational principle and numerical integration is required.Although the present results are very close to the PCM solutions, there are still some differences between these two different numerical methods.
From the tip deflection, the title angle of the mirror can be calculated by (46).The title angles varying with different applied voltages for 51 × 5 regular nodes are shown in Figure 17.The results obtained by the node-based element (NSPE-T3, NSPE-Q4), the cell-based element (SPQ4), and the FEM-T3 [52] are also shown in Figure 17.The linear variation of the angle is in good agreement with the PCM.

Conclusions
In this paper, a meshless radial basis function based on partition of unity method is presented for studying 2D piezoelectric structures.The multiquadric radial basis functions are used for local approximation and Shepard's method is used for the construction of weight functions in the partition of unity method.The weak form equations are derived by the generalized variational principle.The mechanical displacements and electric potential are approximated using the RBF-PUM shape functions which inherit all advantages of the RBF shape functions.The RBF-PUM shape functions have been shown excellent interpolation property and therefore the boundary conditions can be imposed easily without special treatments.
Numerical examples are calculated to illustrate the efficiency of the proposed method for static problems.The effects of some parameters including the shape parameters, the patch size, and nodal density have also been studied in detail.The obtained results are compared with the solutions of the analytical method and other numerical methods such as the RBF method, the finite element method, and the smoothed finite element method.The present results are very accurate and are in excellent agreement with the analytical solutions.
, the shape function and its derivative of point (0, 0) are plotted in interval [−1, 1] × [−1, 1].It can be observed that the RBF-PUM shape functions have interpolation property and local compactly supported property.

Figure 2 :
Figure 2: RBF-PUM shape functions and their derivatives in 1D space.

Figure 4 :
Figure 4: Piezoelectric strip subjected to a uniform load and an applied voltage.

MathematicalFigure 5 :
Figure 5: Node distribution and partition of the piezoelectric strip in shear: (a) regular nodes with regular central points; (b) irregular nodes with regular central points.

Figure 6 :
Figures 11(a) and 11(c), respectively.The variation of the deflection  along the bottom of the piezoelectric strip is depicted in Figure 11(b).It can be observed that the present results match the analytical solutions very well even for irregular central points and the present results are more accurate than the RBF results by comparing with the analytical solutions.4.3.A Parallel Piezoelectric Bimorph Beam.Consider a twolayer parallel piezoelectric bimorph beam as shown in Figure 12.The beam is made of material PVDF with the same thickness and polarization direction and the material properties are summarized as follows: : 2.0 × 10 9 Pa ]: 0.29  31 : 0.046 C/m 2  32 : 0.046 C/m 2  11 : 0.1062 × 10 −9 F/m  33 : 0.1062 × 10 −9 F/m For this problem, a unity voltage (1 V) is applied across the thickness between the surfaces and the intermediate electrode.In this study, a plane stress problem is assumed.For

Figure 7 : 10 (Figure 8 :
Figure 7: The shear deformation of the piezoelectric strip for irregular nodes.Note that the displacements plotted are magnified 100 times.

Figure 9 :
Figure 9: Piezoelectric strip subjected to a linear stress and an applied voltage.

Figure 10 :
Figure 10: Node distribution and partition of the piezoelectric strip in bending: (a) regular nodes with irregular central points; (b) irregular nodes with irregular central points.

Figure 11 :
Figure 11: Comparisons of the present solutions, the RBF solutions, and the analytical solutions of the piezoelectric strip in bending: (a) the variation of the horizontal displacement  along  = ; (b) the variation of the vertical displacement  along  = −ℎ; (c) the variation of the electric potential  along  = .

Figure 13 :
Figure 13: Node distribution and patches of the top beam.

Table 6 :
Numerical results of piezoelectric Cook's membrane and relative errors for different nodal densities.

Table 7 :
Numerical results of piezoelectric optical microscanner device for different nodal densities (m).