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

Authors: Sen Li; Linquan Yao; Shichao Yi; Wei Wang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

A meshless radial basis function based on partition of unity method (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 proposed method 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.


---

## Body

## 1. 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–4], the boundary element method (BEM) [5–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 kernel particle method (RKPM) [13], the h-p clouds method [14], the point interpolation method (PIM) [15–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–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–32] have been widely investigated. Radial basis functions are also used in the meshless collocation method [33–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 Section2. 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.


## 2. 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Ω⊂R2 be a problem domain and {Ωi}i=1M be an open cover of Ω satisfying some mild overlapping conditions such that(1)∀x∈ΩIx=i∣x∈Ωi,1≤cardIx≤M,where I(x) is an index set of the patches in which the interpolation point x locates and M is the total number of all patches.For every patchΩi, we construct a local approximation ui(x). According to the concept of the partition of unity method, the global approximation function uh(x) can be expressed as(2)uhx=∑i=1Mwixuix,where {wi(x)}i=1M are weight functions. The weight functions are compactly supported, nonnegative, and continuous and satisfy the partition of unity:(3)supp⁡wix=Ωii=1,…,M,∑i=1Mwix=1onΩ.Equation (3) shows that wi(x)=0 for i∉I(x). Therefore, (2) can be rewritten as(4)uhx=∑i∈Ixwixuix.


### 2.2. 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{Xi}i=1M be a set of central points of patches {Ωi}i=1M and {ri}i=1M be the radiuses as shown in Figure 1. Nonnegative and compactly supported weight functions {wi}i=1M can be constructed using Shepard’s method [44] as follows:(5)wix=φix∑k=1Mφkx,i=1,…,M,where function φi(x) is given by(6)φix=φx-Xiri.In the present study, Wendland function [47] is used for function φ(r):(7)φr=1-r41+4r0≤r≤1,0r>1,in which(8)r=x-Xiri.Figure 1
Definition of interpolation pointx, central point (red ∗), patch (red circle), and support radius (or patch size). The patches for interpolating point x are defined as Ω1={1,2,4,5,6,8} and Ω2={7,8,9,10,11,12}. The support domain of interpolation point x is defined as Ω^=Ω1∪Ω2={1,2,4,5,6,7,8,9,10,11,12}.


### 2.3. Radial Basis Function Method

In the present method, radial basis functions with polynomial terms are used to construct the local approximation of the partition of unity method. Using radial basis functions with polynomial functions, the local approximationui(x) is given by(9)uix,y=∑i=1niRix,yai+∑j=1mpjx,ybj=Rx,ya+px,yb,where R(x,y) and p(x,y) represent radial basis and polynomial functions, respectively. ni is the total number of nodes in patch Ωi and m represents the number of polynomial terms. a and b are unknown vectors to be determined.In the radial basis functionRi(x,y), the variable is only related to the space distance between interpolation point x and node xi 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–32]. For the present method, the MQ-RBF is adopted and given by(10)Rix,y=1+η2di2q,where di=(x-xi)2+(y-yi)2 and η and q are two shape parameters.Enforcing (9) to pass through all node values in patch Ωi, the following equations are obtained:(11)ui=R0a+Pmb,where(12)aT=a1,a2,…,ani,bT=b1,b2,…,bm,uiT=u1,u2,…,uni,R0=R1x1,y1R2x1,y1⋯Rnix1,y1R1x2,y2R2x2,y2⋯Rnix2,y2⋮⋮⋱⋮R1xni,yniR2xni,yni⋯Rnixni,yni,PmT=11⋯1x1x2⋯xniy1y2⋯yni⋮⋮⋱⋮pmx1pmx2⋯pmxni.Adding polynomial terms is an extra requirement. In order to guarantee unique approximation [48], the following conditions are usually imposed:(13)∑i=1nipjxi,yiai=0,j=1,2,…,m,which can be rewritten in matrix form(14)PmTa=0.Vectorsa and b can be obtained by (11) and (14), and then substitute into (9), the local approximation function ui(x,y) is eventually expressed as(15)uix,y=Φiui,where(16)Φi=Rx,ySa+Px,ySb,Sb=PmTR0Pm-1PmTR0-1,Sa=R0-1I-PmSb.The existence of inverse matrixR0-1 for arbitrary scattered nodes has been given by Wendland [30]. Finally, the shape functions associated to all nodes in the patch Ωi can be expressed as(17)Φix=ϕi1x,ϕi2x,…,ϕinix.


### 2.4. RBF-PUM Shape Functions and Mathematical Properties

The radial basis functions are selected as local approximation and the global approximationuh(x) in (4) can be expressed as(18)uhx=∑i∈IxwixΦixui=Φ~xU~s,where U~s is nodal displacement vector composed of all the nodes in the union of the patches Ω^=⋃i∈I(x)Ωi and Φ~(x) is the vector of the RBF-PUM shape functions given by(19)Φ~x=∑i∈IxwixΦix=ϕ~1x,ϕ~2x,…,ϕ~n^x,where n^ is the total number of nodes in domain Ω^. It is noted that the dimension of the RBF shape functions Φi(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:(20)DαΦ~x=∑i∈Ix∑β≤ααβDβwixDα-βΦix,where α,β∈N02 are multi-indices and Dα=∂|α|/∂xα is a spatial differential operator.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:(21)ϕ~ixj=δij.(2)
The RBF-PUM shape functions possess the property of partition of unity:(22)∑i=1n^ϕ~ix=1.(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,(23)∑i=1n^ϕ~ixxi=x.(4)
Local compactly supported property.(5)
Especially forcard(I(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 Figures2 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, 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
RBF-PUM shape functions and their derivatives in 1D space.Figure 3
RBF-PUM shape function and its derivative of point(0,0) in 2D space.


## 2.1. The Partition of Unity Method

This section recalls briefly the partition of unity method for 2D problem. LetΩ⊂R2 be a problem domain and {Ωi}i=1M be an open cover of Ω satisfying some mild overlapping conditions such that(1)∀x∈ΩIx=i∣x∈Ωi,1≤cardIx≤M,where I(x) is an index set of the patches in which the interpolation point x locates and M is the total number of all patches.For every patchΩi, we construct a local approximation ui(x). According to the concept of the partition of unity method, the global approximation function uh(x) can be expressed as(2)uhx=∑i=1Mwixuix,where {wi(x)}i=1M are weight functions. The weight functions are compactly supported, nonnegative, and continuous and satisfy the partition of unity:(3)supp⁡wix=Ωii=1,…,M,∑i=1Mwix=1onΩ.Equation (3) shows that wi(x)=0 for i∉I(x). Therefore, (2) can be rewritten as(4)uhx=∑i∈Ixwixuix.


## 2.2. 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{Xi}i=1M be a set of central points of patches {Ωi}i=1M and {ri}i=1M be the radiuses as shown in Figure 1. Nonnegative and compactly supported weight functions {wi}i=1M can be constructed using Shepard’s method [44] as follows:(5)wix=φix∑k=1Mφkx,i=1,…,M,where function φi(x) is given by(6)φix=φx-Xiri.In the present study, Wendland function [47] is used for function φ(r):(7)φr=1-r41+4r0≤r≤1,0r>1,in which(8)r=x-Xiri.Figure 1
Definition of interpolation pointx, central point (red ∗), patch (red circle), and support radius (or patch size). The patches for interpolating point x are defined as Ω1={1,2,4,5,6,8} and Ω2={7,8,9,10,11,12}. The support domain of interpolation point x is defined as Ω^=Ω1∪Ω2={1,2,4,5,6,7,8,9,10,11,12}.


## 2.3. Radial Basis Function Method

In the present method, radial basis functions with polynomial terms are used to construct the local approximation of the partition of unity method. Using radial basis functions with polynomial functions, the local approximationui(x) is given by(9)uix,y=∑i=1niRix,yai+∑j=1mpjx,ybj=Rx,ya+px,yb,where R(x,y) and p(x,y) represent radial basis and polynomial functions, respectively. ni is the total number of nodes in patch Ωi and m represents the number of polynomial terms. a and b are unknown vectors to be determined.In the radial basis functionRi(x,y), the variable is only related to the space distance between interpolation point x and node xi 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–32]. For the present method, the MQ-RBF is adopted and given by(10)Rix,y=1+η2di2q,where di=(x-xi)2+(y-yi)2 and η and q are two shape parameters.Enforcing (9) to pass through all node values in patch Ωi, the following equations are obtained:(11)ui=R0a+Pmb,where(12)aT=a1,a2,…,ani,bT=b1,b2,…,bm,uiT=u1,u2,…,uni,R0=R1x1,y1R2x1,y1⋯Rnix1,y1R1x2,y2R2x2,y2⋯Rnix2,y2⋮⋮⋱⋮R1xni,yniR2xni,yni⋯Rnixni,yni,PmT=11⋯1x1x2⋯xniy1y2⋯yni⋮⋮⋱⋮pmx1pmx2⋯pmxni.Adding polynomial terms is an extra requirement. In order to guarantee unique approximation [48], the following conditions are usually imposed:(13)∑i=1nipjxi,yiai=0,j=1,2,…,m,which can be rewritten in matrix form(14)PmTa=0.Vectorsa and b can be obtained by (11) and (14), and then substitute into (9), the local approximation function ui(x,y) is eventually expressed as(15)uix,y=Φiui,where(16)Φi=Rx,ySa+Px,ySb,Sb=PmTR0Pm-1PmTR0-1,Sa=R0-1I-PmSb.The existence of inverse matrixR0-1 for arbitrary scattered nodes has been given by Wendland [30]. Finally, the shape functions associated to all nodes in the patch Ωi can be expressed as(17)Φix=ϕi1x,ϕi2x,…,ϕinix.


## 2.4. RBF-PUM Shape Functions and Mathematical Properties

The radial basis functions are selected as local approximation and the global approximationuh(x) in (4) can be expressed as(18)uhx=∑i∈IxwixΦixui=Φ~xU~s,where U~s is nodal displacement vector composed of all the nodes in the union of the patches Ω^=⋃i∈I(x)Ωi and Φ~(x) is the vector of the RBF-PUM shape functions given by(19)Φ~x=∑i∈IxwixΦix=ϕ~1x,ϕ~2x,…,ϕ~n^x,where n^ is the total number of nodes in domain Ω^. It is noted that the dimension of the RBF shape functions Φi(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:(20)DαΦ~x=∑i∈Ix∑β≤ααβDβwixDα-βΦix,where α,β∈N02 are multi-indices and Dα=∂|α|/∂xα is a spatial differential operator.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:(21)ϕ~ixj=δij.(2)
The RBF-PUM shape functions possess the property of partition of unity:(22)∑i=1n^ϕ~ix=1.(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,(23)∑i=1n^ϕ~ixxi=x.(4)
Local compactly supported property.(5)
Especially forcard(I(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 Figures2 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, 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
RBF-PUM shape functions and their derivatives in 1D space.Figure 3
RBF-PUM shape function and its derivative of point(0,0) in 2D space.


## 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(24)∇sTσ+f=0,∇TD=0,where σ, D are stress and electric displacement vector, respectively. f is body force density. The operators ∇s and ∇ in x-z plane are given by(25)∇s=∂∂x0∂∂z0∂∂z∂∂xT,∇=∂∂x∂∂zT.The boundary conditions are given by(26)σn=t¯,onΓσu=u¯,onΓuDn=-q¯,onΓqϕ=ϕ¯,onΓϕ,where u¯ and ϕ¯ denote prescribed displacements and electric potential, respectively. n is the unit outward normal vector. t¯ and q¯ are the surface traction and charge, respectively.The strain vector related to displacement is(27)ε=∇su=∇sNuu^=Buu^and the electric field vector is(28)E=-∇ϕ=-∇Nϕϕ^=-Bϕϕ^,where(29)u^=u1w1u2w2⋯un^wn^Tϕ^=ϕ1ϕ2⋯ϕn^TNu=ϕ~10⋯ϕ~n^00ϕ~1⋯0ϕ~n^Nϕ=ϕ~1ϕ~2⋯ϕ~n^Bu=ϕ~1,x0⋯ϕ~n^,x00ϕ~1,z⋯0ϕ~n^,zϕ~1,zϕ~1,x⋯ϕ~n^,zϕ~n^,xBϕ=ϕ~1,x⋯ϕ~n^,xϕ~1,z⋯ϕ~n^,z.The constitutive relations for piezoelectric structures are given by(30)σ=cε-eTED=eε+κE,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:(31)σxσzτxz=c11c130c13c33000c55εxεzγxz-0e310e33e150ExEz,DxDz=00e15e31e330εxεzγxz+κ1100κ33ExEz.The mechanical and electrical constitutive relations can also be given by(32)εxεzγxz=s11s130s13s33000s55σxσzτxz-0d310d33d150ExEz,DxDz=00d15d31d330σxσzτxz+ξ1100ξ33ExEz,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:(33)Π=∫Ω12εTσ-12DTE-uTfdΩ-∫ΓuTt¯+ϕTq¯dΓ.Using the generalized variational principle, the variational form of piezoelectric structures can be derived:(34)δΠ=∫ΩδεTσ-δETD-δuTfdΩ-∫ΓδuTt¯+δϕTq¯dΓ=0.Substituting (27)–(30) into (34), the linear piezoelectric equations can be obtained:(35)KuuKuϕKuϕTKϕϕu^ϕ^=FQ,where(36)Kuu=∫ΩBuTcBudΩ,Kuϕ=-∫ΩBuTeTBϕdΩ,Kϕϕ=∫ΩBϕTκBϕdΩ,F=∫ΩNuTfdΩ+∫ΓσNuTt¯dΓ,Q=-∫ΓqNϕTq¯dΓ.


## 4. Numerical Examples

Numerical examples for piezoelectric structures are used to study the RBF-PUM. The choice of the shape parametersq and η is an open issue and has been discussed in the RBF-based methods [18, 31]. When shape parameter q is chosen as an integer, the moment matrix is ill-conditioned. Therefore, shape parameter q 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 q=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(37)H=supx∈Ω⁡min1≤i≤M⁡x-Xi,which indicates how well the central points {Xi}i=1N 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(38)ri=αr×H,where αr is a nondimensional coefficient. The patch size usually controls the number of nodes in each patch and therefore parameter αr 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(39)Er=uex-unu2uex2,where “ex” and “nu” represent exact and numerical solutions, respectively.


### 4.1. Shear Deformation of Piezoelectric Strip

The shear deformation of a square piezoelectric strip polarized inz-direction is taken into consideration. The piezoelectric strip is made of material PZT-5 and its material properties are given as follows:s11:16.4×10-6 (mm)2/Ns13:-7.22×10-6 (mm)2/Ns33:18.8×10-6 (mm)2/Ns55:47.5×10-6 (mm)2/Nσ0:-5.0 N/(mm)2σ1:20.0 N/(mm)2L:1.0 mmd31:-172×10-9 mm/Vd33:374×10-9 mm/Vd15:584×10-9 mm/Vξ11:1.53105×10-8 N/V2ξ33:1.505×10-7 N/V2V0:1000 Vh:0.5 mmThe strip is subjected a uniform loadσ0 on the top and bottom boundaries and an applied voltage V0 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.Figure 4
Piezoelectric strip subjected to a uniform load and an applied voltage.The following boundary conditions of the piezoelectric strip are given by(40)ϕ,zx,±h=0,σzx,±h=σ0,τxzL,z=0,τxzx,±h=0,ϕL,z=-V0,σxL,z=0,ϕ0,z=+V0,u0,z=0,w0,0=0.The analytical solutions for this shear problem are given by [49](41)u=s13σ0,w=d15V0xh+s33σ0z,ϕ=V01-2xL.The problem domain is represented by both11×11 regular nodes and 121 irregular nodes in Figure 5. The PU cover consists of M=m2(5×5) regular circular patches and the coefficient αr=1.2 is taken here. 10×10 square background cells are used for numerical integration.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.
(a) (b)The obtained results of the present method are plotted in Figure6, in which symbols R and IR denote regular and irregular nodes, respectively. Figure 6(a) shows the horizontal displacement u along z=0 and the deflection w along z=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.Figure 6
Comparisons of the present solutions, the RBF solutions, and the analytical solutions of the piezoelectric strip in shear: (a) the variation of the horizontal displacementu along z=0; (b) the variation of the vertical displacement w along z=0; (c) the variation of the electric potential ϕ along z=0.
(a) (b) (c)The computed mechanical deformation results of the present method for irregular nodes are plotted in Figure7, which seems quite reasonable for the shear deformation. The convergence curves and convergence rates of displacements u and w 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.Figure 7
The shear deformation of the piezoelectric strip for irregular nodes. Note that the displacements plotted are magnified 100 times.Figure 8
Convergence property of displacementsu and w and electric potential ϕ for regular nodes.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 pointA (x=3L/4, z=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 αr may change with the patch fill distance and the coefficient αr can be adjusted according to the patch fill distance.Table 1
Numerical results of displacementsuA×10-5 (mm) and wA×10-4 (mm) and electric potential ϕA(V) at point A(x=3L/4,z=0) for different patch sizes.
α r 1.1 1.2 1.3 1.4 1.5 1.6 1.8 2.0 Exact u A R 2.70774 2.70762 2.70747 2.70740 2.70753 2.70752 2.70752 2.70752 2.70750 IR 2.70523 2.70372 2.70676 2.70783 2.70844 2.70735 2.70760 2.70804 w A R 8.76021 8.76018 8.76018 8.76025 8.76045 8.76068 8.76078 8.76089 8.76000 IR 8.76367 8.76041 8.75980 8.75963 8.76029 8.75998 8.75991 8.75998 ϕ A R −499.92 −499.92 −499.93 −499.94 −500.01 −500.00 −500.02 −500.03 −500.00 IR −500.79 −500.05 −499.96 −499.94 −500.03 −499.97 −499.98 −499.98


### 4.2. 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 Figure9. The analytical solutions of displacements and electric potential are given by [50](42)u=s11σ0-d31V0s11hx+s111-d312s11ξ33σ1xz,w=s13σ0-d31V0s11hz+s131-d33d31s13ξ33-s111-d312s11ξ33σ1x22,ϕ=V0zh-d31σ12ξ33h2-z2,and the boundary conditions for this problem are given by(43)ϕ,zx,±h=±V0,σzx,±h=0,τxzx,z±h=0,ϕ,xL,z=0,σxL,z=σ0+σ1z,τxzL,z=0,ϕ,x0,z=0,u0,z=0,w0,0=0.Figure 9
Piezoelectric strip subjected to a linear stress and an applied voltage.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 Figure10, in which the patch fill distance is about H=0.1938 and the coefficient αr=1.2 is also taken here.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.
(a) (b)The solutions of the RBF-PUM, the RBF method, and the analytical method are shown in Figure11. The displacement u and electric potential ϕ along x=L are plotted in Figures 11(a) and 11(c), respectively. The variation of the deflection w 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.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 displacementu along x=L; (b) the variation of the vertical displacement w along z=-h; (c) the variation of the electric potential ϕ along x=L.
(a) (b) (c)


### 4.3. A Parallel Piezoelectric Bimorph Beam

Consider a two-layer parallel piezoelectric bimorph beam as shown in Figure12. The beam is made of material PVDF with the same thickness and polarization direction and the material properties are summarized as follows:E:2.0×109 Paν:0.29e31:0.046 C/m2e32:0.046 C/m2κ11:0.1062×10-9 F/mκ33:0.1062×10-9 F/mFigure 12
A parallel piezoelectric bimorph beam.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 an applied voltageV, the tip deflection δ can be approximated as [51](44)δ=2L2Vd31H2,where L is the length of the beam, H is the total thickness, and d31 is the piezoelectric strain constant. The size of the beam is taken as L=5 mm and H=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 of21×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 αr=1.2 is taken.Figure 13
Node distribution and patches of the top beam.The tip deflections for different node densities are listed in Table2 and the present results are very close to the results given by the ABAQUS [19]. The effect of the patch size on 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 q and η are also studied in Tables 4 and 5, respectively. It can be observed that increasing the parameter q 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].Table 2
The tip deflection of the piezoelectric beam for different nodes whenq=1.95, η=2.0, αr=1.2 (×10-8 m).
Background cell 40 × 4 50 × 4 60 × 4 70 × 4 80 × 4 100 × 4 CPS4E [19] Nodal density 41 × 5 51 × 5 61 × 5 71 × 5 81 × 5 101 × 5  Tip deflection 1.0715 1.0719 1.0705 1.0737 1.0736 1.0719 1.073Table 3
Effects of the patch size on the tip deflection for51×5 nodes when q=1.95, η=2.0 (×10-8 m).
α r 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8  Tip deflection 1.0494 1.0719 1.0721 1.0725 1.0726 1.0738 1.0738 1.0738Table 4
Effects of parameterq on the tip deflection for 51×5 nodes when η=2.0, αr=1.2 (×10-8 m).
q 0.95 0.98 1.03 1.05 1.08 1.92 1.95 1.98  Tip deflection 1.0709 1.0709 1.0710 1.0710 1.0710 1.0719 1.0719 1.0719Table 5
Effects of parameterη on the tip deflection for 51×5 nodes when q=1.95, αr=1.2 (×10-8 m).
η 0.5 0.75 1.0 1.25 1.5 2.0 2.5 3.0  Tip deflection 1.0719 1.0719 1.0719 1.0719 1.0719 1.0719 1.0719 1.0719


### 4.4. Piezoelectric Cook’s Membrane

The piezoelectric Cook’s membrane problem is solved using the present method. The geometry and applied unity tip loadF 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:c11:139×103 N/mm2c13:74.3×103 N/mm2c33:113×103 N/mm2c55:25.6×103 N/mm2e15:13.44×106 pC/mm2e31:-6.98×106 pC/mm2e33:13.84×106 pC/mm2κ11:6.00×109 pC/GVmmκ33:5.47×109 pC/GVmmFigure 14
Cook’s membrane.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 A, B, C are (45)wA=2.109×10-4mm,ϕA=1.732×10-8GVσ1B=0.21613N/mm2,DC=22.409pC/mm2.The problem is represented by7×4 circular patches and 25×7 regular nodes as shown in Figure 15. The patch fill distance is about H=5.341 and the coefficient is αr=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 wA and the first principal stress σ1B are in very good accord with the results given by the FEM. Although the errors of the electric potential ϕA and the electric flux density DC are slightly lager, the present method still achieves reasonable prediction.Table 6
Numerical results of piezoelectric Cook’s membrane and relative errors for different nodal densities.
Background cell Node density w A × 1 0 - 4 (mm) ϕ A × 1 0 - 8 (GV) σ 1 B (N/mm2) D c (pC/mm2) 12 × 6 13 × 7 2.1132 1.7104 0.21582 23.2896 ( 0.199 % ) ( - 1.247 % ) ( - 0.143 % ) ( 3.930 % ) 24 × 6 25 × 7 2.1099 1.6401 0.21764 23.1817 ( 0.043 % ) ( - 5.306 % ) ( 0.699 % ) ( 3.448 % ) 36 × 6 37 × 7 2.1234 1.6173 0.21763 23.2457 ( 0.683 % ) ( - 6.622 % ) ( 0.694 % ) ( 3.734 % ) 48 × 6 49 × 7 2.1342 1.6306 0.21718 22.9610 ( 1.195 % ) ( - 5.855 % ) ( 0.486 % ) ( 2.463 % ) 12 × 12 13 × 13 2.1431 1.7593 0.21740 24.1962 ( 1.617 % ) ( 1.576 % ) ( 0.588 % ) ( 7.975 % ) 24 × 12 25 × 13 2.1356 1.7779 0.21637 23.9918 ( 1.261 % ) ( 2.650 % ) ( 0.111 % ) ( 7.063 % ) 36 × 12 37 × 13 2.1370 1.7935 0.21600 24.0339 ( 1.328 % ) ( 3.551 % ) ( - 0.060 % ) ( 7.251 % ) 48 × 12 49 × 13 2.1377 1.7944 0.21612 24.0447 ( 1.361 % ) ( 3.603 % ) ( - 0.005 % ) ( 7.299 % ) Nguyen-Van et al. [8] 2.100 1.703 0.2109 20.754 ( - 0.379 % ) ( - 1.674 % ) ( - 2.406 % ) ( - 7.385 % ) Long et al. [42] 2.109 1.732 0.21613 22.409Figure 15
Node distribution and patches.


### 4.5. 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 Figure16. Each bimorph has length L=10 μm and height h=0.5 μm and the length of the mirror is λ=1 μm. The two-layer piezoelectric bimorph is made of PVDF material and its material properties are summarized as follows:c11:2.18×10-3 N/μm2c13:6.33×10-4 N/μm2c33:2.18×10-3 N/μm2c55:7.75×10-4 N/μm2e31:4.6×10-8 N/Vμme33:4.6×10-8 N/Vμmκ11:1.062×10-10 N/V2κ33:1.062×10-10 N/V2Figure 16
A piezoelectric bimorph optical microscanner device.When a voltage is applied on the surfaces, the bimorphs deflect in opposite direction so that the mirror moves up and down, resulting in the change of the tilt angleβ. The direction of the reflected light can be changed with applied different voltages and the tip deflection δ can be calculated by the present method. For sufficient small rotation, the angle β can be approximated from the tip deflection(46)β=2δλ.The boundary conditions of the bottom layer are(47)ϕ1x,-h=V0,σz1x,-h=0,τxz1x,-h=0,ϕ1x,0=0,σz1x,0=σz2x,0,τxz1x,0=τxz2x,0,ϕ,x10,z=0,u10,z=0,w10,z=0,ϕ,x1L,z=0,σx1L,z=0,τxz1L,z=0.And the boundary conditions of the top layer are(48)ϕ2x,0=ϕ1x,0,u2x,0=u1x,0,w2x,0=w1x,0,ϕ2x,h=V0,σz2x,h=0,τxz2x,h=0,ϕ,x20,z=0,u20,z=0,w20,z=0,ϕ,x2L,z=0,σx2L,z=0,τxz2L,z=0.Considering the geometry of the optical microscanner device, node arrangements and patches are very similar to the parallel piezoelectric bimorph beam. The central points are at the upper and down quarter through the thickness direction.21×4 patches and αr=1.2 are used for this problem.The tip displacements calculated for several applied voltages with different nodal densities are listed in Table7, 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.Table 7
Numerical results of piezoelectric optical microscanner device for different nodal densities (μm).
Applied voltage (V) 40 × 4 50 × 4 60 × 4 80 × 4 100 × 4 PCM [43] 41 × 5 51 × 5 61 × 5 81 × 5 101 × 5 1.00 4.8742 E - 03 4.8766 E - 03 4.8631 E - 03 4.8754 E - 03 4.8716 E - 03 4.9360 E - 03 2.00 9.7484 E - 03 9.7533 E - 03 9.7262 E - 03 9.7508 E - 03 9.7431 E - 03 9.8720 E - 03 5.00 2.4371 E - 02 2.4383 E - 02 2.4316 E - 02 2.4377 E - 02 2.4358 E - 02 2.4681 E - 02 10.00 4.8742 E - 02 4.8766 E - 02 4.8631 E - 02 4.8754 E - 02 4.8716 E - 02 4.9362 E - 02 15.00 7.3113 E - 02 7.3149 E - 02 7.2947 E - 02 7.3131 E - 02 7.3074 E - 02 7.4043 E - 02 20.00 9.7484 E - 02 9.7533 E - 02 9.7262 E - 02 9.7508 E - 02 9.7431 E - 02 9.8724 E - 02 25.00 1.2186 E - 01 1.2192 E - 01 1.2158 E - 01 1.2188 E - 01 1.2179 E - 01 1.2341 E - 01 50.00 2.4371 E - 01 2.4383 E - 01 2.4316 E - 01 2.4377 E - 01 2.4358 E - 01 2.4681 E - 01From 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.Figure 17
The variation of the angle under different applied voltages.


## 4.1. Shear Deformation of Piezoelectric Strip

The shear deformation of a square piezoelectric strip polarized inz-direction is taken into consideration. The piezoelectric strip is made of material PZT-5 and its material properties are given as follows:s11:16.4×10-6 (mm)2/Ns13:-7.22×10-6 (mm)2/Ns33:18.8×10-6 (mm)2/Ns55:47.5×10-6 (mm)2/Nσ0:-5.0 N/(mm)2σ1:20.0 N/(mm)2L:1.0 mmd31:-172×10-9 mm/Vd33:374×10-9 mm/Vd15:584×10-9 mm/Vξ11:1.53105×10-8 N/V2ξ33:1.505×10-7 N/V2V0:1000 Vh:0.5 mmThe strip is subjected a uniform loadσ0 on the top and bottom boundaries and an applied voltage V0 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.Figure 4
Piezoelectric strip subjected to a uniform load and an applied voltage.The following boundary conditions of the piezoelectric strip are given by(40)ϕ,zx,±h=0,σzx,±h=σ0,τxzL,z=0,τxzx,±h=0,ϕL,z=-V0,σxL,z=0,ϕ0,z=+V0,u0,z=0,w0,0=0.The analytical solutions for this shear problem are given by [49](41)u=s13σ0,w=d15V0xh+s33σ0z,ϕ=V01-2xL.The problem domain is represented by both11×11 regular nodes and 121 irregular nodes in Figure 5. The PU cover consists of M=m2(5×5) regular circular patches and the coefficient αr=1.2 is taken here. 10×10 square background cells are used for numerical integration.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.
(a) (b)The obtained results of the present method are plotted in Figure6, in which symbols R and IR denote regular and irregular nodes, respectively. Figure 6(a) shows the horizontal displacement u along z=0 and the deflection w along z=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.Figure 6
Comparisons of the present solutions, the RBF solutions, and the analytical solutions of the piezoelectric strip in shear: (a) the variation of the horizontal displacementu along z=0; (b) the variation of the vertical displacement w along z=0; (c) the variation of the electric potential ϕ along z=0.
(a) (b) (c)The computed mechanical deformation results of the present method for irregular nodes are plotted in Figure7, which seems quite reasonable for the shear deformation. The convergence curves and convergence rates of displacements u and w 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.Figure 7
The shear deformation of the piezoelectric strip for irregular nodes. Note that the displacements plotted are magnified 100 times.Figure 8
Convergence property of displacementsu and w and electric potential ϕ for regular nodes.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 pointA (x=3L/4, z=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 αr may change with the patch fill distance and the coefficient αr can be adjusted according to the patch fill distance.Table 1
Numerical results of displacementsuA×10-5 (mm) and wA×10-4 (mm) and electric potential ϕA(V) at point A(x=3L/4,z=0) for different patch sizes.
α r 1.1 1.2 1.3 1.4 1.5 1.6 1.8 2.0 Exact u A R 2.70774 2.70762 2.70747 2.70740 2.70753 2.70752 2.70752 2.70752 2.70750 IR 2.70523 2.70372 2.70676 2.70783 2.70844 2.70735 2.70760 2.70804 w A R 8.76021 8.76018 8.76018 8.76025 8.76045 8.76068 8.76078 8.76089 8.76000 IR 8.76367 8.76041 8.75980 8.75963 8.76029 8.75998 8.75991 8.75998 ϕ A R −499.92 −499.92 −499.93 −499.94 −500.01 −500.00 −500.02 −500.03 −500.00 IR −500.79 −500.05 −499.96 −499.94 −500.03 −499.97 −499.98 −499.98


## 4.2. 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 Figure9. The analytical solutions of displacements and electric potential are given by [50](42)u=s11σ0-d31V0s11hx+s111-d312s11ξ33σ1xz,w=s13σ0-d31V0s11hz+s131-d33d31s13ξ33-s111-d312s11ξ33σ1x22,ϕ=V0zh-d31σ12ξ33h2-z2,and the boundary conditions for this problem are given by(43)ϕ,zx,±h=±V0,σzx,±h=0,τxzx,z±h=0,ϕ,xL,z=0,σxL,z=σ0+σ1z,τxzL,z=0,ϕ,x0,z=0,u0,z=0,w0,0=0.Figure 9
Piezoelectric strip subjected to a linear stress and an applied voltage.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 Figure10, in which the patch fill distance is about H=0.1938 and the coefficient αr=1.2 is also taken here.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.
(a) (b)The solutions of the RBF-PUM, the RBF method, and the analytical method are shown in Figure11. The displacement u and electric potential ϕ along x=L are plotted in Figures 11(a) and 11(c), respectively. The variation of the deflection w 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.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 displacementu along x=L; (b) the variation of the vertical displacement w along z=-h; (c) the variation of the electric potential ϕ along x=L.
(a) (b) (c)


## 4.3. A Parallel Piezoelectric Bimorph Beam

Consider a two-layer parallel piezoelectric bimorph beam as shown in Figure12. The beam is made of material PVDF with the same thickness and polarization direction and the material properties are summarized as follows:E:2.0×109 Paν:0.29e31:0.046 C/m2e32:0.046 C/m2κ11:0.1062×10-9 F/mκ33:0.1062×10-9 F/mFigure 12
A parallel piezoelectric bimorph beam.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 an applied voltageV, the tip deflection δ can be approximated as [51](44)δ=2L2Vd31H2,where L is the length of the beam, H is the total thickness, and d31 is the piezoelectric strain constant. The size of the beam is taken as L=5 mm and H=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 of21×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 αr=1.2 is taken.Figure 13
Node distribution and patches of the top beam.The tip deflections for different node densities are listed in Table2 and the present results are very close to the results given by the ABAQUS [19]. The effect of the patch size on 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 q and η are also studied in Tables 4 and 5, respectively. It can be observed that increasing the parameter q 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].Table 2
The tip deflection of the piezoelectric beam for different nodes whenq=1.95, η=2.0, αr=1.2 (×10-8 m).
Background cell 40 × 4 50 × 4 60 × 4 70 × 4 80 × 4 100 × 4 CPS4E [19] Nodal density 41 × 5 51 × 5 61 × 5 71 × 5 81 × 5 101 × 5  Tip deflection 1.0715 1.0719 1.0705 1.0737 1.0736 1.0719 1.073Table 3
Effects of the patch size on the tip deflection for51×5 nodes when q=1.95, η=2.0 (×10-8 m).
α r 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8  Tip deflection 1.0494 1.0719 1.0721 1.0725 1.0726 1.0738 1.0738 1.0738Table 4
Effects of parameterq on the tip deflection for 51×5 nodes when η=2.0, αr=1.2 (×10-8 m).
q 0.95 0.98 1.03 1.05 1.08 1.92 1.95 1.98  Tip deflection 1.0709 1.0709 1.0710 1.0710 1.0710 1.0719 1.0719 1.0719Table 5
Effects of parameterη on the tip deflection for 51×5 nodes when q=1.95, αr=1.2 (×10-8 m).
η 0.5 0.75 1.0 1.25 1.5 2.0 2.5 3.0  Tip deflection 1.0719 1.0719 1.0719 1.0719 1.0719 1.0719 1.0719 1.0719


## 4.4. Piezoelectric Cook’s Membrane

The piezoelectric Cook’s membrane problem is solved using the present method. The geometry and applied unity tip loadF 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:c11:139×103 N/mm2c13:74.3×103 N/mm2c33:113×103 N/mm2c55:25.6×103 N/mm2e15:13.44×106 pC/mm2e31:-6.98×106 pC/mm2e33:13.84×106 pC/mm2κ11:6.00×109 pC/GVmmκ33:5.47×109 pC/GVmmFigure 14
Cook’s membrane.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 A, B, C are (45)wA=2.109×10-4mm,ϕA=1.732×10-8GVσ1B=0.21613N/mm2,DC=22.409pC/mm2.The problem is represented by7×4 circular patches and 25×7 regular nodes as shown in Figure 15. The patch fill distance is about H=5.341 and the coefficient is αr=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 wA and the first principal stress σ1B are in very good accord with the results given by the FEM. Although the errors of the electric potential ϕA and the electric flux density DC are slightly lager, the present method still achieves reasonable prediction.Table 6
Numerical results of piezoelectric Cook’s membrane and relative errors for different nodal densities.
Background cell Node density w A × 1 0 - 4 (mm) ϕ A × 1 0 - 8 (GV) σ 1 B (N/mm2) D c (pC/mm2) 12 × 6 13 × 7 2.1132 1.7104 0.21582 23.2896 ( 0.199 % ) ( - 1.247 % ) ( - 0.143 % ) ( 3.930 % ) 24 × 6 25 × 7 2.1099 1.6401 0.21764 23.1817 ( 0.043 % ) ( - 5.306 % ) ( 0.699 % ) ( 3.448 % ) 36 × 6 37 × 7 2.1234 1.6173 0.21763 23.2457 ( 0.683 % ) ( - 6.622 % ) ( 0.694 % ) ( 3.734 % ) 48 × 6 49 × 7 2.1342 1.6306 0.21718 22.9610 ( 1.195 % ) ( - 5.855 % ) ( 0.486 % ) ( 2.463 % ) 12 × 12 13 × 13 2.1431 1.7593 0.21740 24.1962 ( 1.617 % ) ( 1.576 % ) ( 0.588 % ) ( 7.975 % ) 24 × 12 25 × 13 2.1356 1.7779 0.21637 23.9918 ( 1.261 % ) ( 2.650 % ) ( 0.111 % ) ( 7.063 % ) 36 × 12 37 × 13 2.1370 1.7935 0.21600 24.0339 ( 1.328 % ) ( 3.551 % ) ( - 0.060 % ) ( 7.251 % ) 48 × 12 49 × 13 2.1377 1.7944 0.21612 24.0447 ( 1.361 % ) ( 3.603 % ) ( - 0.005 % ) ( 7.299 % ) Nguyen-Van et al. [8] 2.100 1.703 0.2109 20.754 ( - 0.379 % ) ( - 1.674 % ) ( - 2.406 % ) ( - 7.385 % ) Long et al. [42] 2.109 1.732 0.21613 22.409Figure 15
Node distribution and patches.


## 4.5. 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 Figure16. Each bimorph has length L=10 μm and height h=0.5 μm and the length of the mirror is λ=1 μm. The two-layer piezoelectric bimorph is made of PVDF material and its material properties are summarized as follows:c11:2.18×10-3 N/μm2c13:6.33×10-4 N/μm2c33:2.18×10-3 N/μm2c55:7.75×10-4 N/μm2e31:4.6×10-8 N/Vμme33:4.6×10-8 N/Vμmκ11:1.062×10-10 N/V2κ33:1.062×10-10 N/V2Figure 16
A piezoelectric bimorph optical microscanner device.When a voltage is applied on the surfaces, the bimorphs deflect in opposite direction so that the mirror moves up and down, resulting in the change of the tilt angleβ. The direction of the reflected light can be changed with applied different voltages and the tip deflection δ can be calculated by the present method. For sufficient small rotation, the angle β can be approximated from the tip deflection(46)β=2δλ.The boundary conditions of the bottom layer are(47)ϕ1x,-h=V0,σz1x,-h=0,τxz1x,-h=0,ϕ1x,0=0,σz1x,0=σz2x,0,τxz1x,0=τxz2x,0,ϕ,x10,z=0,u10,z=0,w10,z=0,ϕ,x1L,z=0,σx1L,z=0,τxz1L,z=0.And the boundary conditions of the top layer are(48)ϕ2x,0=ϕ1x,0,u2x,0=u1x,0,w2x,0=w1x,0,ϕ2x,h=V0,σz2x,h=0,τxz2x,h=0,ϕ,x20,z=0,u20,z=0,w20,z=0,ϕ,x2L,z=0,σx2L,z=0,τxz2L,z=0.Considering the geometry of the optical microscanner device, node arrangements and patches are very similar to the parallel piezoelectric bimorph beam. The central points are at the upper and down quarter through the thickness direction.21×4 patches and αr=1.2 are used for this problem.The tip displacements calculated for several applied voltages with different nodal densities are listed in Table7, 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.Table 7
Numerical results of piezoelectric optical microscanner device for different nodal densities (μm).
Applied voltage (V) 40 × 4 50 × 4 60 × 4 80 × 4 100 × 4 PCM [43] 41 × 5 51 × 5 61 × 5 81 × 5 101 × 5 1.00 4.8742 E - 03 4.8766 E - 03 4.8631 E - 03 4.8754 E - 03 4.8716 E - 03 4.9360 E - 03 2.00 9.7484 E - 03 9.7533 E - 03 9.7262 E - 03 9.7508 E - 03 9.7431 E - 03 9.8720 E - 03 5.00 2.4371 E - 02 2.4383 E - 02 2.4316 E - 02 2.4377 E - 02 2.4358 E - 02 2.4681 E - 02 10.00 4.8742 E - 02 4.8766 E - 02 4.8631 E - 02 4.8754 E - 02 4.8716 E - 02 4.9362 E - 02 15.00 7.3113 E - 02 7.3149 E - 02 7.2947 E - 02 7.3131 E - 02 7.3074 E - 02 7.4043 E - 02 20.00 9.7484 E - 02 9.7533 E - 02 9.7262 E - 02 9.7508 E - 02 9.7431 E - 02 9.8724 E - 02 25.00 1.2186 E - 01 1.2192 E - 01 1.2158 E - 01 1.2188 E - 01 1.2179 E - 01 1.2341 E - 01 50.00 2.4371 E - 01 2.4383 E - 01 2.4316 E - 01 2.4377 E - 01 2.4358 E - 01 2.4681 E - 01From 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.Figure 17
The variation of the angle under different applied voltages.


## 5. 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.



---

*Source: 7632176-2016-12-12.xml*