A meshfree formulation for large deformation analysis of flexoelectric structures accounting for the surface effects
AA meshfree formulation for large deformation analysis of flexoelectricstructures accounting for the surface effects
Xiaoying Zhuang a,b, ∗ , S.S.Nanthakumar b , Timon Rabczuk c, ∗∗ a College of Civil Engineering, Tongji University, 1239 Siping Road, 200092 Shanghai, China. b Institute of Continnum Mechanics, Leibniz University Hannover, Appelstrasse 11A, D-30167 Hannover Germany. c Institute of Structural Mechanics, Bauhaus University Weimar, Marienstrasse 5, 99423 Weimar, Germany.
Abstract
In this work, we present a compactly supported radial basis function (CSRBF) based meshfree method toanalyse geometrically nonlinear flexoelectric nanostructures considering surface effects. Flexoelectricity isthe polarization of dielectric materials due to the gradient of strain, which is different from piezoelectricityin which polarization is dependent linearly on strain. The surface effects gain prominence as the size of thestructure tends to nanoscale and so their consideration is inevitable when flexoelectric nanostructures areanalysed. First, the proposed meshfree formulation is validated and the influence of nonlinear strain terms onthe energy conversion ability of flexoelectric beams made of a non-piezoelectric material like cubic StrontiumTitanate is studied. Subsequently, the meshfree formulation for nonlinear flexoelectricity is extended toinclude nonlinear surface effects. It is determined that the surface effects can have notable influence on theoutput flexoelectric voltage of nano-sized cantilever structures in the nonlinear regime.
Keywords:
Meshfree method; Nolinear Flexoelectricity; Geometric nonlinearity; Surface effects
1. Introduction
Flexoelectricity is the generation of electric polarization under mechanical strain gradient or mechanicaldeformation due to electric field gradient (converse flexo). It is a more general phenomenon than the linearchange in polarization due to stress, known as the piezoelectric effect. Flexoelectric polarization is restrictednot only to non-centrosymmetric crystals eventually opening up possibilities for nontoxic electromechanicalmaterials for biomedical application.Piezoelectricity can be characterised by a third rank tensor and is observed only in non-centrosymmetriccrystals (21 types). In contrast, flexoelectricity can be mathematically defined by a fourth order tensor andcan be observed in materials of any symmetry (32). The reason behind is that the homogeneous strain relieson lack of symmetry of materials for polarization, on the other hand, the strain gradient can break the localcentrosymmetry of materials inducing polarization. The strain gradient scales with the size of specimenleading to the possibility of significant flexoelectric effect at the length scale of nanometers. Piezoelectricityexists only below Curie temperature, while flexoelectricity being symmetry independent does not have atemperature constraint [1]. The high energy conversion ability of piezoelectric materials makes them theprominent constituent in several micro [2] and nano-sized [3] energy harvesters developed. While, the recent ∗ Corresponding author ∗∗ Corresponding author
Email addresses: [email protected] (Xiaoying Zhuang), [email protected] (Timon Rabczuk)
Preprint submitted to Elsevier November 13, 2019 esearches show the possibility of an energy harvester made of non-piezoelectric materials exploiting flexoelec-tricity [4]. Nanoelectromechanical systems like actuators [5] are fabricated using non-piezoelectric materialslike Strontium Titanate and are shown to produce curvature/electric field ratio of 3.33
M V − comparableto the ratio of 5.2 M V − in piezoelectric Lead Zirconium Titanate bimorph. The flexoelectricity also offersthe advantage of choosing Lead-free materials as constituents in sensors, actuators and energy harvesters.The theory of flexoelectricity was first identified way back in the 1960s by Mashkevich and Tolpygo [6],followed later by the work of Tagantsev [7] in which bulk and surface mechanisms that can cause polarizationdue to strain gradient were determined. Meanwhile, Kogan [8] made a theoretical estimate of the flexoelectriccoefficient to be of the order of e/a ≈ − C/m , where e is the electronic charge and a is the lattice parameter.The series of experimental works by Cross and co-workers [9, 10, 11] sparked interest over the potential offlexoelectric materials as a substitute to piezoelectric materials. These experimental studies on ceramicswith cubic symmetry like Barium Strontium Titanate (BST) and Barium Titanate (BTO) revealed highervalues of peak flexoelectric coefficients in the range of 50 µC/m . Atomistically, Maranganti et al . [12]determined the flexoelectric coefficients of several ferroelectric and non-ferroelectric crystals. Though thereis discrepancy between theoretical and experimental flexoelectric coefficients of Barium Titanate (BTO),the theoretical estimations are of the same order compared to the experimental results of Zubko et al . [13]for Strontium Titanate (STO) crystals. Several works are available in literature that presents analyticallyderived electro-elastic field for nanobeams and nanowires having flexoelectric effect and surface effects [14,15, 16]. Nevertheless, the analytical solutions are applicable only to simplified one-dimensional models andso numerical methods to analyse flexoelectric structures are required.Conventional finite element method (FEM) cannot be utilised for analysing flexoelectric structures as thefourth order partial differential equations governing flexoelectricity necessitates C continuity of displacementfield. Phase field modelling of flexoelectricity in an epitaxial thin film made of Barium Titanate is presentedby Chen et al . [17], followed by which analysis of a two phase system is performed [18]. Meshfree shapefunctions offer the advantage of having higher order continuity, making them a favourable class of numericalmethods to analyse flexoelectric structures. A numerical approach to analyse two and three dimensionaltruncated pyramid shaped structure due to flexoelectricity utilising local maximum entropy (LME) meshfreemethod is presented by Abdollahi et al . [19, 20]. Ghasemi et al . [21] proposed an IGA formulation exploitingthe higher order continuity of NURBS shape functions. In [22, 23, 24], mixed FE formulations are proposedfor analysis of two dimensional flexoelectric structures. Though the mixed FE formulation requires only C continuity, the number of nodal DOFs required is much higher. For example, in the flexoelectric elementproposed in [22], degrees of freedom in the corner nodes are two displacement DOFs, four displacementgradient DOFs, one electric potential DOF and four Lagrange multiplier DOFs. It is to be noted that in thework of Nanthakumar et al. [22] flexoelectric nanobeams with surface effects, made of Barium Titanate, areanalysed and optimized using the mixed FE formulation. Also there are computational works available inliterature that specifically analyse nanobeams with surface effects [25, 26, 27], in which extended finite elementmethod is the numerical method adopted. The nonlinear electro-elasticity of soft dielectrics combined withflexoelectricity is analysed by Yvonnet et al . [28], adopting finite element discretization ( C Argyris triangularelements) and consistent linearizations. As a shortcoming the authors have stated that due to instability, theutilised formulation could not simulate the entire nonlinear range.Motivated by all these works on flexoelectricity, in the present work, a compactly supported radial basisfunction (CSRBF) based meshfree formulation is proposed to analyse flexoelectric beams subjected to largedeformation considering surface effects. Though there are works available in literature on analysing flexo-2lectric structures using a meshfree method [19, 20], to the best of our knowledge this is the first work ona meshfree formulation to handle nonlinearity in flexoelectric nanostructures accounting for surface effects.The meshfree shape functions with higher order continuity are advantageous compared to complex mixed FE[22] formulations, mainly because the meshfree formulation requires the discretization of ’only’ displacementand electric potential fields.The outline of the paper is as follows: Section 2 presents the governing equations of flexoelectricity.Section 3 describes the meshfree formulation for flexoelectricity including surface elasticity and surface piezo-electricity. Linearization of the weak form and subsequent meshfree discretization is shown in Section 4.Finally, numerical examples on analysis of two-dimensional flexoelectric structures with surface effects con-sidering geometric nonlinearity are presented in Section 5.
2. Governing equations of flexoelectricity with surface effects
The mathematical modelling of flexoelectricity is based on the extended linear theory of piezoelectricitywith additional strain gradient terms. A general internal energy density function, U involving strain energy,electrostatic energy and terms including strain gradient is presented by Shen et al. [29]. The internal energydensity function, U is as follows, U = U b + U s U b = 12 ε : C : ε − E · e : ε − E · µ ... η − E · κ · E + 12 η ... g ... η U s = U s + τ s : ε s + ω s · E s + 12 ε s : C s : ε s − E s · e s : ε s − E s · κ s · E s (1)where U b and U s are the bulk and surface energy density functions respectively. U s is the surface free energydensity. ε is the linear strain tensor, E is the electric field tensor, ε s and E s are their corresponding surfacecounterparts. τ s and ω s are the residual surface stress and residual surface electric displacements respectively. η is the strain gradient tensor. C and C s are the fourth order bulk and surface stiffness tensors, e and e s arethe third order bulk and surface piezoelectric coupling tensors, κ and κ s are the bulk and surface dielectricpermittivity tensors respectively. g is the sixth order strain gradient elasticity tensor. µ is the fourth orderflexoelectric tensor which represents combination of (a) strain-polarization gradient coupling and (b) straingradient-polarization coupling.The physical stress, σ and electric displacement, D can be obtained from the bulk energy density functionas, σ ij = ∂U b ∂ε ij − (cid:18) ∂U b ∂η ijk (cid:19) ,k = C ijkl ε kl − e ijk E k + µ ijkl E k,l − g ijklmn η lmn,k (2) D i = − ∂U b ∂E i = e ijk ε jk + µ ijkl ε jk,l + κ ij E j (3)The surface mechanical stress, σ s and surface electric displacement, D s can be obtained from the surfaceenergy density function as, σ sij = ∂U s ∂ε sij = τ sij + C sijkl ε skl − e sijk E sk (4) D si = − ∂U s ∂E si = − ω si + e sijk ε sjk + κ sij E sj (5)3he total potential energy, Π can be written in terms of internal energy in the bulk, Π bulk , internal energyin the surface, Π s and work done by external forces, Π ext as, Π = Π bulk + Π s − Π ext (6)where, Π bulk = Z Ω U b d Ω (7) Π s = Z Γ U s d Γ (8) Π ext = Z Γ u u · t d Γ u + Z Ω u · b d Ω − Z Γ φ φq d Γ φ (9)Here, u and φ denote mechanical displacement and electric potential respectively. t is the surface tractionon Γ u , b is the prescribed body force and q is the surface charge density on Γ φ . Γ u and Γ φ are the Neumannboundary for mechanical displacement and electric potential respectively.The weak form of the equilibrium equations can be obtained by finding u ∈ { u = ¯ u on Γ du , u ∈ H (Ω) } and φ ∈ { φ = ¯ φ on Γ dφ , φ ∈ H (Ω) } such that δ Π = 0 = ⇒ Z Ω ε ( δ u ) : C : ε ( u ) d Ω − Z Ω ε ( δ u ) : e · E ( φ ) d Ω − Z Ω E ( δ φ ) · e : ε ( u ) d Ω − Z Ω η ( δ u ) ... µ · E ( φ ) d Ω − Z Ω E ( δ φ ) · µ ... η ( u ) d Ω − Z Ω E ( δ φ ) · κ · E ( φ ) d Ω+ Z Ω η ( δu ) ... g ... η ( u ) d Ω + Z Γ ε s ( δ u ) : τ s d Γ+ Z Γ E s ( δ φ ) · ω s d Γ + Z Γ ε s ( δ u ) : C s : ε s ( u ) d Γ − Z Γ ε s ( δ u ) : e s · E s ( φ ) d Γ − Z Γ E s ( δ φ ) · e s : ε s ( u ) d Γ − Z Γ E s ( δ φ ) · κ s · E s ( φ ) d Γ = Z Γ u δ u · ¯ t d Γ u + Z Ω δ u · b d Ω − Z Γ φ δ φ q d Γ φ (10)for all δ u ∈ { δ u = 0 on Γ du , δ u ∈ H (Ω) } and δ φ ∈ { δ φ = 0 on Γ dφ , δ φ ∈ H (Ω) } . Γ du and Γ dφ are theDirichlet boundary for mechanical displacement and electric potential respectively.
3. Mesh free formulation for flexoelectricity
The numerical discretization of the governing partial differential equations of flexoelectricity requires C continuous basis functions for a Galerkin method. In the present work, we utilize a meshfree methodwith compactly supported radial basis function (CSRBF) shape functions. Popular radial basis functions[30] include the Multi-Quadrics, Gaussian and Thin Plate Splines. These radial basis functions are globallysupported and their accuracy highly depends on the condition number of the collocation matrix. However,4he collocation matrix will be a sparse matrix, well conditioned and compactly supported if we adopt CSRBFshape functions. The Wendland type CSRBFs with C and C continuity proposed in [31] are, f ( r ( x, y )) = max (cid:8) , (1 − r ) (cid:9) (4 r + 1) ∈ C (11) f ( r ( x, y )) = max (cid:8) , (1 − r ) (cid:9) (35 r + 18 r + 3) ∈ C (12)where r ( x, y ) is given by, r i ( x, y ) = d i R = p ( x − x i ) + ( y − y i ) R (13)where d i is the distance of a point of interest ( x, y ) from a knot at ( x i , y i ) . The dimension of support domain, R is given by R = αd c , where α is the shape parameter and d c is the average nodal spacing. It is to be notedthat the value of r i lies between 0 and 1.An approximation for a general function can be written as u h ( x ) = f T ( x ) a + p T ( x ) b (14)where f ( x ) and a denote the vector of CSRBF and expansion coefficients respectively, f T ( x ) = [ f ( x ) , f ( x ) , ...f n ( x )] (15) a T = [ a , a , ...a n ] . (16)In Equations 15 and 16, the variable n stands for the number of nodes in the support domain of the point ofinterest. Here, p ( x ) and b are the vector of polynomial basis functions and coefficients respectively, p T ( x ) = [ p ( x ) , p ( x ) , ...p m ( x )] (17) b T = [ b , b , ...b m ] . (18)In Equations (17) and (18), the variable m stands for the number of terms of polynomial basis. The coefficientvectors a and b can be obtained by solving the following algebraic equations " A P m P Tm ab ) = ( U ) , (19)where U is a vector of nodal values of function u h ( x ) and matrices A and P m are, A = f ( x ) · · · f n ( x ) ... . . . ... f ( x n ) · · · f n ( x n ) (20) P m = p ( x ) · · · p m ( x ) ... . . . ... p ( x n ) · · · p m ( x n ) . (21)5rom Equation (14), interpolation of the nodal function values, U , at any point of interest, x can be writtenas, u h ( x ) = [ f T ( x ) S a + p T ( x ) S b ] U = N ( x ) U (22)As a result, the meshfree CSRBF based shape function, N ( x ) is given by, N ( x ) = f T ( x ) S a + p T ( x ) S b (23)where S a = A − [1 − P m S b ] and S b = [ P Tm A − P m ] − P Tm A − .The polynomial basis functions of linear order are added to the radial basis functions in order to ensurethat the shape functions possess C consistency. The vector p ( x ) given in Equation (17) can be rewrittensuch that m = 3 as, p T ( x ) = [1 x y ] (24)The discrete form of the weak formulation in Equation (10) using the meshfree shape functions is as follows, δu T (cid:18)Z Ω B T u CB u d Ω (cid:19) u + δu T (cid:18)Z Ω B T u e T B φ d Ω (cid:19) φ + δφ T (cid:18)Z Ω B T φ eB u d Ω (cid:19) u + δu T (cid:18)Z Ω H T u µ T B φ d Ω (cid:19) u + δφ T (cid:18)Z Ω B T φ µH u d Ω (cid:19) φ − δφ T (cid:18)Z Ω B T φ κB φ d Ω (cid:19) φ + δu T (cid:18)Z Ω H T u gH u d Ω (cid:19) u + δu T (cid:18)Z Γ B T u M T p τ s d Γ (cid:19) − δφ T (cid:18)Z Γ P T ω s d Γ (cid:19) + δu T (cid:18)Z Γ B T u M T p C s M p B u d Γ (cid:19) u − δu T (cid:18)Z Γ B T u M T p e s T P B φ d Γ (cid:19) φ + δφ T (cid:18)Z Γ B T φ P T e s M p B u d Γ (cid:19) u − δφ T (cid:18)Z Γ B T φ P T κ s P B φ d Γ (cid:19) φ = δu Z Γ u N T ¯ t d Γ u + δu Z Ω N T b d Ω − δφ Z Γ φ N T qd Γ φ (25)where C , e , µ , κ and g are the matrix form of the tensors C ijkl , e ijk , µ ijkl , κ ij and g ijklmn respectively and C s , e s , µ s and κ s are the matrix form of the tensors C sijkl , e sijk , µ sijkl , and κ sij respectively. The gradientand Hessian matrices in Equation 25 are defined as follows, B u = N I,x N I,y N I,y N I,x (26) B φ = − " N I,x N I,y (27)6 u = N I,xx N I,yx N I,yx N I,xx N I,xy N I,yy N I,yy N I,xy (28)where I = 1 , ....n , n is the number of nodes in the support domain of the point of interest and this numbercan be different for different points of interest. The projection matrix is denoted as M P M P = P P P P P P P P P P P P P + P P (29)where the entries of M P are from P , the tangential projection tensor given by I − n ⊗ n . Here, I refers toidentity matrix of rank 2 and n is the outward normal vector to the surface, Γ . The final system of equationscan be written as follows, " K uu + K suu K uφ + K suφ K φu + K sφu K φφ + K sφφ uφ = " F u + F su F φ + F sφ (30)where K uu = Z Ω B T u C B u d Ω + Z Ω H T u gH u d Ω K uφ = Z Ω B T u e T B φ d Ω + Z Ω H T u µ T B φ d Ω = K T φu K φφ = − Z Ω B T φ κB φ d Ω K suu = Z Γ B T u M T p C s M p B u d Γ K sφu = Z Γ B T φ P T e s M p B u d Γ = K suφ T K sφφ = − Z Γ B T φ P T κ s P B φ d Γ F su = Z Γ B T u M T p τ s d Γ F sφ = Z Γ B T φ P T ω s d Γ F u = δu Z Γ u N T ¯ t d Γ u + δu Z Ω N T b d Ω F φ = δφ Z Γ φ N T q d Γ φ (31)7 = C C C C
00 0 C µ = " µ µ µ µ µ µ e = " e e e κ = " κ κ (32) g = l C C C C C C C
00 0 0 C C
00 0 0 0 0 C (33)In Equation 33, the term l is the length scale representing the size dependency of strain gradient effects.
4. Mesh free formulation for flexoelectricity including geometric nonlinearity
In this section, the proposed meshfree formulation is extended to handle geometric nonlinearities in flex-oelectric structures considering surface elasticity. Saint Venant-Kirchhoff material model is considered forflexoelectric solids, the internal energy density given in Equation 1 is modified as, U = 12 S : G + 12 ˜ S ... ˜ G − D · E + 12 S s : G s (34)The total potential energy, Π is given by, Π = Z Ω U d Ω − Z Γ u u · t d Γ u − Z Ω u · b d Ω+ Z Γ φ φq d Γ φ (35)where S is the second Piola-Kirchhoff tensor, ˜ S is the double stress tensor and D is the electric displacementvector; S s is the second Piola-Kirchhoff surface stress tensor; G is the Green Lagrange strain tensor and ˜ G is the gradient of the Green Lagrange strain tensor; E is the electric field vector and G s is the GreenLagrange surface strain tensor. u and φ are mechanical displacement and electric potential respectively. t isthe surface traction on Γ u , b is the prescribed body force and q is the surface charge density on Γ φ . Γ u and Γ φ are the Neumann boundary for mechanical displacement and electric potential respectively.The constitutive equations of the assumed Saint-Venant Kirchhoff material model are as follows S = C : G − e · E ˜ S = − µ · E + g ... ˜ GD = e : G + µ ... ˜ G + κ · E (36)8 aking the first variation of the total potential energy in Equation 35 yields, δ Π = Z Ω S : δG d Ω + Z Ω ˜ S ... δ ˜ G d Ω − Z Ω D · δE d Ω + Z Γ S s : δG s d Γ − Z Γ u δu · t d Γ u − Z Ω δu · b d Ω + Z Γ φ δ φ q d Γ φ = 0 (37)Each term in Equation 37 has to be linearized. The final expression obtained after linearizing each termin Equation 37 are subsequently presented. The intermediate steps are detailed in Appendix A. A totalLagrangian formulation is presented such that all the integrals are performed on the undeformed configurationand derivatives are with respect to the material coordinates.The linearization of the term R Ω S : δG d Ω in Equation 37 can be obtained as L (cid:20)Z Ω S : δG d Ω (cid:21) = Z Ω ¯ S : δ ¯ G d Ω + Z Ω ∆( S : δG ) d Ω (38) Z Ω ∆( S : δG ) d Ω = Z Ω S : ∆( δG ) d Ω + Z Ω δG : ∆ S d Ω= Z Ω S : ∆ ( δG ) d Ω + Z Ω δG : C : ∆ G d Ω − Z Ω δG : e · ∆ E d Ω= Z Ω S : [( ∇ δu ) T ( ∇ (∆ u ))] d Ω + Z Ω δG : C : ∆ G d Ω − Z Ω δG : e · ∆ E d Ω . (39)The linearization of the term R Ω ˜ S ... δ ˜ G d Ω in Equation 37 can be derived as follows, L (cid:20)Z Ω ˜ S ... δ ˜ G d Ω (cid:21) = Z Ω ¯˜ S ... δ ¯˜ G d Ω + Z Ω ∆( ˜ S ... δ ˜ G ) d Ω (40) Z Ω ∆( ˜ S ... δ ˜ G ) d Ω = Z Ω ˜ S ... ∆( δ ˜ G ) d Ω + Z Ω δ ˜ G ... ∆ ˜ S d Ω= Z Ω ˜ S ... ∆( δ ˜ G ) d Ω − Z Ω δ ˜ G ... µ · ∆ E d Ω= Z Ω ˜ S ... [( ∇ δu )( ∇ ∆ u ) + ( ∇ ∆ u )( ∇ δu )] d Ω − Z Ω δ ˜ G ... µ · ∆ E d Ω (41)The linearization of the term R Ω D · δE d Ω in Equation 37 is as follows, L (cid:20)Z Ω D · δE d Ω (cid:21) = Z Ω ¯ D · δ ¯ E d Ω + Z Ω ∆( D · δE ) d Ω (42) Z Ω ∆( D · δE ) d Ω = Z Ω D · ∆ δE d Ω + Z Ω ∆ D · δE d Ω= Z Ω δ E · µ ... ∆ ˜ G d Ω + Z Ω δE · e : ∆ G d Ω + Z Ω δE · κ · ∆ E d Ω (43)The linearization of the term R Γ S s : δG s d Γ in Equation 37 is as follows, L (cid:20)Z Γ S s : δG s d Γ (cid:21) = Z Γ ¯ S s : δ ¯ G s d Γ + Z Γ ∆( S s : δG s ) d Γ (44)9 Γ ∆( S s : δG s ) d Γ = Z Γ S s : ∆( δG s ) d Γ + Z Ω δG s : ∆ S s d Γ= Z Γ S s : ∆ ( δG s ) d Γ + Z Γ δG s : C s : ∆ G s d Γ= Z Ω S s : P · [( ∇ δu ) T ( ∇ (∆ u ))] · P d Γ + Z Γ δG s : C : ∆ G s d Γ (45)The algebraic forms of Equations 38,40,42 and 44 are as follows, L (cid:20)Z Ω S : δG d Ω (cid:21) = δu (cid:18)Z Ω B T ˆ R d Ω (cid:19) + δu (cid:18)Z Ω B T C B d Ω (cid:19) ∆ u + δu (cid:18)Z Ω B T e B φ d Ω (cid:19) ∆ φ + δu (cid:18)Z Ω H T RH d Ω (cid:19) ∆ u (46) L (cid:20)Z Ω ˜ S ... δ ˜ G d Ω (cid:21) = δu (cid:18)Z Ω H TD ˆ R D d Ω (cid:19) + δu (cid:18)Z Ω H TD µ T B φ d Ω (cid:19) ∆ φ + δu (cid:18)Z Ω H T R TD H d Ω (cid:19) ∆ u + δu (cid:18)Z Ω H T R D H d Ω (cid:19) ∆ u (47) L (cid:20)Z Ω D · δE d Ω (cid:21) = − δφ Z Ω B T φ ˆ D d Ω − δφ (cid:18)Z Ω B T φ µH u d Ω (cid:19) ∆ u − δφ (cid:18)Z Ω B T φ eB d Ω (cid:19) ∆ u + δφ (cid:18)Z Ω B Tφ κB φ d Ω (cid:19) ∆ φ (48) L (cid:20)Z Γ S s : δG s d Γ (cid:21) = δu (cid:18)Z Γ B T M T p ˆ R s d Γ (cid:19) + δu (cid:18)Z Γ B T M T p C s M p B d Ω (cid:19) ∆ uδu (cid:18)Z Γ H T P Tn R s P n H d Γ (cid:19) ∆ u (49)All the matrices involved in Equations 46, 47, 48, 49 are presented in Appendix B. The final algebraic formof linearization of Equation 37 can be written as, K ∆ U = F ext − F int (50)where, K = " K uu K uφ K φu K φφ ∆ U = " ∆ u ∆ φ (51)10 able 1: Electromechanical properties of STO Elastic Constants Dielectric constants Flexoelectric constants [32] C =310 GP a κ =2.66 C/ ( GV − m ) µ =-0.26 nC/mC =115 GP a κ =2.66 C/ ( GV − m ) µ =-3.74 nC/mC =310 GP a µ =-3.56 nC/mC =54 GP a K uu = Z Ω B T CB d Ω + Z Ω H T RH d Ω+ Z Ω H T R TD H d Ω + Z Ω H T R D H d Ω+ Z Γ B T M pT C s M p B d Γ + Z Γ H T P n R s P n H d Γ (52) K uφ = Z Ω B T e T B φ d Ω + Z Ω H TD µ T B φ d Ω = K T φu (53) K φφ = − Z Ω B Tφ κB φ d Ω (54) F int = " F u F φ F u = Z Ω B T ˆ R d Ω + Z Ω H TD ˆ R D d Ω + Z Γ B T M T p ˆ R s d Γ F φ = Z Ω B T φ ˆ D d Ω (55)The nonlinear equation 50 is solved by using the Newton-Raphson iterative scheme. Solving this equation,gives the deflection and voltage responses of flexoelectric structures that undergo large deformations.
5. Numerical Examples
In this section, the proposed meshfree formulation is utilised to analyse flexoelectric cantilever beamsaccounting for surface effects and also to study the influence of geometric nonlinearity on their voltageoutput. As an initial step, the meshfree formulation is validated by determining the energy conversion factor(ECF) of a cantilever beam. The ECF is given by the ratio between stored electrical energy and mechanicalenergy in the flexoelectric structure.
A cantilever beam subjected to a mechanical point load at the free end is analysed in order to validatethe proposed meshfree formulation. The cantilever beam has an aspect ratio of 6. The Young’s modulus,Y is assumed to be 100 GPa. For validation, only the flexoelectric constant µ and dielectric constant κ are considered non-zero and they are assumed to be 10 nC/m and 1 nC/V m respectively. The beam isdiscretized by 121 ×
21 nodes with uniform spacing and the background mesh for numerical integration is120 ×
20. The variation of the energy conversion factor (ECF) with decreasing depth of the 1D flexoelectricbeam is shown in Figure 1. The figure shows good agreement between the numerical and analytical k values.11 epth of beam, d (nm)
20 40 60 80 100 E C F AnalyticalMeshfree
Figure 1: The variation of ECF, k with beam depth for a one dimensional beam model. The analytical energy conversion factor, k for a 1D model is given in [33] as, ECF anl = k = r κY (cid:16) µd (cid:17) (56) In order to further validate the proposed formulation, a flexoelectric tube made of STO is analysed withplane strain assumption. The tube with an inner radius, r i of 10 µm and outer radius, r o of 20 µm issubjected to a radial displacement of 0.045 µm and 0.05 µm at r i and r o respectively [23, 24]. The tube isgrounded along the inner face and a voltage of 1 V is applied along the outer face. The nodal distribution ofthe quarter model is as shown in Figure 2(a). The distribution of electric potential obtained for a length scale, l of 2 µm is shown in figure 2(b). The variation of electric potential along the thickness of the tube is shownin figure 2(c). The analytical results for the flexoelectric tube model with assumed material parameters isderived by Mao et al . [23]. The results presented in Figure 2(c) shows good agreement between the analyticaland numerical results. A cantilever beam made of cubic STO is analysed in this section. The material properties of STO isgiven in Table 1. The bottom and top face of the beam are coated with electrodes. The bottom electrode isgrounded and the top electrode is free to have a potential value.The dimension of the beam is 1200 ×
100 nm. The beam is subjected to a point load of 10 nN at the freeend. The length scale, l of g tensor is taken as 0 (i.e.) g is not considered in this analysis. The beam hasan almost linear variation of potential along the beam depth as shown in Figure 3. The potential obtainedat the top face is 30 mV. Now if we fix the aspect ratio to be 12, and reduce the beam depth for instance to40 nm, then the potential at the top face is 74 mV. If the tensor, g is included in the analysis with a length12 , nm × Y , n m × (a) X, nm × Y , n m × -0.200.20.40.60.81 (b) Radius r, in µ m
10 12 14 16 18 20 E l e c t r i c po t en t i a l , φ i n V -0.4-0.200.20.40.60.81 AnalyticalNumerical (c)Figure 2: (a) Nodal distribution for the quarter tube model, (b) Electric potential distribution across the tube cross section, (c)The variation in output voltage along the radius of the tube. X Y Figure 3: The potential distribution in a flexoelectric beam made of STO. epth of beam, nm
40 50 60 70 80 90 100 E C F × -4 g=0g 0 Figure 4: The variation of ECF with depth of beam including and excluding strain gradient tensor, g . The length scale is takenas 5 nm. scale, l of 5 nm, then the potential obtained at the top face of the 40 nm beam depth reduces to 52 mV.The change in energy conversion factor with depth of beam for an aspect ratio of 12, excluding and including g tensor is shown in Figure 4. The energy conversion factor for 40 nm beam depth with and without including g tens or are 1.1e-4 and 1.61e-4 respectively. The inclusion of strain gradient elasticity increases the stiffnessof the beam, reduces the voltage obtained and as a result reduces the ECFs. Note that the difference betweenthe energy conversion factors with and without the g tensor increases with reduction in depth of the beam.On the other hand, the influence of the length scale, l on the output voltage is shown in Figure 5. It can beseen that, when the length scale is increased from 0 to 5 nm, the output voltage of the 100 nm beam reducesfrom 30 mV to 27.5 mV (6.3 % reduction). While for the same increase in length scale, the output voltage ofthe 40 nm beam reduces from 73.6 mV to 51.7 mV (28.6 % reduction). In this section, we analyse a Zinc Oxide nano cantilever beam. The interplay between piezoelectric,surface elastic, surface piezoelectric and flexoelectric effects is studied. Zinc Oxide is the ideal material forperforming this study because it is widely used in several nanoscale energy harvesters [34, 3] and studies onsurface properties of Zinc Oxide [35] is available.The cantilever ZnO beam is of length, 120 nm and width, 15 nm. A point load of 10 nN is applied inx-direction at the mid-point of the top face. The beam is fixed at the bottom and free at the top. The beamis poled along the length (y-direction). The bottom end of the beam is grounded. The elastic, piezoelectric,surface elastic and surface piezoelectric properties of ZnO are given in Table 2 and Table 3. The residualstress, τ s and residual electric displacement, ω s are not considered in the study. The flexoelectric constantof ZnO, µ , µ and µ are assumed to be 2 nC/m, 2 nC/m and 0.5 nC/m respectively.14 eam depth, nm
40 50 60 70 80 90 100 V o l t age , m V l =0l =1 nml =2.5 nml =5 nm Figure 5: Output voltage for varying beam depths for internal length scales of 0,1,2.5 and 5 nm.Table 2: Material properties of bulk ZnO
Elastic Constants Piezoelectric constants Dielectric constants C =206 GP a e =-0.58 C/m κ =0.0811 C/ ( GV − m ) C =117 GP a e =1.55 C/m κ =0.112 C/ ( GV − m ) C =211 GP a e =-0.48 C/m C =44.3 GP a
Table 3: Material properties of ZnO surface
Elastic Constants Piezoelectric constants C s =44.2 N/m e s =-0.216 nC/mC s =14.2 N/m e s =0.451 nC/mC s =35 N/m e s =-0.253 nC/mC s =11.7 N/m a) (b)Figure 6: The potential distribution in the ZnO beam considering (a) Flexoelectricity, piezoelectricity and surface effects, (b)Pure Flexoelectricity. The potential distribution across the beam width is shown in Figure 6. The combination of bulk piezoelec-tricity and surface elastic effect results in a potential of +1.18 V to -1.18 V at the top face. The combinationof bulk piezoelectricity and surface piezoelectric effect leads to a potential varying from +1.5 V to -1.5 Vat the top face of the beam. Finally, the combination of bulk piezoelectricity, bulk flexoelectricity, surfaceelasticity and surface piezoelectricity results in a potential of +1.7 V to -1.7 V at the top face as shown inFigure 6(a). The variation of electric potential considering only flexoelectric effect is +0.3 V to -0.3 V at thetop face (Figure 6(b)). The relative influence of the different phenomenon on the output voltage is shown inFigure 7. The contribution of flexoelectric effect to output voltage is higher compared to the contribution ofsurface effects and the difference in the contributions to output voltage increases as the beam width decreases.For a 40 nm wide beam, the flexoelectric and surface effect contributions are 7 % and 1 % respectively. Whilefor a width of 15 nm, the difference between the contributions is higher, the flexoelectric and surface effectcontributions are 18 % and 7 % respectively.The change in ECF with width of beam is shown in Figure 8. The pattern is similar to the one obtainedfor output voltage. The ECF for pure Piezoelectricity and Piezo + Surface effects for 15 nm wide beamare 0.0067 and 0.00728 respectively. While for the same beam width, the ECF considering Piezo + Surfaceeffects + Flexoelectricity is 0.014. The percentage contribution of flexoelectricity and surface effects to the totalECF are 48 % and 8 % respectively. There is a discrepancy in the percentage contribution of flexoelectricityto ECF and output voltage. This is because the potential due to flexoelectricity reaches its peak near thefixed end and reduces significantly along the length. So, though the flexoelectric contribution to total ECFis 48 % , the contribution of flexoelectricity to total voltage measured at the top face (y=120 nm) is only 18 % . In this section, the flexoelectric response in the nonlinear regime is studied. The nonlinearity emerges dueto large deformation of the flexoelectric cantilever beam. The flexoelectric beam is assumed to be made ofSTO, in addition to flexoelectricity, the surface elasticity of STO is also considered. The material properties16 idth of beam, nm
20 40 60 80 100 V o l t age , V e 0,e s =0,C s =0, µ =0e 0,e s s µ =0e 0,e s s µ Figure 7: The variation of output voltage for ZnO beam with width for three cases, Piezoelectricity, Piezoelectricity+surfaceeffects and Piezoelectricity+Surface effects+Flexoelectricity.
Width of beam, nm
20 40 60 80 100 E C F e 0, C s = 0 , e s = 0 , µ = 0e 0, C s s µ = 0e 0 C s s µ Figure 8: The variation of ECF for ZnO beam with width for three cases, Piezoelectricity, Piezoelectricity+Surface effects andPiezoelectricity+Surface effects+Flexoelectricity.
17f STO are given in Table 1. The surface elastic constants of STO are not available in literature and areassumed to be, C s = C s = 310 GP a , C s = 115 GP a and C s = 54 GP a . The beam is subjected tomechanical deformation by a point load of 10 nN at the free end. The load is applied in increments of 1 nN.
In each increment, the tangent stiffness matrix is determined and the Newton-Raphson method is adoptedto minimize the residual. Two different beams each of thickness 50 nm and 30 nm respectively are analysed.The aspect ratio of the beams is fixed as 12. The fixed end of the cantilever is grounded.The variation of maximum output voltage with load steps is shown in Figure 9. As shown in Figure 9(a),at the end of ten load steps the voltage due to flexoelectricity and surface elasticity for 50 nm thick beam is7.8 mV. The ratio between nonlinear and linear voltage is 0.9 (i.e.) the nonlinear voltage deviates from thelinear response by . If surface effects are not considered, then the final output voltage is 8.5 mV. In caseof 30 nm thick beam, as shown in Figure 9(b), the final output flexoelectric voltage is 11 mV and 12 mVconsidering and ignoring surface effects respectively.The variation of free end deflection with load steps is shown in Figure 10. The free end deflection of 50nm thick beam after ten load steps is 68 nm (Figure 10(a)). The ratio between nonlinear and linear deflectionfor a load value of 10nN for 50 nm thick beam is 0.89 (i.e.) a deviation of nonlinear displacement is from linear displacement. In the absence of surface effects, the final output deflection of 50 nm thick beamis 75 nm. For the case of 30 nm thick beam, as shown in Figure 10(b), the final deflection is 56 nm and63.5 nm considering and ignoring surface effects respectively.From Figures 9 and 10, it is clear that the beambecomes stiffer in the presence of surface elasticity, which in turn leads to reduction in output voltage. Incase of 50 nm thick beam, due to surface elasticity the absolute value of final voltage reduces from 8.5 mV to7.8 mV. While for 30 nm thick beam, the voltage reduction is from 12 mV to 11 mV. The variation of energyconversion factor with load steps is shown in Figure 13. The rate of increase in ECF with load steps is higherfor a 30 nm beam compared to the 50 nm beam (Aspect ratio = 12). It is to be noted that, in the absenceof geometric nonlinearity, the energy conversion factor is a constant value and remains independent of theapplied loads, whereas the consideration of geometric nonlinearity has led to change in ECF value with loadincrements.The contour plot showing the normal strain in x-direction is given in Figures 11 and 12. Figure 11 showsthat at the first load step, the nonlinear strain terms are negligible and the Green-Lagrange strain contourand linear strain contour are quite similar. While the strain contour at the 10 th load step given in Figure 12shows that nonlinear terms play significant role in determining the Green-Lagrange strain. Consequently,the linear strain contour and Green Lagrange strain contours become dissimilar. The studies performed inthis section show that the surface elasticity can reduce the output flexoelectric voltage. For instance, it isobserved that in case of a 30 nm thick beam the reduction in the final voltage is 1 mV due to surface elasticeffects (Figure 9). The surface elastic effects can increase with increase in the surface area of the beam.Therefore it is worthwhile to study a tapered beam model with inclined top and bottom faces and comparetheir response with the rectangular beams studied previously.Considering an average thickness of 30 nm, a tapered cantilever beam, T B , of length 360 nm and ofdepth, d l =40 nm at x=0, d r =20 nm at x=360 nm is subjected to a load increment of 1 nN over ten loadsteps. At the tenth load step, the output flexoelectric voltage in T B reduces from 15 mV to 13 mV due tosurface elastic effects as shown in Figure 15(a). Secondly, a tapered cantilever beam T B of length 360 nmand of depth, d l =45 nm at x=0, d r =15 nm at x=360 nm is considered. The final voltage reduces from 12mV to 10.4 mV as shown in Figure 15(b). In both the tapered beams T B and T B with average thicknessof 30 nm, the voltage reduces by . due to surface elastic effects. This is higher than the reduction of18 oad steps E l e c t r i c po t en t i a l , φ i n V -0.01-0.008-0.006-0.004-0.0020 nonlinear ( µ
0, C s µ
0, C s =0) (a) Load steps E l e c t r i c po t en t i a l , φ i n V -0.012-0.01-0.008-0.006-0.004-0.0020 nonlinear ( µ
0, C s µ
0, C s =0) (b)Figure 9: The variation in voltage with load increments for flexoelectric beam made of STO of depth (a) 50 nm and (b) 30 nm(Aspect ratio=12). Load steps D e f l e c t i on , n m nonlinear( µ
0, C s µ
0, C s =0) (a) Load steps D e f l e c t i on , n m nonlinear( µ
0, C s µ
0, C s =0) (b)Figure 10: The load deflection curve for flexoelectric beam made of STO of depth (a) 50 nm and (b) 30 nm (Aspect ratio=12). . (from 11.8 mV to 10.8 mV) in 30 nm thick rectangular beam. Therefore, as expected the influence ofsurface effects increases with increase in surface area of the beam and the negative influence on flexoelectricvoltage is higher for a tapered beam compared to a rectangular beam of same volume.In summary, the results obtained in this section show that considering nonlinear terms in strain andgradient of strain, is more essential as the influence of flexoelectricity on output voltage gets significant innanoscale. We conclude that the geometric nonlinearity cannot be ignored if one analyses flexoelectric beamsof dimensions of under 100 nanometers when subjected to loads in the range of 10 nNs. It is to be notedthat the flexoelectric material, STO used in this example has a lesser flexoelectric constant of only 1.4 V,while flexoelectric constant of a dielectric material can even be upto 10 V based on the theoretical upperlimit estimated by Kogan et al. [8]. 19 a)(b)(c)Figure 11: The strain contour ( G , ε , η ) for flexoelectric beam made of STO of depth 50 nm at load step 1 (a) GreenLagrange strain (b) Linear strain (c) Nonlinear part of strain.(a)(b)(c)Figure 12: The strain contour ( G , ε , η ) for flexoelectric beam made of STO of depth 50 nm at load step 10 (a) GreenLagrange strain (b) Linear strain (c) Nonlinear part of strain. oad steps E C F × -6 d = 50 nmd = 30 nm Figure 13: The variation in ECF with load increments for flexoelectric beam made of STO of depth 50 nm and 30 nm (Aspectratio =12). d l d r l = 600 nm Figure 14: Tapered flexoelectric beam model.
Load steps E l e c t r i c po t en t i a l , φ i n V -0.015-0.01-0.0050 nonlinear( µ
0, C s µ
0, C s =0) (a) Load steps E l e c t r i c po t en t i a l , φ i n V -0.012-0.01-0.008-0.006-0.004-0.0020 nonlinear( µ
0, C s µ
0, C s =0) (b)Figure 15: The variation in voltage with load increments for flexoelectric tapered beam made of STO of type (a) T B and (b) T B . . Conclusion A CSRBF based meshfree formulation is presented in this paper to handle geometric nonlinearity in flex-oelectric structures. In addition to flexoelectricity, the surface effects are also considered in the analysis ofnano-sized two dimensional structures. Flexoelectric beams made of cubic STO, which is non-piezoelectricand ZnO, which is piezoelectric are analysed. The meshfree analysis shows that for ZnO, the contribution ofsurface effects to the output voltage of a nanosized cantilever structure (of width 15 nm) is smaller comparedto the contribution of flexoelectricity.The analysis of flexoelectric nanostructures undergoing large deformation shows that the difference betweennonlinear and linear flexoelectric voltage increases with reduction in beam depth. The surface elastic effectsstiffen the beam leading to reduction in output flexoelectric voltage. The influence of surface elasticity ishigher in tapered beams compared to rectangular beams. In future, the presented formulation will be ex-tended to study nonlinear flexoelectricity under dynamic excitations and the influence of nonlinearity onresponse bandwidth of nanosized flexoelectric energy harvesters will be investigated.
Appendix A. Intermediate Steps in the derivation of nonlinear meshfree formulation for flex-oelectricity
The terms G , δG , ∆ δG and ∆ S in Equation 38 are as follows, G = 12 ( u i,j + u j,i + u k,i u k,j ) (A.1) δG = 12 ( δu i,j + δu j,i + δu k,i u k,j + u k,i δu k,j ) (A.2) ∆ δG = 12 ( δu k,i ∆ u k,j + ∆ u k,i δu k,j )= δu k,i ∆ u k,j (A.3) S = C : G − e · E ∆ S = C : ∆ G − e · ∆ E (A.4)The terms ˜ G , δ ˜ G , ∆ δ ˜ G and ˜ S in Equation 40 are as follows, ˜ G = 12 ( u i,jk + u j,ik + u k,ij u k,j + u k,i u k,ji ) δ ˜ G = 12 ( δu i,jk + δu j,ik + δu k,ij u k,j + u k,ij δu k,j + δu k,i u k,ji + u k,i δu k,ji ) ∆ δ ˜ G = 12 ( δu k,ij ∆ u k,j + ∆ u k,ij δu k,j + δu k,i ∆ u k,ji + ∆ u k,i δu k,ji )= δu k,ij ∆ u k,j + ∆ u k,ij δu k,j ˜ S = − µ · E (A.5)22he terms ∆ D , δE and ∆ δE in Equation 42 are as follows, D = e : G + µ ... ˜ G + κ · E ∆ D = e : ∆ G + µ ... ∆ ˜ G + κ · ∆ E E i = − φ ,i δE i = − δφ ,i ∆ δE i = 0 (A.6)The terms G s , δG s , ∆ δG s and ∆ S s in Equation 44 are as follows, G s = 12 P ij ( u i,j + u j,i + u k,i u k,j ) P ji (A.7) δG s = 12 P ij ( δu i,j + δu j,i + δu k,i u k,j + u k,i δu k,j ) P ji (A.8) ∆ δG s = 12 P ji ( δu k,i ∆ u k,j + ∆ u k,i δu k,j ) P ji = P ji δu k,i ∆ u k,j P ji (A.9) S s = C s : G s ∆ S s = C s : ∆ G s (A.10) Appendix B. Matrices - nonlinear meshfree formulation for flexoelectricity
The expressions defining the matrices B , B φ , H , H , H u , H D , ˆ R , R , ˆ R D , R D , ˆ D , R s , ˆ R s and P n are as follows, B = N I,x N I,y N I,y N I,x + A H (B.1) A = ∂u Ix ∂x ∂u Iy ∂x ∂u Ix ∂y ∂u Iy ∂y∂u Ix ∂y ∂u Ix ∂x ∂u Iy ∂y ∂u Iy ∂x (B.2) H = N I,x N I,y N I,x N I,y (B.3) B φ = " N I,x N I,y (B.4) H D = H u + A D H (B.5)23 D = ∂u Ix ∂x ∂u Iy ∂x ∂u Ix ∂x ∂u Iy ∂x ∂u Ix ∂y ∂u Iy ∂y
00 0 0 ∂u Ix ∂y ∂u Iy ∂y∂u Ix ∂y ∂u Ix ∂x ∂u Iy ∂y ∂u Iy ∂x ∂u Ix ∂y ∂u Ix ∂x ∂u Iy ∂y ∂u Iy ∂x (B.6) H = N I,xx N I,xy N I,yx N I,yy N I,xx N I,xy N I,yx N I,yy (B.7) R = S S S S S S S S (B.8) ˆ R = h S S S i T (B.9) R D = S S S S S S S S S S S S S S S S (B.10) ˆ R D = h S S S S S S i T (B.11) ˆ D = " D D (B.12) ˆ R s = h τ s + S s τ s + S s τ s + S s i T (B.13) R s = τ s + S s τ s + S s τ s + S s τ s + S s τ s + S s τ s + S s τ s + S s τ s + S s (B.14)24 n = P P P P P P P P (B.15) Acknowledgments
X. Zhuang acknowledges the support from State Key Laboratory of Structural Analysis for IndustrialEquipment (GZ1607) and National Science Foundation of China (11772234).
References [1] R. Mbarki, J. Haskins, A. Kinaci, T. Cagin, Temperature dependence of flexoelectricity in
BaT iO and SrT iO perovskite nanostructures, Physics Letters A 378 (2014) 2181–2183.[2] S. Priya, H.-C. Song, Y. Zhou, R. Varghese, A. Chopra, S.-G. Kim, I. Kanno, L. Wu, D. S. Ha, J. Ryu,et al., A review on piezoelectric energy harvesting: materials, methods, and circuits, Energy Harvestingand Systems 4 (2017) 3–39.[3] Z. L. Wang, J. Song, Piezoelectric nanogenerators based on zinc oxide nanowire arrays, Science 312(2006) 242–246.[4] N. Sharma, R. Maranganti, P. Sharma, On the possibility of piezoelectric nanocomposites without usingpiezoelectric materials, Journal of the Mechanics and Physics of Solids 55 (2007) 2328–2350.[5] U. K. Bhaskar, N. Banerjee, A. Abdollahi, Z. Wang, D. G. Schlom, G. Rijnders, G. Catalan, A flexo-electric microelectromechanical system on silicon, Nature nanotechnology 11 (2016) 263–266.[6] V. Mashkevich, Electrical, optical, and elastic properties of diamond-type crystals ii. lattice vibrationswith calculation of atomic dipole moments, Soviet Physics Jetp 5 (1957).[7] A. K. Tagantsev, Piezoelectricity and flexoelectricity in crystalline dielectrics, Physical Review B 34(1986) 5883.[8] S. M. Kogan, Piezoelectric effect under an inhomogeneous strain and acoustic scattering of carriers incrystals, Fiz. Tverd Tela 5 (1963) 2829–2831.[9] W. Ma, L. E. Cross, Large flexoelectric polarization in ceramic Lead Magnesium Niobate, AppliedPhysics Letters 79 (2001) 4420–4422.[10] W. Ma, L. E. Cross, Flexoelectric polarization of Barium Strontium Titanate in the paraelectric state,Applied Physics Letters 81 (2002) 3440–3442.[11] W. Ma, L. E. Cross, Flexoelectricity of Barium Titanate, Applied Physics Letters 88 (2006) 232902.[12] R. Maranganti, N. D. Sharma, P. Sharma, Electromechanical coupling in nonpiezoelectric materials dueto nanoscale nonlocal size effects: Green’s function solutions and embedded inclusions, Physical ReviewB 74 (2006) 014110. 25
13] P. Zubko, G. Catalan, A. Buckley, P. Welche, J. Scott, Strain-gradient-induced polarization in
SrT iO single crystals, Physical Review Letters 99 (2007) 167601.[14] Z. Yan, L. Y. Jiang, Flexoelectric effect on the electroelastic responses of bending piezoelectricnanobeams, J. Appl. Phys. 113 (2013) 194102.[15] C. C. Liu, S. L. Hu, S. P. Shen, Effect of flexoelectricity on electrostatic potential in a bent piezoelectricnanowire, Smart Mater. Struct. 21 (2012) 115024.[16] S. Zhang, H. Yao, W. Fan, Y. Hao, X. Wu, D. Hou, Effects of flexoelectricity and surface elasticityon piezoelectric potential in a bent ZnO nanowire, in: IOP Conference Series: Materials Science andEngineering, volume 167, IOP Publishing, p. 012023.[17] H. Chen, A. Soh, Y. Ni, Phase field modeling of flexoelectric effects in ferroelectric epitaxial thin films,Acta Mechanica 225 (2014) 1323–1333.[18] H. Chen, S. Zhang, A. Soh, W. Yin, Phase field modeling of flexoelectricity in solid dielectrics, Journalof Applied Physics 118 (2015) 034106.[19]
A. Abdollahi, C. Peco, D. Millán, M. Arroyo, I. Arias, Computational evaluation of the flexoelectric effect in dielectric solids, Journal of Applied Physics 116 (2014) 093502.[20]
A. Abdollahi, D. Millán, C. Peco, M. Arroyo, I. Arias, Revisiting pyramid compression to quantifyflexoelectricity: A three-dimensional simulation study, Physical Review B 91 (2015) 104103.[21] H. Ghasemi, H. S. Park, T. Rabczuk, A level-set based IGA formulation for topology optimization offlexoelectric materials, Computer Methods in Applied Mechanics and Engineering 313 (2017) 239–258.[22] S. S. Nanthakumar, X. Zhuang, H. S. Park, T. Rabczuk, Topology optimization of flexoelectric struc-tures, Journal of the Mechanics and Physics of Solids 105 (2017) 217 – 234.[23] S. Mao, P. K. Purohit, N. Aravas, Mixed finite-element formulations in piezoelectricity and flexoelec-tricity, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences472 (2016).[24] F. Deng, Q. Deng, W. Yu, S. Shen, Mixed finite elements for flexoelectric solids, Journal of AppliedMechanics 84 (2017) 081004.[25]