A meshfree formulation for large deformation analysis of flexoelectric structures accounting for the surface effects

A meshfree form ulation for large deformation analysis of flexo electric structures accoun ting for the surface effects Xiao ying Zh uang a,b, ∗ , S.S.Nan thakumar b , Timon Rabczuk c, ∗∗ a Col le ge of Civil Engine ering, T ongji University, 1239 Siping R o ad, 200092 Shanghai, China. b Institute of Continnum Me chanics, L eibniz University Hannover, App elstr asse 11A, D-30167 Hannover Germany. c Institute of Structur al Me chanics, Bauhaus University W eimar, Marienstr asse 5, 99423 W eimar, Germany. Abstract In this w ork, we presen t a compactly supp orted radial basis function (CSRBF) based meshfree method to analyse geometrically nonlinear flexo electric nanostructures considering surface effects. Flexo electricit y is the p olarization of dielectric materials due to the gradien t of strain, which is different from piezoelectricity in which polarization is dep endent linearly on strain. The surface effects gain prominence as the size of the structure tends to nanoscale and so their consideration is inevitable when flexo electric nanostructures are analysed. First, the proposed meshfree formulation is v alidated and the influence of nonlinear strain terms on the energy con version ability of flexo electric b eams made of a non-piezoelectric material like cubic Strontium Titanate is studied. Subsequen tly , the meshfree form ulation for nonlinear flexo electricit y is extended to include nonlinear surface effects. It is determined that the surface effects can ha v e notable influence on the output flexoelectric voltage of nano-sized can tilev er structures in the nonlinear regime. Keyw ords: Meshfree metho d; Nolinear Flexo electricit y; Geometric nonlinearit y; Surface effects 1. In tro duction Flexo electricit y is the generation of electric polarization under mechanical strain gradien t or mechanical deformation due to electric field gradien t (con verse flexo). It is a more general phenomenon than the linear c hange in p olarization due to stress, kno wn as the piezo electric effect. Flexo electric p olarization is restricted not only to non-cen trosymmetric crystals ev en tually op ening up possibilities for non toxic electromec hanical materials for biomedical application. Piezo electricit y can be c haracterised b y a third rank tensor and is observed only in non-centrosymmetric crystals (21 t yp es). In contrast, flexo electricity can b e mathematically defined b y a fourth order tensor and can b e observ ed in materials of an y symmetry (32). The reason b ehind is that the homogeneous strain relies on lack of symmetry of materials for p olarization, on the other hand, the strain gradient can break the local cen trosymmetry of materials inducing p olarization. The strain gradient scales with the size of sp ecimen leading to the possibility of significan t flexoelectric effect at the length scale of nanometers. Piezo electricity exists only b elo w Curie temp erature, while flexo electricity being symmetry indep enden t do es not hav e a temp erature constrain t [1]. The high energy conv ersion abilit y of piezo electric materials makes them the prominen t constituent in sev eral micro [2] and nano-sized [3] energy harv esters developed. While, the recen t ∗ Corresponding author ∗∗ Corresponding author Email addr esses: zhuang@ikm.uni-hannover.de (Xiaoying Zh uang), timon.rabczuk@uni-weimar.de (Timon Rab czuk) Pr eprint submitte d to Elsevier Novemb er 13, 2019 researc hes sho w the p ossibilit y of an energy harvester made of non-piezo electric materials exploiting flexo elec- tricit y [4]. Nano electromechanical systems like actuators [5] are fabricated using non-piezo electric materials lik e Stron tium Titanate and are shown to pro duce curv ature/electric field ratio of 3.33 M V − 1 comparable to the ratio of 5.2 M V − 1 in piezo electric Lead Zirconium Titanate bimorph. The flexo electricity also offers the adv antage of c ho osing Lead-free materials as constituen ts in sensors, actuators and energy harvesters. The theory of flexo electricity w as first identified w ay back in the 1960s by Mashkevic h and T olpygo [6], follo wed later b y the w ork of T agan tsev [7] in which bulk and surface mec hanisms that can cause p olarization due to strain gradient were determined. Mean while, Kogan [8] made a theoretical estimate of the flexo electric co efficien t to be of the order of e/a ≈ 10 − 9 C /m , where e is the electronic charge and a is the lattice parameter. The series of exp erimen tal w orks b y Cross and co-work ers [9, 10, 11] sparked interest ov er the potential of flexo electric materials as a substitute to piezo electric materials. These exp erimental studies on ceramics with cubic symmetry lik e Barium Stron tium Titanate (BST) and Barium Titanate (BTO) revealed higher v alues of p eak flexo electric co efficients in the range of 50 µC /m . A tomistically , Maranganti et al . [12] determined the flexoelectric coefficients of sev eral ferro electric and non-ferro electric crystals. Though there is discrepancy b etw een theoretical and exp erimental flexo electric co efficients of Barium Titanate (BTO), the theoretical estimations are of the same order compared to the exp erimen tal results of Zubko et al . [13] for Strontium Titanate (STO) crystals. Several w orks are av ailable in literature that presents analytically deriv ed electro-elastic field for nanobeams and nanowires ha ving flexo electric effect and surface effects [14, 15, 16]. Nevertheless, the analytical solutions are applicable only to simplified one-dimensional mo dels and so n umerical methods to analyse flexoelectric structures are required. Con ven tional finite elemen t metho d (FEM) cannot b e utilised for analysing flexo electric structures as the fourth order partial differential equations go verning flexoelectricity necessitates C 1 con tinuit y of displacement field. Phase field modelling of flexo electricit y in an epitaxial thin film made of Barium Titanate is presented b y Chen et al . [17], follo wed by which analysis of a t w o phase system is p erformed [18]. Meshfree shape functions offer the adv antage of having higher order contin uit y , making them a fav ourable class of numerical metho ds to analyse flexoelectric structures. A numerical approac h to analyse tw o and three dimensional truncated p yramid shap ed structure due to flexo electricit y utilising lo cal maximum entrop y (LME) meshfree metho d is presented by Abdollahi et al . [19, 20]. Ghasemi et al . [21] proposed an IGA formulation exploiting the higher order con tinuit y of NURBS shape functions. In [22, 23, 24], mixed FE formulations are prop osed for analysis of t wo dimensional flexo electric structures. Though the mixed FE formulation requires only C 0 con tinuit y , the num b er of no dal DOF s required is m uc h higher. F or example, in the flexo electric element prop osed in [22], degrees of freedom in the corner no des are tw o displacement DOF s, four displacement gradien t DOF s, one electric p otential DOF and four Lagrange multiplier DOF s. It is to b e noted that in the w ork of Nan thakumar et al. [22] flexo electric nanob eams with surface effects, made of Barium Titanate, are analysed and optimized using the mixed FE formulation. Also there are computational works a v ailable in literature that specifically analyse nanobeams with surface effects [25, 26, 27], in which extended finite element metho d is the n umerical metho d adopted. The nonlinear electro-elasticit y of soft dielectrics combined with flexo electricit y is analysed b y Y vonnet et al . [28], adopting finite ele men t discretization ( C 1 Argyris triangular elemen ts) and consistent linearizations. As a shortcoming the authors hav e stated that due to instability , the utilised form ulation could not simulate the en tire nonlinear range. Motiv ated by all these works on flexo electricity , in the presen t w ork, a compactly supp orted radial basis function (CSRBF) based meshfree formulation is prop osed to analyse flexoelectric b eams sub jected to large deformation considering surface effects. Though there are w orks av ailable in literature on analysing flexo- 2 electric structures using a meshfree metho d [19, 20], to the b est of our knowledge this is the first work on a meshfree formulation to handle nonlinearit y in flexo electric nanostructures accounting for surface effects. The meshfree shape functions with higher order con tinuit y are adv an tageous compared to complex mixed FE [22] formulations, mainly b ecause the meshfree form ulation requires the discretization of ’only’ displacemen t and electric p otential fields. The outline of the pap er is as follows: Section 2 presents the gov erning equations of flexo electricit y . Section 3 describ es the meshfree form ulation for flexo electricit y including surface elasticity and surface piezo- electricit y . Linearization of the weak form and subsequen t meshfree discretization is sho wn in Section 4. Finally , n umerical examples on analysis of tw o-dimensional flexoelectric structures with surface effects con- sidering geometric nonlinearity are presented in Section 5. 2. Go v erning equations of flexo electricit y with surface effects The mathematical mo delling of flexoelectricity is based on the extended linear theory of piezo electricity with additional strain gradien t terms. A general internal energy densit y function, U inv olving strain energy , electrostatic energy and terms including strain gradient is presented by Shen et al. [29]. The in ternal energy densit y function, U is as follows, U = U b + U s U b = 1 2 ε : C : ε − E · e : ε − E · µ . . . η − 1 2 E · κ · E + 1 2 η . . . g . . . η U s = U s 0 + τ s : ε s + ω s · E s + 1 2 ε s : C s : ε s − E s · e s : ε s − 1 2 E s · κ s · E s (1) where U b and U s are the bulk and surface energy density functions resp ectiv ely . U s 0 is the surface free energy densit y . ε is the linear strain tensor, E is the electric field tensor, ε s and E s are their corresp onding surface coun terparts. τ s and ω s are the residual surface stress and residual surface electric displacements resp ectiv ely . η is the strain gradien t tensor. C and C s are the fourth order bulk and surface stiffness tensors, e and e s are the third order bulk and surface piezo electric coupling tensors, κ and κ s are the bulk and surface dielectric p ermittivit y tensors resp ectiv ely . g is the sixth order strain gradient elasticit y tensor. µ is the fourth order flexo electric tensor which represen ts com bination of (a) strain-p olarization gradient coupling and (b) strain gradien t-p olarization coupling. The physical stress, σ and electric displacement, D can b e obtained from the bulk energy density function as, σ ij = ∂ U b ∂ ε ij −  ∂ U b ∂ η ij k  ,k = C ij kl ε kl − e ij k E k + µ ij kl E k,l − g ij klmn η lmn,k (2) D i = − ∂ U b ∂ E i = e ij k ε j k + µ ij kl ε j k,l + κ ij E j (3) The surface mechanical stress, σ s and surface electric displacement, D s can b e obtained from the surface energy densit y function as, σ s ij = ∂ U s ∂ ε s ij = τ s ij + C s ij kl ε s kl − e s ij k E s k (4) D s i = − ∂ U s ∂ E s i = − ω s i + e s ij k ε s j k + κ s ij E s j (5) 3 The total potential energy , Π can be written in terms of in ternal energy in the bulk, Π bulk , internal energy in the surface, Π s and w ork 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 displacemen t and electric p otential resp ectiv ely . t is the surface traction on Γ u , b is the prescrib ed b o dy force and q is the surface charge densit y on Γ φ . Γ u and Γ φ are the Neumann b oundary for mechanical displacement and electric potential resp ectively . The weak form of the equilibrium equations can be obtained by finding u ∈ { u = ¯ u on Γ d u , u ∈ H 2 (Ω) } and φ ∈ { φ = ¯ φ on Γ d φ , φ ∈ H 2 (Ω) } suc h 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 Γ d u , δ u ∈ H 2 (Ω) } and δ φ ∈ { δ φ = 0 on Γ d φ , δ φ ∈ H 2 (Ω) } . Γ d u and Γ d φ are the Diric hlet b oundary for mechanical displacement and electric potential resp ectively . 3. Mesh free form ulation for flexo electricit y The numerical discretization of the gov erning partial differential equations of flexo electricit y requires C 1 con tinuous basis functions for a Galerkin method. In the presen t work, we utilize a meshfree method with compactly supp orted radial basis function (CSRBF) shap e functions. Popular radial basis functions [30] include the Multi-Quadrics, Gaussian and Thin Plate Splines. These radial basis functions are globally supp orted and their accuracy highly dep ends on the condition n um b er of the collo cation matrix. Ho wev er, 4 the collo cation matrix will be a sparse matrix, well conditioned and compactly supp orted if we adopt CSRBF shap e functions. The W endland t yp e CSRBF s with C 2 and C 4 con tinuit y prop osed in [31] are, f ( r ( x, y )) = max  0 , (1 − r ) 4  (4 r + 1) ∈ C 2 (11) f ( r ( x, y )) = max  0 , (1 − r ) 6  (35 r 2 + 18 r + 3) ∈ C 4 (12) where r ( x, y ) is giv en b y , r i ( x, y ) = d i R = p ( x − x i ) 2 + ( y − y i ) 2 R (13) where d i is the distance of a p oin t of interest ( x, y ) from a knot at ( x i , y i ) . The dimension of supp ort domain, R is giv en b y R = αd c , where α is the shap e parameter and d c is the av erage nodal spacing. It is to b e noted that the v alue of r i lies betw een 0 and 1. An appro ximation for a general function can b e 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 co efficien ts resp ectiv ely , f T ( x ) = [ f 1 ( x ) , f 2 ( x ) , ...f n ( x )] (15) a T = [ a 1 , a 2 , ...a n ] . (16) In Equations 15 and 16, the v ariable n stands for the num b er of no des in the supp ort domain of the point of in terest. Here, p ( x ) and b are the vector of polynomial basis functions and co efficients resp ectiv ely , p T ( x ) = [ p 1 ( x ) , p 2 ( x ) , ...p m ( x )] (17) b T = [ b 1 , b 2 , ...b m ] . (18) In Equations (17) and (18), the v ariable m stands for the n umber of terms of p olynomial basis. The co efficient v ectors a and b can b e obtained by solving the follo wing algebraic equations " A P m P T m 0 # ( a b ) = ( U 0 ) , (19) where U is a vector of no dal v alues of function u h ( x ) and matrices A and P m are, A =     f 1 ( x 1 ) · · · f n ( x 1 ) . . . . . . . . . f 1 ( x n ) · · · f n ( x n )     (20) P m =     p 1 ( x 1 ) · · · p m ( x 1 ) . . . . . . . . . p 1 ( x n ) · · · p m ( x n )     . (21) 5 F rom Equation (14), in terp olation of the no dal function v alues, U , at any point of interest, x can be written as, 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 shap e function, N ( x ) is giv en b y , N ( x ) = f T ( x ) S a + p T ( x ) S b (23) where S a = A − 1 [1 − P m S b ] and S b = [ P T m A − 1 P m ] − 1 P T m A − 1 . The p olynomial basis functions of linear order are added to the radial basis functions in order to ensure that the shap e functions possess C 1 consistency . The vector p ( x ) given in Equation (17) can b e rewritten suc h that m = 3 as, p T ( x ) = [1 x y ] (24) The discrete form of the weak form ulation in Equation (10) using the meshfree shap e functions is as follo ws, δ u T  Z Ω B T u C B u d Ω  u + δ u T  Z Ω B T u e T B φ d Ω  φ + δ φ T  Z Ω B T φ eB u d Ω  u + δ u T  Z Ω H T u µ T B φ d Ω  u + δ φ T  Z Ω B T φ µH u d Ω  φ − δ φ T  Z Ω B T φ κB φ d Ω  φ + δ u T  Z Ω H T u g H u d Ω  u + δ u T  Z Γ B T u M T p τ s d Γ  − δ φ T  Z Γ P T ω s d Γ  + δ u T  Z Γ B T u M T p C s M p B u d Γ  u − δ u T  Z Γ B T u M T p e s T P B φ d Γ  φ + δ φ T  Z Γ B T φ P T e s M p B u d Γ  u − δ φ T  Z Γ B T φ P T κ s P B φ d Γ  φ = δ u Z Γ u N T ¯ t d Γ u + δ u Z Ω N T b d Ω − δ φ Z Γ φ N T q d Γ φ (25) where C , e , µ , κ and g are the matrix form of the tensors C ij kl , e ij k , µ ij kl , κ ij and g ij klmn resp ectiv ely and C s , e s , µ s and κ s are the matrix form of the tensors C s ij kl , e s ij k , µ s ij kl , and κ s ij resp ectiv ely . The gradien t and Hessian matrices in Equation 25 are defined as follo ws, B u =    N I ,x 0 0 N I ,y N I ,y N I ,x    (26) B φ = − " N I ,x N I ,y # (27) 6 H u =            N I ,xx 0 0 N I ,y x N I ,y x N I ,xx N I ,xy 0 0 N I ,y y N I ,y y N I ,xy            (28) where I = 1 , 2 ....n , n is the num b er of no des in the support domain of the p oin t of interest and this num ber can be different for differen t p oin ts of interest. The pro jection matrix is denoted as M P M P =    P 2 11 P 2 12 P 11 P 12 P 2 12 P 2 22 P 12 P 22 2 P 11 P 12 2 P 12 P 22 P 2 12 + P 11 P 22    (29) where the en tries of M P are from P , the tangential pro jection tensor giv en by I − n ⊗ n . Here, I refers to iden tity matrix of rank 2 and n is the outw ard normal v ector to the surface, Γ . The final system of equations can be written as follo ws, " K uu + K s uu K uφ + K s uφ K φu + K s φu K φφ + K s φφ # " u φ # = " F u + F s u F φ + F s φ # (30) where K uu = Z Ω B T u C B u d Ω + Z Ω H T u g H 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 s uu = 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 s uφ T K s φφ = − Z Γ B T φ P T κ s P B φ d Γ F s u = 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 11 C 12 0 C 12 C 22 0 0 0 C 66    µ = " µ 11 µ 12 0 0 0 µ 44 0 0 µ 44 µ 12 µ 11 0 # e = " 0 0 e 15 e 31 e 33 0 # κ = " κ 11 0 0 κ 22 # (32) g = l 0 2            C 11 C 12 0 0 0 0 C 12 C 22 0 0 0 0 0 0 C 66 0 0 0 0 0 0 C 11 C 12 0 0 0 0 C 12 C 22 0 0 0 0 0 0 C 66            (33) In Equation 33, the term l 0 is the length scale represen ting the size dep endency of strain gradient effects. 4. Mesh free form ulation for flexo electricit y including geometric nonlinearit y In this section, the prop osed meshfree formulation is extended to handle geometric nonlinearities in flex- o electric structures considering surface elasticit y . Saint V enant-Kirc hhoff material mo del is considered for flexo electric solids, the in ternal energy densit y giv en in Equation 1 is modified as, U = 1 2 S : G + 1 2 ˜ S . . . ˜ G − 1 2 D · E + 1 2 S s : G s (34) The total p otential energy , Π is giv en b y , Π = 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 displacement v ector; S s is the second Piola-Kirc hhoff 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 Green Lagrange surface strain tensor. u and φ are mechanical displacement and electric p oten tial resp ectively . t is the surface traction on Γ u , b is the prescrib ed bo dy force and q is the surface charge density on Γ φ . Γ u and Γ φ are the Neumann b oundary for mechanical displacement and electric p oten tial resp ectively . The constitutiv e equations of the assumed Sain t-V enant Kirchhoff material mo del are as follows S = C : G − e · E ˜ S = − µ · E + g . . . ˜ G D = e : G + µ . . . ˜ G + κ · E (36) 8 T aking the first v ariation of the total p oten tial 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) Eac h term in Equation 37 has to be linearized. The final expression obtained after linearizing each term in Equation 37 are subsequently presented. The in termediate steps are detailed in App endix A. A total Lagrangian formulation is presented such that all the integrals are p erformed on the undeformed configuration and deriv atives are with resp ect to the material co ordinates. The linearization of the term R Ω S : δ G d Ω in Equation 37 can b e obtained as L  Z Ω S : δ G d Ω  = 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 : [( ∇ 0 δ u ) T ( ∇ 0 (∆ 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 follo ws, L  Z Ω ˜ S . . . δ ˜ G d Ω  = 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 . . . [( ∇ 2 0 δ u )( ∇ 0 ∆ u ) + ( ∇ 2 0 ∆ u )( ∇ 0 δ u )] d Ω − Z Ω δ ˜ G . . . µ · ∆ E d Ω (41) The linearization of the term R Ω D · δ E d Ω in Equation 37 is as follows, L  Z Ω D · δ E d Ω  = 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  Z Γ S s : δ G s d Γ  = Z Γ ¯ S s : δ ¯ G s d Γ + Z Γ ∆( S s : δ G s ) d Γ (44) 9 Z Γ ∆( 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 · [( ∇ 0 δ u ) T ( ∇ 0 (∆ 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  Z Ω S : δ G d Ω  = δ u  Z Ω B T ˆ R d Ω  + δ u  Z Ω B T C B d Ω  ∆ u + δ u  Z Ω B T e B φ d Ω  ∆ φ + δ u  Z Ω H T 1 RH 1 d Ω  ∆ u (46) L  Z Ω ˜ S . . . δ ˜ G d Ω  = δ u  Z Ω H T D ˆ R D d Ω  + δ u  Z Ω H T D µ T B φ d Ω  ∆ φ + δ u  Z Ω H T 1 R T D H 2 d Ω  ∆ u + δ u  Z Ω H T 2 R D H 1 d Ω  ∆ u (47) L  Z Ω D · δ E d Ω  = − δ φ Z Ω B T φ ˆ D d Ω − δ φ  Z Ω B T φ µH u d Ω  ∆ u − δ φ  Z Ω B T φ eB d Ω  ∆ u + δ φ  Z Ω B T φ κB φ d Ω  ∆ φ (48) L  Z Γ S s : δ G s d Γ  = δ u  Z Γ B T M T p ˆ R s d Γ  + δ u  Z Γ B T M T p C s M p B d Ω  ∆ u δ u  Z Γ H T 1 P T n R s P n H 1 d Γ  ∆ u (49) All the matrices inv olv ed in Equations 46, 47, 48, 49 are presented in App endix B. The final algebraic form of linearization of Equation 37 can b e written as, K ∆ U = F ext − F int (50) where, K = " K uu K uφ K φu K φφ # ∆ U = " ∆ u ∆ φ # (51) 10 T able 1: Electromec hanical prop erties of STO Elastic Constan ts Dielectric constants Flexoelectric constants [32] C 11 =310 GP a κ 11 =2.66 C / ( GV − m ) µ 11 =-0.26 nC /m C 12 =115 GP a κ 33 =2.66 C / ( GV − m ) µ 12 =-3.74 nC /m C 22 =310 GP a µ 44 =-3.56 nC /m C 66 =54 GP a K uu = Z Ω B T C B d Ω + Z Ω H T 1 RH 1 d Ω + Z Ω H T 1 R T D H 2 d Ω + Z Ω H T 2 R D H 1 d Ω + Z Γ B T M p T C s M p B d Γ + Z Γ H T 1 P n R s P n H 1 d Γ (52) K uφ = Z Ω B T e T B φ d Ω + Z Ω H T D µ 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 T D ˆ 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 iterativ e sc heme. Solving this equation, giv es the deflection and voltage resp onses of flexo electric structures that undergo large deformations. 5. Numerical Examples In this section, the prop osed meshfree formulation is utilised to analyse flexo electric cantilev er b eams accoun ting for surface effects and also to study the influence of geometric nonlinearity on their voltage output. As an initial step, the meshfree formulation is v alidated by determining the energy conv ersion factor (ECF) of a cantilev er b eam. The ECF is giv en b y the ratio b etw een stored electrical energy and mechanical energy in the flexo electric structure. 5.1. V alidation: ECF A can tilev er beam sub jected to a mechanical p oin t load at the free end is analysed in order to v alidate the prop osed meshfree formulation. The cantilev er b eam has an asp ect ratio of 6. The Y oung’s mo dulus, Y is assumed to b e 100 GPa. F or v alidation, only the flexo electric constan t µ 12 and dielectric constant κ 22 are considered non-zero and they are assumed to be 10 nC /m and 1 nC /V m respectively . The b eam is discretized b y 121 × 21 no des with uniform spacing and the background mesh for n umerical in tegration is 120 × 20. The v ariation of the energy con version factor (ECF) with decreasing depth of the 1D flexo electric b eam is shown in Figure 1. The figure shows go o d agreement b etw een the numerical and analytical k 2 v alues. 11 Depth of beam, d (nm) 20 40 60 80 100 ECF 0 0.02 0.04 0.06 0.08 0.1 0.12 Analytical Meshfree Figure 1: The v ariation of ECF, k 2 with b eam depth for a one dimensional b eam mo del. The analytical energy conv ersion factor, k 2 for a 1D mo del is given in [33] as, E C F anl = k = r 12 κY  µ d  2 (56) 5.2. V alidation: T ub e mo del In order to further v alidate the proposed formulation, a flexo electric tub e made of STO is analysed with plane strain assumption. The tub e with an inner radius, r i of 10 µm and outer radius, r o of 20 µm is sub jected to a radial displacemen t of 0.045 µm and 0.05 µm at r i and r o resp ectiv ely [23, 24]. The tube is grounded along the inner face and a voltage of 1 V is applied along the outer face. The no dal distribution of the quarter mo del is as shown in Figure 2(a). The distribution of electric p otential obtained for a length scale, l 0 of 2 µm is shown in figure 2(b). The v ariation of electric p oten tial along the thic kness of the tub e is sho wn in figure 2(c). The analytical results for the flexo electric tube mo del with assumed material parameters is deriv ed by Mao et al . [23]. The results presented in Figure 2(c) shows go od agreemen t b etw een the analytical and n umerical results. 5.3. Pur e flexo ele ctricity A can tilev er b eam made of cubic STO is analysed in this section. The material prop erties of STO is giv en in T able 1. The b ottom and top face of the b eam are coated with electro des. The bottom electro de is grounded and the top electro de is free to ha v e a potential v alue. The dimension of the beam is 1200 × 100 nm. The beam is sub jected to a p oint load of 10 nN at the free end. The length scale, l 0 of g tensor is tak en as 0 (i.e.) g is not considered in this analysis. The beam has an almost linear v ariation of p oten tial along the b eam depth as sho wn in Figure 3. The p otential obtained at the top face is 30 m V. No w if we fix the asp ect ratio to b e 12, and reduce the b eam depth for instance to 40 nm, then the p oten tial at the top face is 74 m V. If the tensor, g is included in the analysis with a length 12 X, nm × 10 4 0 0.5 1 1.5 2 Y, nm × 10 4 0 0.5 1 1.5 2 (a) X, nm × 10 4 0.5 1 1.5 2 Y, nm × 10 4 0 0.5 1 1.5 2 -0.2 0 0.2 0.4 0.6 0.8 1 (b) Radius r, in µ m 10 12 14 16 18 20 Electric potential, φ in V -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Analytical Numerical (c) Figure 2: (a) No dal distribution for the quarter tub e mo del, (b) Electric p otential distribution across the tub e cross section, (c) The v ariation in output voltage along the radius of the tub e. X 0 200 400 600 800 1000 1200 Y 0 50 100 0 0.01 0.02 Figure 3: The p oten tial distribution in a flexo electric b eam made of STO. 13 Depth of beam, nm 40 50 60 70 80 90 100 ECF × 10 -4 0 0.5 1 1.5 2 g=0 g 0 Figure 4: The v ariation of ECF with depth of b eam including and excluding strain gradient tensor, g . The length scale is taken as 5 nm. scale, l 0 of 5 nm, then the potential obtained at the top face of the 40 nm b eam depth reduces to 52 m V. The change in energy conv ersion factor with depth of b eam for an asp ect ratio of 12, excluding and including g tensor is shown in Figure 4. The energy conv ersion factor for 40 nm b eam depth with and without including g tens or are 1.1e-4 and 1.61e-4 respectively . The inclusion of strain gradien t elasticit y increases the stiffness of the b eam, reduces the voltage obtained and as a result reduces the ECF s. Note that the difference b et ween the energy conv ersion factors with and without the g tensor increases with reduction in depth of the b eam. On the other hand, the influence of the length scale, l 0 on the output voltage is sho wn in Figure 5. It can b e seen that, when the length scale is increased from 0 to 5 nm, the output voltage of the 100 nm b eam reduces from 30 m V to 27.5 m V (6.3 % reduction). While for the same increase in length scale, the output voltage of the 40 nm b eam reduces from 73.6 m V to 51.7 m V (28.6 % reduction). 5.4. Flexo ele ctricity and surfac e effe cts In this section, we analyse a Zinc Oxide nano cantilev er b eam. The interpla y b etw een piezo electric, surface elastic, surface piezo electric and flexo electric effects is studied. Zinc Oxide is the ideal material for p erforming this study b ecause it is widely used in several nanoscale energy harvesters [34, 3] and studies on surface properties of Zinc Oxide [35] is a v ailable. The can tilever ZnO beam is of length, 120 nm and width, 15 nm. A p oint load of 10 nN is applied in x-direction at the mid-p oin t of the top face. The beam is fixed at the b ottom and free at the top. The b eam is poled along the length (y-direction). The b ottom end of the b eam is grounded. The elastic, piezo electric, surface elastic and surface piezoelectric properties of ZnO are given in T able 2 and T able 3. The residual stress, τ s and residual electric displacement, ω s are not considered in the study . The flexoelectric constant of ZnO, µ 11 , µ 12 and µ 44 are assumed to b e 2 nC/m, 2 nC/m and 0.5 nC/m resp ectiv ely . 14 Beam depth, nm 40 50 60 70 80 90 100 Voltage, mV 20 30 40 50 60 70 80 l 0 =0 l 0 =1 nm l 0 =2.5 nm l 0 =5 nm Figure 5: Output voltage for varying b eam depths for internal length scales of 0,1,2.5 and 5 nm. T able 2: Material prop erties of bulk ZnO Elastic Constan ts Piezo electric constants Dielectric constants C 11 =206 GP a e 31 =-0.58 C /m 2 κ 11 =0.0811 C / ( GV − m ) C 12 =117 GP a e 33 =1.55 C /m 2 κ 33 =0.112 C / ( GV − m ) C 22 =211 GP a e 15 =-0.48 C /m 2 C 66 =44.3 GP a T able 3: Material prop erties of ZnO surface Elastic Constan ts Piezo electric constants C s 11 =44.2 N/m e s 31 =-0.216 nC /m C s 12 =14.2 N/m e s 33 =0.451 nC /m C s 22 =35 N/m e s 15 =-0.253 nC /m C s 66 =11.7 N/m 15 (a) (b) Figure 6: The potential distribution in the ZnO b eam considering (a) Flexoelectricity , piezo electricit y and surface effects, (b) Pure Flexo electricit y . The p oten tial distribution across the b eam width is shown in Figure 6. The combination of bulk piezo elec- tricit y and surface elastic effect results in a p oten tial of +1.18 V to -1.18 V at the top face. The combination of bulk piezoelectricity and surface piezoelectric effect leads to a p oten tial v arying from +1.5 V to -1.5 V at the top face of the b eam. Finally , the combination of bulk piezo electricity , bulk flexo electricit y , surface elasticit y and surface piezoelectricity results in a p oten tial of +1.7 V to -1.7 V at the top face as sho wn in Figure 6(a). The v ariation of electric p oten tial considering only flexo electric effect is +0.3 V to -0.3 V at the top face (Figure 6(b)). The relativ e influence of the different phenomenon on the output v oltage is sho wn in Figure 7. The con tribution of flexo electric effect to output voltage is higher compared to the contribution of surface effects and the difference in the contributions to output voltage increases as the b eam width decreases. F or a 40 nm wide b eam, the flexo electric and surface effect con tributions are 7 % and 1 % resp ectiv ely . While for a width of 15 nm, the difference b etw een the contributions is higher, the flexo electric and surface effect con tributions are 18 % and 7 % resp ectiv ely . The change in ECF with width of b eam is shown in Figure 8. The pattern is similar to the one obtained for output voltage. The ECF for pure Piezo electricity and Piezo + Surface effects for 15 nm wide b eam are 0.0067 and 0.00728 resp ectively . While for the same beam width, the ECF considering Piezo + Surface effects + Flexo electricit y is 0.014. The p ercen tage contribution of flexo electricity and surface effects to the total ECF are 48 % and 8 % respectively . There is a discrepancy in the p ercentage con tribution of flexoelectricity to ECF and output voltage. This is b ecause the potential due to flexo electricity reaches its p eak near the fixed end and reduces significan tly along the length. So, though the flexo electric contribution to total ECF is 48 % , the contribution of flexo electricity to total voltage measured at the top face (y=120 nm) is only 18 % . 5.5. Flexo ele ctricity and surfac e effe cts: Ge ometric nonline arity In this section, the flexo electric resp onse in the nonlinear regime is studied. The nonlinearity emerges due to large deformation of the flexoelectric cantilev er b eam. The flexoelectric beam is assumed to be made of STO, in addition to flexoelectricity , the surface elasticity of STO is also considered. The material properties 16 Width of beam, nm 20 40 60 80 100 Voltage, V 1.2 1.3 1.4 1.5 1.6 1.7 1.8 e 0,e s =0,C s =0, µ =0 e 0,e s 0,C s 0, µ =0 e 0,e s 0,C s 0, µ 0 Figure 7: The v ariation of output voltage for ZnO b eam with width for three cases, Piezoelectricity , Piezo electricity+surface effects and Piezo electricity+Surface effects+Flexo electricit y . Width of beam, nm 20 40 60 80 100 ECF 0.006 0.008 0.01 0.012 0.014 e 0, C s = 0 , e s = 0 , µ = 0 e 0, C s 0 , e s 0 , µ = 0 e 0 C s 0 , e s 0 , µ 0 Figure 8: The v ariation of ECF for ZnO beam with width for three cases, Piezo electricity , Piezo electricity+Surface effects and Piezoelectricity+Surface effects+Flexo electricit y . 17 of STO are giv en in T able 1. The surface elastic constan ts of STO are not av ailable in literature and are assumed to be, C s 11 = C s 22 = 310 GP a , C s 12 = 115 GP a and C s 66 = 54 GP a . The b eam is sub jected to mec hanical deformation by a p oin t load of 10 nN at the free end. The load is applied in increments of 1 nN. In eac h increment, the tangent stiffness matrix is determined and the Newton-Raphson method is adopted to minimize the residual. T w o different b eams each of thickness 50 nm and 30 nm resp ectiv ely are analysed. The aspect ratio of the b eams is fixed as 12. The fixed end of the cantilev er is grounded. The v ariation 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 v oltage due to flexo electricit y and surface elasticit y for 50 nm thick b eam is 7.8 m V. The ratio b et w een nonlinear and linear v oltage is 0.9 (i.e.) the nonlinear voltage deviates from the linear resp onse by 10% . If surface effects are not considered, then the final output voltage is 8.5 m V. In case of 30 nm thic k b eam, as sho wn in Figure 9(b), the final output flexoelectric voltage is 11 m V and 12 m V considering and ignoring surface effects respectively . The v ariation of free end deflection with load steps is shown in Figure 10. The free end deflection of 50 nm thick b eam after ten load steps is 68 nm (Figure 10(a)). The ratio b etw een nonlinear and linear deflection for a load v alue of 10nN for 50 nm thic k b eam is 0.89 (i.e.) a deviation of nonlinear displacemen t is 11% from linear displacemen t. In the absence of surface effects, the final output deflection of 50 nm thic k b eam is 75 nm. F or the case of 30 nm thic k beam, as shown in Figure 10(b), the final deflection is 56 nm and 63.5 nm considering and ignoring surface effects resp ectiv ely .F rom Figures 9 and 10, it is clear that the b eam b ecomes stiffer in the presence of surface elasticit y , whic h in turn leads to reduction in output v oltage. In case of 50 nm thick b eam, due to surface elasticity the absolute v alue of final voltage reduces from 8.5 m V to 7.8 m V. While for 30 nm thick beam, the voltage reduction is from 12 m V to 11 m V. The v ariation of energy con version factor with load steps is shown in Figure 13. The rate of increase in ECF with load steps is higher for a 30 nm b eam compared to the 50 nm beam (Asp ect ratio = 12). It is to b e noted that, i n the absence of geometric nonlinearit y , the energy conv ersion factor is a constan t v alue and remains independent of the applied loads, whereas the consideration of geometric nonlinearity has led to change in ECF v alue with load incremen ts. The contour plot showing the normal strain in x-direction is given in Figures 11 and 12. Figure 11 shows that at the first load step, the nonlinear strain terms are negligible and the Green-Lagrange strain con tour and linear strain con tour are quite similar. While the strain contour at the 10 th load step given in Figure 12 sho ws that nonlinear terms play significan t role in determining the Green-Lagrange strain. Consequently , the linear strain con tour and Green Lagrange strain contours b ecome dissimilar. The studies p erformed in this section sho w that the surface elasticit y can reduce the output flexo electric v oltage. F or instance, it is observ ed that in case of a 30 nm thick b eam the reduction in the final voltage is 1 m V due to surface elastic effects (Figure 9). The surface elastic effects can increase with increase in the surface area of the beam. Therefore it is w orthwhile to study a tap ered b eam model with inclined top and bottom faces and compare their response with the rectangular b eams studied previously . Considering an av erage thickness of 30 nm, a tapered can tilever b eam, T B 1 , of length 360 nm and of depth, d l =40 nm at x=0, d r =20 nm at x=360 nm is sub jected to a load increment of 1 nN o ver ten load steps. A t the tenth load step, the output flexo electric voltage in T B 1 reduces from 15 m V to 13 m V due to surface elastic effects as sho wn in Figure 15(a). Secondly , a tapered cantilev er b eam T B 2 of length 360 nm and of depth, d l =45 nm at x=0, d r =15 nm at x=360 nm is considered. The final voltage reduces from 12 m V to 10.4 m V as shown in Figure 15(b). In b oth the tap ered beams T B 1 and T B 2 with a verage thickness of 30 nm, the v oltage reduces b y 13 . 3% due to surface elastic effects. This is higher than the reduction of 18 Load steps 2 4 6 8 10 Electric potential, φ in V -0.01 -0.008 -0.006 -0.004 -0.002 0 nonlinear ( µ 0, C s 0) nonlinear ( µ 0, C s =0) (a) Load steps 2 4 6 8 10 Electric potential, φ in V -0.012 -0.01 -0.008 -0.006 -0.004 -0.002 0 nonlinear ( µ 0, C s 0) nonlinear ( µ 0, C s =0) (b) Figure 9: The v ariation in voltage with load incremen ts for flexo electric b eam made of STO of depth (a) 50 nm and (b) 30 nm (Aspect ratio=12). Load steps 2 4 6 8 10 Deflection, nm 0 20 40 60 80 nonlinear( µ 0, C s 0) nonlinear( µ 0, C s =0) (a) Load steps 2 4 6 8 10 Deflection, nm 0 10 20 30 40 50 60 70 nonlinear( µ 0, C s 0) nonlinear( µ 0, C s =0) (b) Figure 10: The load deflection curve for flexo electric b eam made of STO of depth (a) 50 nm and (b) 30 nm (Asp ect ratio=12). 8 . 5% (from 11.8 m V to 10.8 m V) in 30 nm thic k rectangular beam. Therefore, as exp ected the influence of surface effects increases with increase in surface area of the b eam and the negativ e influence on flexo electric v oltage is higher for a tap ered b eam compared to a rectangular be am of same v olume. In summary , the results obtained in this section show that considering nonlinear terms in strain and gradien t of strain, is more essential as the influence of flexo electricit y on output voltage gets significant in nanoscale. W e conclude that the geometric nonlinearit y cannot b e ignored if one analyses flexo electric b eams of dimensions of under 100 nanometers when sub jected to loads in the range of 10 nNs. It is to b e noted that the flexo electric material, STO used in this example has a lesser flexo electric constan t of only 1.4 V, while flexo electric constan t of a dielectric material can even b e upto 10 V based on the theoretical upp er limit estimated by Kogan et al. [8]. 19 (a) (b) (c) Figure 11: The strain con tour ( G 11 , ε 11 , η 11 ) for flexo electric beam made of STO of depth 50 nm at load step 1 (a) Green Lagrange strain (b) Linear strain (c) Nonlinear part of strain. (a) (b) (c) Figure 12: The strain con tour ( G 11 , ε 11 , η 11 ) for flexoelectric beam made of STO of depth 50 nm at load step 10 (a) Green Lagrange strain (b) Linear strain (c) Nonlinear part of strain. 20 Load steps 2 4 6 8 10 ECF × 10 -6 0 1 2 3 4 5 6 d = 50 nm d = 30 nm Figure 13: The v ariation in ECF with load increments for flexo electric b eam made of STO of depth 50 nm and 30 nm (Asp ect ratio =12). d l d r l = 600 nm Figure 14: T apered flexoelectric beam model. Load steps 2 4 6 8 10 Electric potential, φ in V -0.015 -0.01 -0.005 0 nonlinear( µ 0, C s 0) nonlinear( µ 0, C s =0) (a) Load steps 2 4 6 8 10 Electric potential, φ in V -0.012 -0.01 -0.008 -0.006 -0.004 -0.002 0 nonlinear( µ 0, C s 0) nonlinear( µ 0, C s =0) (b) Figure 15: The v ariation in v oltage with load incremen ts for flexo electric tap ered b eam made of STO of type (a) T B 1 and (b) T B 2 . 21 6. Conclusion A CSRBF based meshfree formulation is presented in this pap er to handle geometric nonlinearity in flex- o electric structures. In addition to flexoelectricity , the surface effects are also considered in the analysis of nano-sized t w o dimensional structures. Flexo electric b eams made of cubic STO, whic h is non-piezo electric and ZnO, which is piezoelectric are analysed. The meshfree analysis sho ws that for ZnO, the con tribution of surface effects to the output v oltage of a nanosized cantilev er structure (of width 15 nm) is smaller compared to the contribution of flexo electricit y . The analysis of flexo electric nanostructures undergoing large deformation shows that the difference betw een nonlinear and linear flexo electric voltage increases with reduction in b eam depth. The surface elastic effects stiffen the beam leading to reduction in output flexo electric voltage. The influence of surface elasticity is higher in tap ered b eams compared to rectangular b eams. In future, the presented form ulation will b e ex- tended to study nonlinear flexo electricit y under dynamic excitations and the influence of nonlinearity on resp onse bandwidth of nanosized flexoelectric energy harv esters will be inv estigated. App endix A. Intermediate Steps in the deriv ation of nonlinear meshfree form ulation for flex- o electricit y The terms G , δ G , ∆ δ G and ∆ S in Equation 38 are as follows, G = 1 2 ( u i,j + u j,i + u k,i u k,j ) (A.1) δ G = 1 2 ( δ u i,j + δ u j,i + δ u k,i u k,j + u k,i δ u k,j ) (A.2) ∆ δ G = 1 2 ( δ 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 follo ws, ˜ G = 1 2 ( u i,j k + u j,ik + u k,ij u k,j + u k,i u k,j i ) δ ˜ G = 1 2 ( δ u i,j k + δ u j,ik + δ u k,ij u k,j + u k,ij δ u k,j + δ u k,i u k,j i + u k,i δ u k,j i ) ∆ δ ˜ G = 1 2 ( δ u k,ij ∆ u k,j + ∆ u k,ij δ u k,j + δ u k,i ∆ u k,j i + ∆ u k,i δ u k,j i ) = δ u k,ij ∆ u k,j + ∆ u k,ij δ u k,j ˜ S = − µ · E (A.5) 22 The 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 follo ws, G s = 1 2 P ij ( u i,j + u j,i + u k,i u k,j ) P j i (A.7) δ G s = 1 2 P ij ( δ u i,j + δ u j,i + δ u k,i u k,j + u k,i δ u k,j ) P j i (A.8) ∆ δ G s = 1 2 P j i ( δ u k,i ∆ u k,j + ∆ u k,i δ u k,j ) P j i = P j i δ u k,i ∆ u k,j P j i (A.9) S s = C s : G s ∆ S s = C s : ∆ G s (A.10) App endix B. Matrices - nonlinear meshfree formulation for flexoelectricity The expressions defining the matrices B , B φ , H 1 , H 2 , 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 0 0 N I ,y N I ,y N I ,x    + A H 1 (B.1) A =    ∂ u I x ∂ x 0 ∂ u I y ∂ x 0 0 ∂ u I x ∂ y 0 ∂ u I y ∂ y ∂ u I x ∂ y ∂ u I x ∂ x ∂ u I y ∂ y ∂ u I y ∂ x    (B.2) H 1 =       N I ,x 0 N I ,y 0 0 N I ,x 0 N I ,y       (B.3) B φ = " N I ,x N I ,y # (B.4) H D = H u + A D H 2 (B.5) 23 A D =            ∂ u I x ∂ x 0 0 0 ∂ u I y ∂ x 0 0 0 0 ∂ u I x ∂ x 0 0 0 ∂ u I y ∂ x 0 0 0 0 ∂ u I x ∂ y 0 0 0 ∂ u I y ∂ y 0 0 0 0 ∂ u I x ∂ y 0 0 0 ∂ u I y ∂ y ∂ u I x ∂ y 0 ∂ u I x ∂ x 0 ∂ u I y ∂ y 0 ∂ u I y ∂ x 0 0 ∂ u I x ∂ y 0 ∂ u I x ∂ x 0 ∂ u I y ∂ y 0 ∂ u I y ∂ x            (B.6) H 2 =                 N I ,xx 0 N I ,xy 0 N I ,y x 0 N I ,y y 0 0 N I ,xx 0 N I ,xy 0 N I ,y x 0 N I ,y y                 (B.7) R =       S 11 S 12 0 0 S 12 S 22 0 0 0 0 S 11 S 12 0 0 S 12 S 22       (B.8) ˆ R = h S 11 S 22 S 12 i T (B.9) R D =                 S 111 S 121 0 0 S 122 S 212 0 0 S 121 S 122 0 0 S 212 S 222 0 0 0 0 S 111 S 121 0 0 S 211 S 212 0 0 S 121 S 122 0 0 S 212 S 222                 (B.10) ˆ R D = h S 111 S 211 S 122 S 222 S 121 S 212 i T (B.11) ˆ D = " D 1 D 2 # (B.12) ˆ R s = h τ s 11 + S s 11 τ s 22 + S s 22 τ s 12 + S s 12 i T (B.13) R s =       τ s 11 + S s 11 τ s 12 + S s 12 0 0 τ s 12 + S s 12 τ s 22 + S s 22 0 0 0 0 τ s 11 + S s 11 τ s 12 + S s 12 0 0 τ s 12 + S s 12 τ s 22 + S s 22       (B.14) 24 P n =       P 11 P 12 0 0 P 21 P 22 0 0 0 0 P 11 P 12 0 0 P 21 P 22       (B.15) A c knowledgmen ts X. Zh uang ac kno wledges the support from State Key Lab oratory of Structural Analysis for Industrial Equipmen t (GZ1607) and National Science F oundation of China (11772234). References [1] R. Mbarki, J. Haskins, A. Kinaci, T. Cagin, T emp erature dep endence of flexo electricity in B aT iO 3 and S r T iO 3 p ero vskite nanostructures, Ph ysics Letters A 378 (2014) 2181–2183. [2] S. Priy a, H.-C. Song, Y. Zhou, R. V arghese, A. Chopra, S.-G. Kim, I. Kanno, L. W u, D. S. Ha, J. Ryu, et al., A review on piezo electric energy harvesting: materials, metho ds, and circuits, Energy Harvesting and Systems 4 (2017) 3–39. [3] Z. L. W ang, J. Song, Piezo electric nanogenerators based on zinc oxide nano wire arra ys, Science 312 (2006) 242–246. [4] N. Sharma, R. Maranganti, P . Sharma, On the p ossibility of piezoelectric nano composites without using piezo electric materials, Journal of the Mechanics and Ph ysics of Solids 55 (2007) 2328–2350. [5] U. K. Bhask ar, N. Banerjee, A. Abdollahi, Z. W ang, D. G. Schlom, G. Rijnders, G. Catalan, A flexo- electric microelectromechanical system on silicon, Nature nanotec hnology 11 (2016) 263–266. [6] V. Mashkevic h, Electrical, optical, and elastic prop erties of diamond-t yp e crystals ii. lattice vibrations with calculation of atomic dip ole moments, So viet Ph ysics Jetp 5 (1957). [7] A. K. T agantsev, Piezoelectricity and flexo electricity in crystalline dielectrics, Physical Review B 34 (1986) 5883. [8] S. M. Kogan, Piezo electric effect under an inhomogeneous strain and acoustic scattering of carriers in crystals, Fiz. T verd T ela 5 (1963) 2829–2831. [9] W. Ma, L. E. Cross, Large flexo electric polarization in ceramic Lead Magnesium Niobate, Applied Ph ysics Letters 79 (2001) 4420–4422. [10] W. Ma, L. E. Cross, Flexo electric polarization of Barium Strontium Titanate in the paraelectric state, Applied Ph ysics Letters 81 (2002) 3440–3442. [11] W. Ma, L. E. Cross, Flexo electricity of Barium Titanate, Applied Ph ysics Letters 88 (2006) 232902. [12] R. Marangan ti, N. D. Sharma, P . Sharma, Electromec hanical coupling in nonpiezo electric materials due to nanoscale nonlo cal size effects: Green’s function solutions and em b edded inclusions, Ph ysical Review B 74 (2006) 014110. 25 [13] P . Zubk o, G. Catalan, A. Buckley , P . W elche, J. Scott, Strain-gradien t-induced p olarization in S rT iO 3 single crystals, Ph ysical Review Letters 99 (2007) 167601. [14] Z. Y an, L. Y. Jiang, Flexo electric effect on the electro elastic resp onses of b ending piezo electric nanob eams, J. Appl. Phys. 113 (2013) 194102. [15] C. C. L iu, S. L. Hu, S. P . Shen, Effect of flexo electricit y on electrostatic p oten tial in a b ent piezo electric nano wire, Smart Mater. Struct. 21 (2012) 115024. [16] S. Zhang, H. Y ao, W. F an, Y. Hao, X. W u, D. Hou, Effects of flexoelectricity and surface elasticit y on piezo electric potential in a bent ZnO nano wire, in: IOP Conference Series: Materials Science and Engineering, v olume 167, IOP Publishing, p. 012023. [17] H. Chen, A. Soh, Y. Ni, Phase field mo deling of flexo electric effects in ferroelectric epitaxial thin films, A cta Mec hanica 225 (2014) 1323–1333. [18] H. Chen, S. Zhang, A. Soh, W. Yin, Phase field mo deling of flexo electricit y in solid dielectrics, Journal of Applied Physics 118 (2015) 034106. [19] A. Abdollahi, C. Peco, D. Millán, M. Arroy o, I. Arias, Computational ev aluation of the flexo electric effect in dielectric solids, Journal of Applied Ph ysics 116 (2014) 093502. [20] A. Ab dollahi, D. Millán, C. Peco, M. Arroy o, I. Arias, Revisiting p yramid compression to quantify flexo electricit y: A three-dimensional sim ulation study , Physical Review B 91 (2015) 104103. [21] H. Ghasemi, H. S. Park, T. Rab czuk, A lev el-set based IGA form ulation for topology optimization of flexo electric materials, Computer Methods in Applied Mec hanics and Engineering 313 (2017) 239–258. [22] S. S. Nanthakumar, X. Zh uang, H. S. P ark, T. Rab czuk, T op ology optimization of flexoelectric struc- tures, Journal of the Mechanics and Physics of Solids 105 (2017) 217 – 234. [23] S. Mao, P . K. Purohit, N. Ara v as, Mixed finite-elemen t form ulations in piezo electricit y and flexoelec- tricit y , Proceedings of the Roy al So ciety of London A: Mathematical, Physical and Engineering Sciences 472 (2016). [24] F. Deng, Q. Deng, W. Y u, S. Shen, Mixed finite elements for flexoelectric solids, Journal of Applied Mec hanics 84 (2017) 081004. [25] S. S. Nanthakumar, N. V alizadeh, H. S. Park, T. Rab czuk, Surface effects on shap e and top ology optimization of nanostructures, Computational Mec hanics 56 (2015) 97–112. [26] S. S. Nan thakumar, T. Lahmer, X. Zhuang, H. S. Park, T. Rab czuk, T opology optimization of piezo- electric nanostructures, Journal of the Mechanics and Ph ysics of Solids 94 (2016) 316 – 335. [27] M. F arsad, F. J. V ernerey , H. S. Park, An extended finite elemen t/level set metho d to study surface effects on the mec hanical b eha vior and prop erties of nanomaterials, International Journal for Numerical Metho ds in Engineering 84 (2010) 1466–1489. [28] J. Y vonnet, L. Liu, A numerical framework for mo deling flexo electricity and maxwell stress in soft dielectrics at finite strains, Computer Metho ds in Applied Mechanics and Engineering 313 (2017) 450– 482. 26 [29] S. Shen, S. Hu, A theory of flexo electricit y with surface effect for elastic dielectrics, Journal of the Mec hanics and Ph ysics of Solids 58 (2010) 665–677. [30] G.-R. Liu, Meshfree metho ds: moving b eyond the finite elemen t method, T a ylor & F rancis, 2009. [31] H. W endland, Piecewise p olynomial, p ositive definite and compactly supp orted radial functions of minimal degree, A dv ances in computational Mathematics 4.1 (1995) 389–396. [32] R. Marangan ti, P . Sharma, A tomistic determination of flexo electric prop erties of crystalline dielectrics, Ph ysical review B 80 (2009) 054109. [33] M. S. Ma jdoub, P . Sharma, T. Cagin, Enhanced size-dep endent piezo electricity and elasticit y in nanos- tructures due to the flexo electric effect 77 (2008) 125424. [34] G. Zhu, R. Y ang, S. W ang, Z. L. W ang, Flexible high-output nanogenerator based on lateral zno nanowire arra y , Nano letters 10 (2010) 3151–3155. [35] M.-T. Hoang, J. Y vonnet, A. Mitrushc henk ov, G. Chambaud, First-principles based m ultiscale mo del of piezoelectric nanowires with surface effects, Journal of Applied Physics 113 (2013) –. 27