On the consistency of two-phase local/nonlocal piezoelectric integral model

In this paper, we propose general strain- and stress-driven two-phase local/nonlocal piezoelectric integral models, which can distinguish the difference of nonlocal effects on the elastic and piezoelectric behaviors of nanostructures. The nonlocal piezoelectric model is transformed from integral to an equivalent differential form with four constitutive boundary conditions due to the difficulty in solving intergro-differential equations directly. The nonlocal piezoelectric integral models are used to model the static bending of the Euler-Bernoulli piezoelectric beam on the assumption that the nonlocal elastic and piezoelectric parameters are coincident with each other. The governing differential equations as well as constitutive and standard boundary conditions are deduced. It is found that purely strain- and stress-driven nonlocal piezoelectric integral models are ill-posed, because the total number of differential orders for governing equations is less than that of boundary conditions. Meanwhile, the traditional nonlocal piezoelectric differential model would lead to inconsistent bending response for Euler-Bernoulli piezoelectric beam under different boundary and loading conditions. Several nominal variables are introduced to normalize the governing equations and boundary conditions, and the general differential quadrature method (GDQM) is used to obtain the numerical solutions. The results from current models are validated against results in the literature. It is clearly established that a consistent softening and toughening effects can be obtained for static bending of the Euler-Bernoulli beam based on the general strain- and stress-driven local/nonlocal piezoelectric integral models, respectively.


Introduction
Since Pan and his co-authors [1] first reported ZnO piezoelectric nanowires in 2001, piezoelectric nanomaterials and nanostructures have received close attention from the modern sci-ence and technology. Piezoelectric nanostructures were widely used in chemical sensors [2] , micro/nano electromechanical [3] , nanogenerators [4] , piezoelectric field effect transistor [5] and so on due to the outstanding mechanical and electric properties. The experimental tests and atomic simulations [6][7] showed that the piezoelectric nanostructures have a strong size effect at the nanoscale and exhibit different behaviors from the macro-scale.
However, several studies showed that there were some inconsistent results while using the nonlocal elastic differential model to study the bar under a tensile load [27] , the cantilever beam under a concentrated load at the end [28][29][30] , and the clamped beams under a uniformly distributed load [31] . Therefore, the original nonlocal model in integral form restarts to attract the attention from the scientific community. Romano et al. [32] found that the transformation for nonlocal elastic model from integral to differential forms should be equipped with two extra constraints from integral constitutive relation, and further pointed out the purely nonlocal integral model would lead to an ill-posed problem, since there was a conflict between the two constitutive constraints and the equilibrium equations [33] .
Wang et al. [34][35] transformed the two-phase local/nonlocal integral model proposed initially by Eringen [36] into the equivalent differential form with two constitutive constraints and used it to study the static bending of straight Euler-Bernoulli and Timoshenko beams under different boundary and loading conditions, respectively. With the help of Laplace transform, Zhang et al. [37] and Zhang and Qing [38] utilized the two-phase local/nonlocal integral model directly to study the static bending of curved Euler-Bernoulli and Timoshenko beams. Meng et al. [39] deduced the semi-analytic solution for static bending Euler-Bernoulli beam with the axial force. A consistent softening effect was obtained while using the two-phase strain-driven local/nonlocal integral model.
In order to overcome the ill-posedness of the purely strain-driven nonlocal integral model, Romano and Barretta [33] proposed a stress-driven nonlocal integral model by swapping the interpretation of source and output fields in Eringen's nonlocal integral model. With the application of stress-derived nonlocal integral model, Oskouie et al. [40][41] utilized the Ritz-method to bypass constitutive constraints and studied the bending response of Euler-Bernoulli and Timoshenko beams. Based on the Laplace transform, Zhang et al. [42] deduced the exact and asymptotic bending analysis of straight Euler-Bernoulli and Timoshenko beams, and Zhang et al. [43] derived the exact solutions for bending of curved functionally graded Timoshenko beam. Barretta et al. [44][45] used the stress-driven based purely nonlocal and two-phase local/nonlocal integral models to study the static bending of Timoshenko beam, respectively. A consistent toughening effect was obtained while using two-phase stress-driven local/nonlocal integral model.
To the best knowledge of the authors, there is no literature about the application of the nonlocal integral model to the nanoscale piezoelectric structures. The plan for this paper is as follows. In Section 2, we present the Helmholtz-kernel based two-phase local/nonlocal piezoelectric integral models, and transform it from integral equivalently differential forms with constitutive boundary conditions. In Section 3, the differential governing equation and standard boundary conditions for a piezoelectric Euler-Bernoulli beam are derived, and the equivalently nonlocal differential constitutive relations with constitutive boundary conditions are presented. In Section 4, nominal variables are introduced to normalize the governing equations as well as constitutive and standard boundary conditions, and the general differential quadrature method is used to solve the differential equations numerically. In Section 5, several benchmark exam-ples for the nonlocal piezoelectric Euler-Bernoulli beam under different boundary and loading conditions are illustrated numerically and discussed.
2 Two-phase local/nonlocal piezoelectric integral model

Nonlocal piezoelectric integral model
According to the classic piezoelectric model [46] , the constitutive relations in stress-charge and strain-voltage forms can be expressed respectively as where S ij and ∆ i are components for the local stress tensor and the electric displacement vector, respectively. c ijkl , e kij , and µ ik are the elastic, piezoelectric, and dielectric constants, respectively. s ijkl , g kij , and β ik are the elastic compliance constants, alternative forms of the piezoelectric constant, and the permittivity constants, respectively. ε kl and E i are components of strain tensor and electric filed vector, which are defined as where u i and ϕ are the displacement component and electric potential, respectively. From the idea of Romano and Barretta [33] , one-dimensional general strain-and stress-driven two-phase local/nonlocal piezoelectric integral models are expressed respectively as and where σ ij and D i are components for the nonlocal stress tensor and the electric displacement vector, and κ i is the nonlocal parameter. ξ i is the volume fraction of nonlocal effect. It degenerates into the classic piezoelectric model when ξ i = 0, and the model becomes a purely nonlocal integral model when ξ i = 1.

2.2
Transformation of nonlocal piezoelectric model from integral to differential forms It is well-known that the integro-differential equations are generally more difficult to solve than the differential equations. The nonlocal integral equations (4) and (5) are transformed into equivalent differential forms with corresponding constitutive boundary conditions in the following.
Proposition 1 The nonlocal piezoelectric integral model can be generally expressed as which can be transformed into an equivalent differential equation as with constitutive boundary conditions at x = a and in which (·) = d(·) dx . Proof Removing the absolute sign in Eq. (6), we get Performing the first, the second, the third, and the fourth derivatives of Eq. (10), we obtain Combination of Eqs. (6), (10), and (12) gives differential equation (6).
3 Static bending analysis of nonlocal piezoelectric Euler-Bernoulli nanobeam 3.1 Governing equation of Euler-Bernoulli beam Figure 1 shows a piezoelectric nanobeam with length L and cross-sectional area A = bh (thickness h and width b) under the action of a distributed load q(x) as well as the concentrated force Q and the bending moment M . The polarization direction of the piezoelectric structure is parallel to the positive z-axis. The displacement field for an Euler-Bernoulli beam theory is

Fig. 1 Schematic illustration of a piezoelectric nanobeam
According to Wang [47] , the electrical potential of piezoelectric beams, which should satisfy the Maxwell equation, can be assumed to be where β = π/h, and φ(x) is the electric potential in the x-direction. Combining Eqs. (3) with (22)- (23), we express the nonzero components for the strain and electric filed as Combination of Eqs. (1) and (24) gives Therefore, the variation of the strain energy of the piezoelectric beams is given as Substituting Eq. (24) into Eq. (26), after a lengthy simplification, we get where The variation of external work δW e due to external load is where q is the distributed load, and M and Q are the bending moment and force applied at beam ends, respectively. According to the minimum total potential energy principle, the government equations and the standard boundary conditions can be expressed as 3.2 Two-phase strain-driven local/nonlocal piezoelectric integral model Based on simplified strain-driven local/nonlocal piezoelectric integral model, we can express relation between nonlocal stress/electric displacement and strain and electric filed as Combining Eqs. (28) and (32) where According to Eqs. (19)-(21), the integral constitutive equations can be transformed equivalently into differential forms with the corresponding constitutive constraints Combining Eqs. (30a) with (36a) as well as Eqs. (30b) with (35a), we get Substituting P and M from Eqs.

Mathematical formulation in normalized form
The following nominal quantities are introduced to normalize the differential equations and boundary conditions as well as constitutive constraints: For the strain-driven nonlocal piezoelectric integral model, the differential governing equations (37a) and (40) where (·) η = d(·) dη , and The constitutive constraints are normalized as By combining Eqs. (31) with (38) and (39), the normalized standard boundary conditions can be expressed as where Q = L 2 EI1 Q, and M = L EI1 M .
For the stress-driven nonlocal piezoelectric integral model, the differential governing equations (37a) and (40) The corresponding normalized constitutive and standard boundary conditions are (55) It should be noticed that for both strain-and stress-driven two-phase local/nonlocal piezoelectric integral models, Eqs. (50) and (54) are about sixth-order W , fourth-order Φ, and secondorder R ordinary differential equations. There are totally 12 boundary conditions (six standard boundary conditions and six constitutive boundary conditions), which have the same number of total orders for the differential governing equations. Therefore, the static bending of piezoelectric Euler-Bernoulli beam is well-posed for the two-phase local/nonlocal piezoelectric integral model.
However, if ξ = 1, Eq. (50) turns into fourth-order W , second-order Φ, and second-order R ordinary differential equations for strain-driven purely nonlocal piezoelectric integral model, and Eq. (54) turns into sixth-order W and fourth-order Φ ordinary differential equations for the stress-driven purely nonlocal piezoelectric integral model. Therefore, the total number of differential orders for governing equations is less than that of boundary conditions. Therefore, the static bending of piezoelectric Euler-Bernoulli beam is ill-posed for the purely stress-driven nonlocal piezoelectric integral models.
Special attention is paid to the traditional nonlocal piezoelectric differential model when we set ξ = 1 and neglect the constitutive boundary conditions, in which the total number of differential orders for governing equations is the same as the number of standard boundary conditions. However, it can be seen that the differential governing equations are independent of the nonlocal parameter for a concentrated force and a constantly or linearly distributed load. Taking into account the boundary conditions, we obtain that no size-dependent bending appears while using traditional nonlocal piezoelectric differential model under following cases.
(iii) Concentrated force at beam end: clamped-free, clamped-guided, and simply supportedguided boundary conditions.
In other words, the traditional nonlocal piezoelectric differential model would lead to inconsistent bending response for Euler-Bernoulli piezoelectric beam under different boundary and loading conditions.

Numerical solution using GDQM
The GDQM shows great advantages in solving differential equations, which is used to solve Eqs. (50) and (54). On the basis of the GDQM and following the Chebyshev-Gauss-Lobatto rule, the beam domain is discretized with n grid nodes as where N is the grid number. The functions R, Φ, and W can be approximated by [48][49][50] where ηi−η k is the Lagrange interpolation function, and where j = 1, N . Performing derivative with respect to η on Eqs. (58)-(60), we obtain where X (i) , Y (i) , and Z (i) are the weighting coefficients of the ith-order derivative, which are defined explicitly [47][48][49] as Following the standard GDQM procedure, we obtain the discrete form of the governing differential equations as well as approximate constitutive and standard boundary conditions. The static bending problem for a piezoelectric Euler-Bernoulli beam based on nonlocal integral model turns a set of 3N +6 linear equations about ∆ R j , ∆ Φ j , and ∆ W j .

Numerical results
In this section, the effect of nonlocal parameters ξ and κ on deflections of piezoelectric Euler-Bernoulli beam with L/h = 20 is investigated on the basis of strain-and stress-driven nonlocal integral models. It is assumed that the piezoelectric beam is made of the PZT-4 with c 1111 = 132 GPa, e 311 = −4.1 C · m −2 , µ 11 = 5.841 × 10 −9 C · (V · m) −1 , and µ 33 = 7.124 × 10 −9 C · (V · m) −1 [51] . tributed load (Q = 1), respectively, and CFC, CGC, and SGC indicate clamped-free, clampedguided, and simply-supported-guided beams under a unit concentrated force ( Q = 1), respectively. In Figs. 2-3, the solid lines represent the results from current models and the scatters indicate the results for corresponding nonlocal Euler-Bernoulli beams from Wang et al. [34] and Zhang et al. [42] . It can be seen from Figs. 2-3 that the results based on current model agree well with those based on nonlocal integral elastic model while piezoelectric and dielectric constants are very small. In other words, the nonlocal piezoelectric integral model can be degenerated to the nonlocal integral elastic model while setting piezoelectric and dielectric constants as zeros.   Figures 6-9 show the contour plots of the normalized maximum deflection (NMD), which is defined as the ratio between maximum deflection from nonlocal piezoelectric model and that from classic piezoelectric model, in which λ and ξ range from 0.05 to 1 and 0 to 0.95, respectively. It can be seen that, under different boundary and loading conditions, a consistent softening and toughening response can be observed for strain-and stress-driven nonlocal piezoelectric models, respectively. Meanwhile, the NMDs increase and decrease consistently with the increase of λ and ξ for strain-and stress-driven nonlocal piezoelectric models, respectively. When ξ is small, both the increase and decrease rates of NMDs decrease with the increase in λ for strain-and stress-driven nonlocal piezoelectric integral models, respectively, which is also shown in Fig. 4(a).

Effect of nonlocal parameters on Euler-Bernoulli beam bending
When λ is small, the increase and decrease rates of NMDs decrease with the increase of ξ for strain-and stress-driven nonlocal piezoelectric integral models, respectively. In addition, along the line λ = ξ, the increase and decrease rates of NMDs increase with the increase in λ for both strain-and stress-driven nonlocal piezoelectric integral models.

Conclusions
General strain-and stress-driven two-phase local/nonlocal piezoelectric integral models are proposed in this paper, which can distinguish the nonlocal effect on elastic and piezoelectric behaviors. The integral nonlocal model is transformed equivalently into differential form with four constitutive boundary conditions. The simplified strain-and stress-driven local/nonlocal piezoelectric integral models can be obtained when the nonlocal parameters for elastic and piezoelectric behaviors are the same. On the basis of the minimum total potential energy principle, the governing differential equations as well as standard and constitutive boundary conditions are derived for static bending of Euler-Bernoulli beam for general strain-and stressdriven two-phase local/nonlocal piezoelectric integral models.
It is found that the traditional piezoelectric differential model, purely strain-and stressdriven nonlocal piezoelectric integral models are not applicable for static bending analysis of Euler-Bernoulli beam, because that the numbers of differential order of governing equations and boundary conditions are different. Several nominal variables are introduced to normalize the expression for governing equations and boundary conditions, and the GDQM with different order interpolations are used to obtain the numerical solution.
The numerical study shows that, by setting e 311 , µ 11 , and µ 33 as very small value, the strain-driven local/nonlocal piezoelectric integral model can be degenerated to strain-driven local/nonlocal integral elastic model, and stress-driven local/nonlocal piezoelectric integral model can be generated to purely stress-driven nonlocal elastic model for ξ = 0. A consistent softening and toughening effect can be obtained for strain-and stress-driven local/nonlocal piezoelectric integral models, respectively. For a fixed value of λ(ξ), both the increase and decrease rates of NMDs decrease with the increase in ξ(λ) for strain-and stress-driven nonlocal piezoelectric integral models, respectively.