A gradient electromechanical theory for thin dielectric curved beams considering direct and converse flexoelectric effects

This work presents the development of a unified gradient electromechanical theory for thin flexoelectric beams considering both direct and converse flexoelectric effects. The two-way coupled electromechanical theory is developed starting from 3D variationals formulation by considering an electric field-strain-based free energy function. The formulation incorporates mechanical as well as electrical size effects. The coupled 3D theory is specialized to isotropic materials, and a 1D beam theory for composite flexoelectric curved beams is derived using the classical Kirchhoff assumptions. The beam theory is solved using a novel C2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^2$$\end{document} continuous finite element framework for different loading and boundary conditions. Our finite element results are verified with analytical solutions for a simply supported flexoelectric beam operating in both actuator and sensor modes. The results are also compared with existing literature for the special case of a passive micro-beam. Our computational framework is subsequently used to perform various parametric studies to analyze the effect of electrical and mechanical length scale parameters, geometric parameters like the radius of curvature, flexoelectric layer thickness etc., on the response of the beam. Also, contribution of converse flexoelectricity in the overall response of the flexoelectric beam is compared with that of the direct effect. Our simulation results predict that the converse effect is significant (≈\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx $$\end{document} 10–25% of the direct effect) for a wide range of thickness and length scale parameter values. It is also observed that the effective electromechanical coupling coefficient, calculated in terms of the voltage developed across the flexoelectric layer thickness, is higher in flexoelectric materials compared to piezoelectric materials at smaller length scales (thickness of the order of a few microns). Our simulation results also agree well with the trends observed in recent experiments [1].


Introduction
Flexoelectric materials exhibit spontaneous polarization in the presence of strain gradient. This electromechanical phenomenon, known as the direct flexoelectric effect, is observed in dielectric materials. Conversely, these materials also exhibit deformation in the presence of polarization gradient. This effect is known as the converse flexoelectric effect. Flexoelectric effect is often compared with the piezoelectric effect which is found in non-centrosymmetric dielectric materials. Piezoelectric materials exhibit large electromechanical coupling at macro-scale, whereas flexoelectric effect is significant only in thin structures and is comparable to piezoelectric coupling at micro/nanoscales [2,3]. Unlike the piezoelectric effect, flexoelectric effect is also exhibited by centro-symmetric crystalline materials [4]. Thus, flexoelectric materials are used as micro-transducers in MEMS and NEMS applications [5,6]. Specifically, they have been explored in applications like micro-actuators [7,8], structural health monitoring [9], vibration control [10], curvature sensors [11], and energy harvesting [12][13][14], to name a few.
Flexoelectric effect has been studied extensively in the recent literature [15][16][17][18]. Majority of the existing literature deals with modeling and applications of the direct flexoelectric effect. Converse flexoelectric effect is typically ignored as it is assumed to be insignificant compared to the direct effect. However, Page 2 of 33 Y. S. Joshan and S. Santapuri ZAMP recent studies have shown that at micro/nanoscales converse flexoelectric effect gives rise to a large localized electroelastic response in all dielectric materials. Abdollahi et al. [1] found that neglecting converse flexoelectricity results in an incorrect estimation of piezoelectric coefficient using piezoresponse force microscopy. Moreover, analogous to the mechanical size effects that result in a stiffer material at microscale, recent studies have demonstrated an increase in the electrical permittivity of the material at smaller scales [19]. Motivated by these recent studies, here we aim to develop a generalized two-way coupled size-dependent flexoelectric framework, that can explain these observation and enable modeling of a wider range of actuator and sensor devices by incorporating both direct and converse flexoelectric effects.
Continuum strain-gradient theories are typically developed using variational principles wherein the choice of independent variables affects the mathematical form of the governing equations, boundary conditions and the constitutive equations [20]. As one of the earliest works, a three-dimensional continuum flexoelectric theory was developed by Sahin and Dost [21] in which polarization and polarization gradient were chosen as the independent electrical variables along with the deformation gradient and double deformation gradient as the mechanical variables. The theory was linearized for the case of infinitesimal strains and small polarization charges. Extending this work, Maranganti et al. [22] presented a linear flexoelectric theory for centro-symmetric dielectrics assuming free energy as a function of double deformation gradient and polarization. This theory considers two material length scale parameters in constitutive equations accounting for dilatational and rotational gradients. Majdoub et al. [4] developed a linear size-dependent piezo-flexoelectric theory for nano-structures. This work analyzed a cantilever beam using molecular dynamics and found that the effective electroelastic coefficient increases at nanoscales due to the flexoelectric effect. In the above-mentioned developments of continuum flexoelectric polarization-strain gradient theory, the electric forces due to Maxwell stresses were ignored. Hu and Shen [23] considered Maxwell stresses in the flexoelectric theory for nano-dielectrics using the variational principle. A micromorphic continuum theory to model electro-magneto-elastic solids was developed by Romeo [24]. This theory was used to analyze polarization in dielectric medium [25]. By restricting the analysis to micropolar isotropic dielectrics, the flexoelectric contribution to polarization vector was derived.
One-dimensional beam structures are the building blocks of MEMS and NEMS devices and have been extensively studied in the last decade. Liu et al. [15] analyzed the effect of flexoelectricity on cantilever piezoelectric nano-wires. It was found that flexoelectricity significantly affects the electrostatic potential at nanoscales, even for small values of flexoelectric coefficients. Yan and Jiang [3] presented Euler-Bernoulli beam model for piezo-flexoelectric nano-beams. This work analyzed the effect of different boundary conditions on the bending of flexoelectric nano-beams. Deng et al. [16] developed an Euler-Bernoulli beam model for flexoelectric energy harvesters. It was shown that the conversion efficiency of the flexoelectric energy harvesters increases with decrease in the sample size. Wang and Wang [18] extended the flexoelectric energy harvester model to piezo-flexoelectric unimorphs. The authors found that the contribution of flexoelectric effect on the power output is larger than that of piezoelectric effect at nanoscale. Yue et al. [26] used Timoshenko beam model to analyze bending and free vibration response of simply supported piezoelectric beams considering flexoelectricity and surface effects. Qi et al. [27] developed a size-dependent model considering both direct and converse flexoelectric effects using a simplified strain gradient elasticity theory. They analyzed a bilayer cantilever beam and found that the strain gradients and flexoelectricity significantly influence the bending response, when beam thickness is comparable to length scale parameters. Wang and Wang [28] developed a nonlinear energy harvester model for piezo-flexoelectric micro and nano-beams, and concluded that the voltage output contributed by flexoelectric effect is five times higher than that by piezoelectric effect at nanoscale. Chu et al. [29] analyzed the effect of gradation of material properties on static and free vibration response of functionally graded flexoelectric nano-beams using modified strain gradient theory. They found that flexoelectric response depends on the volume fraction and gradation of the material. Rahmati et al. [30] investigated nonlinear bending response of soft electrets considering flexoelectricity and piezoelectricity. The modeling results were used to discuss design of possible electret structures that can generate large transverse piezoelectricity and flexoelectricity.
The beam models mentioned above are based on polarization-strain gradient flexoelectricity theory which considers free energy as a function of polarization and strain gradient. However, compared to polarization, electric field is easier to control and measure in practical applications. Thus, electric fieldstrain gradient-based flexoelectric theories have recently been developed in literature [31,32]. Liu et al. [33] developed a electromechanical theory considering electric field and electric field gradient as independent variables in the energy formulation and analyzed effective piezoelectric response considering flexoelectric effect under a scanning probe using mixed finite element method. They observed that flexoelectric effect contributes around 20% to the effective electroelastic response. Experimental and computational analysis of flexoelectric structures incorporating converse flexoelectricity has been studied in recent literature [34][35][36][37]. Furthermore, electric field-based formulations have been used to analyze flexoelectric beams [38][39][40][41] for different transducer applications. However, the existing electric field-strain gradient flexoelectric theories do not consider the higher-order effects due to electrical field gradient in the constitutive relations and also ignore the electrical size effects that contribute to an increase in electric permittivity of the material at lower scales [19]. Recent experimental studies suggest the need for a generalized gradient electroelastic theory that can corroborate the experimental observations and enable design of a wide range of sensor and actuator devices.
To this end, in this work, a generalized linear flexoelectric theory is developed considering infinitesimal strain, strain gradient, electric field, and electric field gradient variables, starting from the 3D variational principle. The generalized framework is specialized to a two-way coupled linear flexoelectric theory which incorporates the often ignored converse flexoelectric effect. Further specializing to isotropic dielectric materials, the mechanical size effects have been modeled using three length scale parameters. Relationship between the length scale parameters and higher-order constitutive equations is obtained as discussed in Lam et al. [42]. Two additional electrical length scale parameters are introduced in this work, to model the electrical size effects. The theory is applied to a composite flexoelectric curved beam and solved as a 1D electroelastic boundary value problem. Displacement of the beam is modeled using the classical beam theory assumptions. Also, a quadratic variation of the electrostatic potential is assumed across the thickness of the flexoelectric layer. Coupled 1D governing equations and boundary conditions for the curved beam are derived starting from the 3D variational formulation. An electromechanically coupled C 2 continuous finite element (FE) framework is subsequently developed to analyze the beam structure for different boundary and loading conditions. The flexoelectric curved beam is analyzed in two different modes: (i) sensor mode, wherein an electrostatic potential is generated due to applied mechanical loads, and (ii) actuator mode, wherein the beam deforms in response to an applied electric field. Our FE results are verified with analytical solution available in literature for a passive beam. For both sensor and actuator modes, analytical solution for a simply supported flexoelectric curved beam is developed using the Fourier series method to assess the accuracy of our coupled finite element.
Various parametric studies are performed to analyze the effect of radius of curvature, layer thickness, mechanical and electrical size effects on the response of the flexoelectric beam. The contribution of converse flexoelectric effect in the beam response is compared with that of the direct effect. The converse effect gives rise to a larger flexoelectric response, which is also noted in recent experimental studies [1]. The effective piezoelectric coupling coefficient, calculated in terms of the voltage generated across the thickness of flexoelectric material, is compared with the coupling coefficient of piezoelectric materials. Our simulations show that this value increases drastically for smaller thickness of flexoelectric materials and surpasses the coupling coefficient of piezoelectric materials for very thin structures (thickness of the order of a few micron). Along with the mechanical size effects, our model is able to capture electrical size effects at lower scales, as observed by Ren and Sun [19]. The finite element framework developed in this work is generalized and its application to sensor and actuator device design has been demonstrated through various examples. Y. S. Joshan and S. Santapuri ZAMP

Mathematical formulation
In this section, the development of a two-way coupled linear flexoelectric theory, incorporating both direct and converse flexoelectric effects, is presented. A majority of existing literature deals with analysis of direct flexoelectric effect (i.e., presence of strain gradient giving rise to electric field) and its applications [3,16,18,43]. More recent literature has focused on converse flexoelectric effect in these materials where mechanical stress is induced due to electric polarization gradient [34,36,37]. This effect is found to be predominant at nanoscale and often gives rise to a large electric response even in non-piezoelectric materials. This has also resulted in incorrect piezoelectric coefficient calculation using piezoresponse force microscopy at nanoscale [1]. Motivated by these recent findings, in this work we develop a comprehensive flexoelectric theory which incorporates a two-way coupling of electroelastic fields and their gradients. This generalized framework is implemented towards modeling of composite flexoelectric curved beams. An electromechanically coupled C 2 continuous finite element, based on classical beam theory, is developed to analyze the flexoelectric curved beams. The strain gradient and electric filed gradients are larger in curved beams as compared to straight beams due to curvature effects and it enhances direct and converse flexoelectric effects in curved structures. Simulations describing both the sensing and actuating modes of a curved beam are presented to demonstrate its possible applications. A 1D Fourier series analytical solution is also developed for a simply supported beam, which is used to verify the accuracy of our finite element results.
We note that a majority of the existing flexoelectric theories in literature are developed in terms of polarization and polarization gradient as independent variables [16,44]. However, the present formulation is developed considering electric field and its gradient as independent variables. Electric displacement vector and a higher-order electric displacement tensor are the corresponding dependent variables. It is convenient to obtain the governing equations and the corresponding weak form by representing electric field as gradient of electrostatic potential, thus reducing the number of electrical dependent variables. The electrical boundary conditions, i.e., ground, open/closed circuit conditions used in experimental investigations can be easily defined in terms of electrostatic potential and charge displacement as compared to polarization.

3D governing equations and boundary conditions
The equilibrium governing equations and boundary conditions describing a flexoelectric solid are obtained using the principle of minimum potential energy. For a general electroelastic system, the total potential energy Π is defined as where V 0 represents the total volume occupied by the flexoelectric solid, dV is the infinitesimal volume element, ω denotes the work done by the externally applied mechanical and electric loads. The free energy function Ω is expressed in terms of the infinitesimal strain tensor S ij , strain gradient G ijk , electric field E i , and the electric field gradient K ij . The infinitesimal strain tensor and its gradient are defined as where u i represent the components of the mechanical displacement vector. Similarly, the electric field vector E i , and electric field gradient tensor K ij are given by where φ represents the scalar electrostatic potential. It can be observed that the electric field gradient tensor K ij is symmetric and the strain gradient tensor obeys G ijk = G jik , owing to the symmetry of strain tensor. The governing equations are obtained by setting the first variation of potential energy (1) to zero, i.e., where δΩ denotes the first variation of the free energy function given by where Here, T ij and H ijk represent the Cartesian components of the Cauchy stress tensor and the higher-order stress tensor, respectively. D i denotes the components of electric displacement vector and B ij denotes the higher-order electric displacement tensor components. The variation of the free energy function δΩ can be expressed in terms of strain and electric potential variables using (2)-(3) as whereĤ ijk = 1 2 (H ijk + H jik ). Also, the first variation of the work done by external electroelastic loads, δω, is defined as where ∂V t 0 represents a subset of the boundary ∂V 0 on which the applied traction t a i is prescribed. Similarly, ∂V e 0 represents the subset on which the surface charge density ρ s is prescribed. Also, f b i represents the components of the body force and ρ v denotes the applied volumetric charge density.
The first term in Eq. (4) is obtained by integrating (7) across the material volume as follows: where each term is rewritten using integration by parts as Y. S. Joshan and S. Santapuri ZAMP Substituting these expressions in (9) and applying the Gauss's divergence theorem, we get where n i denotes the unit normal to the boundary surface ∂V 0 . The terms in boundary integral involving gradients of δu i and δφ are further simplified by decomposing the gradient operator as follows: where Δ j (·) = (δ jk − n j n k )(·) ,k and (·) = n k (·) ,k , denote the surface gradient and the normal gradient (derivative along the surface normal) operators, respectively. Now, the boundary integral term ∂V0Ĥ ijk n j δu i,k dA is decomposed in the following manner: The surface ∂V 0 of the electroelastic solid can be assumed to be divided into finite number of surfaces S m , each surface with edge boundary c m . Using the Stoke's surface divergence theorem each surface S m [42], Eq. (17) is further simplified as follows: where n k is the surface unit vector normal to the edge c m , m s represents the total number of sharp edges of the body, and [[·]] denotes the jump in the value of the enclosed quantity across the sharp edge. Following the same procedure, the boundary integral term Substituting Eqs. (18)- (19) in Eq. (14), we get Finally, the governing equations and boundary conditions are obtained by substituting (20) and (8) in (4) by setting the coefficients of δφ and δu i zero in the volumetric and surface terms. The 3D governing equations, applicable at every point within V 0 , are given by The boundary conditions, specified at any point on ∂V 0 , maybe prescribed in terms of applied traction/charges or displacement/electric potential and their derivatives as follows: Finally, the edge boundary conditions, prescribed on the sharp edges of the boundary, are given by [[n jĤijk n k ]] = 0 or u i = 0.

Linear flexoelectric constitutive modeling
Flexoelectric effect, exhibited by dielectric materials at micro/nanoscales, describes the coupling of strain gradient and polarization. Conversely, flexoelectric materials also exhibit a coupling of electric field gradient and strain. This effect is significant at lower dimensions as the strain and electric field gradient terms scale up at micro/nanoscales [44]. The electroelastic free energy function Ω corresponding to a linear flexoelectric material is given by where S ij and G ijk denote the Cartesian components of the infinitesimal strain tensor and the strain gradient tensor, respectively. Similarly, E i represents the electric field vector components, and K ij represents the electric field gradient. The components c ijkl represent the stiffness tensor, g ijklmn is the higher-order stiffness tensor, ij is the electric permittivity tensor, k ijkl represents the higher-order electric permittivity tensor, e ijk is the piezoelectric coupling tensor, f ijkl is the direct flexoelectric coupling tensor and h ijkl is the converse flexoelectric coupling tensor. The corresponding linear constitutive equations are given by where T ij and H ijk denote the Cartesian components of the Cauchy stress tensor and the higher-order stress tensor, respectively. Also, D i represents the components of electric displacement vector and B ij denotes the higher-order electric displacement tensor components.

Specialization to isotropic materials.
The linear constitutive Eqs. (30)- (33) are now specialized to isotropic materials by specializing the material constants c ijkl , g ijklmn , ij , k ijkl , f ijkl , and h ijkl in their isotropic forms [45] as follows: wherec i (i = 1, 2, 3),g i (i = 1, 2, ...15), ,k i (i = 1, 2, 3),f i (i = 1, 2, 3) andh i (i = 1, 2, 3) are the scalar material constants corresponding to an isotropic material (to be determined through material characterization experiments), and δ il represents the second-order isotropic Kronecker delta function. The piezoelectric coupling tensor e ijk vanishes for isotropic materials and any material with a centrosymmetric crystal structure [46]. The generalized constitutive equations defined in (30)- (33) are now specialized to isotropic materials, i.e., where λ =c 1 , 2μ =c 2 +c 3 , 2a 1 =g 3 , 2a 2 =g 1 +g 2 +g 6 +g 15 , 2a 3 =g 4 +g 5 +g 7 +g 11 , 2a 4 =g 9 +g 12 , 2a 5 =g 8 +g 10 +g 13 +g 14 , The expression for free energy function (29) is reduced to ZAMP A gradient electromechanical theory Page 9 of 33 178 2.2.2. Representation of gradient effects in terms of length scale parameters. It has been shown in Lam et al. [42] that in a linear elastic isotropic material, the dependence on stain gradient reduces to the following quantities: where γ i represents the dilatation gradient vector in component form, S mm is the trace of the infinitesimal strain tensor, η 1 ijk is the deviatoric stretch tensor, χ s ij denotes the symmetric part of curvature tensor, η ijk = u k,ij is the third-order deformation gradient, and η s ijk = In a strain gradient theory, the material constants a i (i = 1, 2, ..5) describing the relationship between higher-order stress and strain gradient (Eq. (41)), are expressed as functions of characteristic length scale parameters. These length scale parameters are phenomenological material constants that describe the size effects observed in micro/nano-structures. In this work, we use three length scale parameters to explain the size effects. The higher-order elastic constants a i can be expressed in terms of the three length scale parameters l 0 , l 1 and l 2 as follows: where μ represents the elastic Lame's constant appearing in the constitutive Eq. (41). These relationships were derived in Lam et al. [42] using the reduced form of constitutive equations for an linear elastic material (no electric field coupling), i.e., where p i , τ 1 ijk and m s ij are the work conjugate of γ i , η 1 ijk and χ s ij , respectively. We note that there are several possible formulations available in literature [22,27,[47][48][49][50][51].
For a dielectric material, the constitutive equation relating higher-order electric displacement tensor B ij and electric field gradient K ij (neglecting the flexoelectric effect) is given by wherein the electric field gradient tensor is decomposed in two independent components corresponding to higher-order dielectric constants k 1 and k 2 . In this work, we additionally define electrical length scale parameters l e 0 and l e 1 to describe the two higher-order dielectric constants k 1 and k 2 appearing to consider electric size effects at micro/nanoscales as follows: where represents the dielectric permittivity of the material. Symmetry of the electric field gradient tensor allows the representation of the higher-order dielectric coefficients using only two length scale parameters.

Boundary value problem: a flexoelectric curved beam
In what follows, the generalized 3D formulation is specialized to curved flexoelectric beams operating in sensing and actuation modes. In the present formulation, a composite curved beam of radius of curvature Y. S. Joshan and S. Santapuri ZAMP R and thickness h is analyzed in the curvilinear coordinate system (α 1 , α 2 , α 3 ), as shown in Fig. 1. The coordinate α 3 is along the normal to the neutral surface (α 3 = 0), α 2 is along the width of the beam, and α 1 denotes the arc length along the 2D projection of the neutral surface (dotted line shown in Fig. 1). The arc length of the beam along the neutral surface is taken as L and its width is taken as b. The composite beam consisting of an isotropic flexoelectric layer at the top (active layer) and a passive substrate layer at the bottom is considered. The thickness of the flexoelectric layer is denoted by h f .

Displacement field
In what follows, modeling of the curved flexoelectric beam is presented using the classical beam theory (CBT) assumptions [52][53][54]. The displacement field, defined at any point (α 1 , α 2 , α 3 ) on the beam, is given by where u 1 , u 2 and u 3 represent the components of displacement vector; u 0 and w 0 are the components of displacement of any point of the α 1 axis along α 1 and α 3 directions, respectively, on the centroidal axis (along α 1 ) of the beam. Following Kirchhoff hypothesis, i.e., normal to the neutral surface remain straight and normal during deformation, the only nonzero strain component S 11 is obtained as follows [55]: The corresponding nonzero strain gradient components are given by

Variation of the electric potential across the thickness of the flexoelectric beam
The electrostatic potential is often assumed to be linearly varying across the thickness of the flexoelectric beam [44,56,57]. However, higher-order correction terms need to incorporated for thicker beams [58]. These higher-order terms are modeled as quadratic or other non-polynomial functions in literature [59,60]. In the present formulation, a quadratic variation of electrostatic potential φ is assumed as follows: where φ a is the linearly varying potential term applied at the top of the flexoelectric layer, φ 1 is the higher-order correction term, α u and α l are the thickness coordinates corresponding to the top and bottom surfaces of the flexoelectric layer. The coefficient g 1 (α 3 ) is a quadratic interpolation function defined by Similarly, g a (α 3 ) is a linear interpolation function given by The electric field components are obtained using E i = −Grad φ as follows: where (·) denotes derivative of the quantity along the normal, i.e., along α 3 . The nonzero components of the electric field gradient tensor are given by

Two-way coupled flexoelectric beam formulation
The equations governing a curved flexoelectric beam under transverse electromechanical loading are derived using the principle of virtual work defined in Eq. (4). For a composite curved beam, the first variation of work done due to externally applied electromechanical loads is given by where t a 3 represents the resultant traction along the transverse direction on the top and bottom surfaces (α 3 = ±h/2), t a 1 represents the axial traction on the cross sections α 1 = 0, L, and ρ s denotes the surface charge density applied at the top and bottom surfaces (α 3 = ±h/2) of the flexoelectric layer. The volumetric charge density and mechanical body force terms have not been included in our analysis. However, Eq. (64) maybe easily modified to incorporate these terms. For thin beam structures, the traction components are expressed in terms of applied transverse and axial load as follows: where q + 3 and q − 3 are transverse loads per unit length acting on the top and bottom surfaces of the curved beam, q 3 is the resultant transverse load per unit length of the beam, q 0 1 is the axial load per unit cross-section area, and q 1 is the resultant axial load acting at the ends of the beam, i.e., at α 1 = 0 and α 1 = L. Similarly, the applied surface charge density on the flexoelectric layer can be expressed as where ρ L is the applied charge per unit length. Using (65)-(67), Eq. (64) is reduced to Finally, the energy function (44) is reduced to where the nonzero strain and strain gradient components S 11 , G 111 and G 113 are given by (53)-(55); the electric field and electric field gradient components E 1 , E 3 , K 11 , K 33 and K 13 are given by (59) where the axial force resultant N 11 , moment resultant M 11 , higher-order stress resultants N h 111 , N h 113 , higher-order moment resultant M h 111 , charge displacement resultants D g , D g , D a , D a and higher-order charge displacement resultants B g , B g , B g , B a , B a , B a are given by ⎡ Applying integration by parts on the volumetric terms, Eq. (70) is further reduced to The terms corresponding to variation of independent variables (δu 0 , δw 0 , δφ 1 , δφ a ) are grouped and the governing equations are obtained using the fundamental lemma of variation [61]. As δu 0 , δw 0 , δφ 1 and δφ a are arbitrary, the governing equations can be obtained equating their coefficients in the energy expression to zero. The Euler-Lagrange equations for the curved flexoelectric beam are obtained by substituting (68) and (77) in Eq. (4), and the resulting governing equations are given by and the corresponding natural and essential boundary conditions, defined on the two ends of the curved beam α 1 = 0, L, are given by

Finite element framework
In what follows, an electromechanical finite element framework is developed to solve the curved flexoelectric beam problem for different boundary and loading conditions. The principle of virtual work (4) is used to obtain the finite element equations governing the curved flexoelectric composite beam subjected to electric and mechanical loads.

Shape functions.
The flexoelectric beam is discretized into smaller two-noded curved beam elements, each of length L e . Variation of the electroelastic field variables (u 0 , w 0 , φ 1 , φ a ) are expressed in terms of shape functions within each element as follows: where U = [u [1] 0 u [1] 0,1 u [2] 0 u [2] 0,1 ] T , W = [w [1] 0 w [1] 0,1 w [1] 0,11 w [2] 0 w [2] 0,1 w 1,1 ] T and Φ a =[φ [1] a φ [1] a,1 φ [2] a φ [2] a,1 ] T denote the degrees of freedom within each element. The superscript (·) [i] denotes the nodal value of the variable at the ith node of the beam element, and the subscript (·) ,1 denotes the differentiation of the quantity with respect to α 1 . The corresponding shape functions, i.e., N=[N 1 N 2 N 3 N 4 ] such that N i denote C 1 continuous shape functions of cubic order and P=[P 1 P 2 P 3 P 4 P 5 P 6 ] is composed of C 2 continuous shape functions of quintic order. The polynomial shape functions N and P are ZAMP A gradient electromechanical theory Page 15 of 33 178 derived following the standard finite element procedure [62] by taking the continuity requirements of each variable into consideration.

Derivation of finite element governing equations.
Variation of the elemental free energy Ω (l) corresponding to the lth element is obtained using Eq. (70). The field variables (u 0 , w 0 , φ 1 , φ a ) are subsequently expressed in terms of shape functions and elemental degrees of freedom defined in Eq. (83). Performing the integration along the axis (i.e., with respect to α 1 ) within each element, we obtain the elemental free energy in the form where the matrices K where The variation in the total free energy δΩ for the flexoelectric curved beam is obtained by assembling the elemental energies of all the elements to obtain global stiffness matrix given by where K is the global stiffness matrix, d denotes the global field variable vector consisting of the degrees of freedom values at every node, and n e is the total number of elements used for discretization of the beam. Further, elemental work done due to the external load is obtained in terms of the elemental force vector as follows: where F (l) ge denotes the elemental force vector. Again, performing assembly over all the elements, the variation in the total external work done is evaluated in terms of the global force vector as Y. S. Joshan and S. Santapuri ZAMP where P F represents the global force vector, which includes both mechanical and electrical load terms. Finally, using the principle of energy minimization, the finite element governing equation is obtained as follows: Given that δd T can be arbitrarily assigned, (89) is satisfied for all δd T iff These global finite element governing equations are algebraically solved for given boundary conditions to obtain nodal values of mechanical displacements and electrostatic potential. This formulation is used to analyze the flexoelectric curved beam, both in sensor and actuator modes. In sensor mode, electrostatic potential generated across the flexoelectric layer is examined for applied mechanical loads. In actuation mode, the deformation of the curved beam is studied for applied electrostatic loads. An analytical solution is also developed using the Fourier series method for simply supported boundary conditions to verify the accuracy of the developed finite element formulation. The details of the procedure for analytical solution are given in "Appendix B".

Results and discussion
The finite element framework is implemented to analyze the sensor and actuator response of curved flexoelectric composite beams under different loading and boundary conditions. A MATLAB program is written to compute the global stiffness matrix and force vector. The following mechanical boundary conditions are implemented at the ends of the flexoelectric curved beam: Simply supported : The electrical boundary conditions are given by Charge boundary conditions (CBC) : Potential boundary conditions (PBC) : Higher − order potential boundary conditions (HPBC) : boundary conditions, loading conditions, curvature, thickness of the flexoelectric layer etc. Specifically, microscale effects, described in terms of mechanical and electrical length scale parameters, are studied for flexoelectric beams operating in actuator and sensor modes. The effective piezoelectric coefficient is calculated for the curved flexoelectric beam and compared with piezoelectric materials at different length scales. It is observed that both direct and converse flexoelectric effects are prominent at small scales, which makes flexoelectric materials suitable for applications in MEMS and NEMS.

Bending analysis of a passive isotropic micro-beam
A passive cantilever micro-beam is analyzed under a point load (q 0 = 100 µN) applied at the free end (α 1 = L). The material properties of the micro-beam, made of epoxy, are taken as [42,63] Transverse deflection w (in µm) of the micro-beam is calculated using (i) classical theory i.e. l 0 = l 1 = l 2 = 0, (ii) modified couple stress theory i.e. l 0 = l 1 = 0, l 2 =17.6 µm, and (iii) strain gradient theory i.e. l 0 = l 1 = l 2 = 17.6 µm [42]. Variation of transverse deflection w along the neutral axis of the beam is presented in Fig. 2. Our finite element results agree well with analytical results presented in Kong et al. [63].

Bending analysis of flexoelectric curved beams under applied mechanical load: sensor mode
In this section, a unimorph beam consisting of a flexoelectric layer at the top and a passive substrate layer at the bottom (shown in Fig. 1), is analyzed for different boundary and loading conditions. The following abbreviations are used for different types of mechanical loads: 1. SSL: Sinusoidal load q 3 = q 0 sin(πα 1 /L), 2. UDL: Uniformly distributed load q 3 = q 0 , 3. HSL: Hydrostatic load q 3 = q 0 (α 1 /L). Y. S. Joshan and S. Santapuri ZAMP Here, q 0 represents the magnitude of the applied mechanical load. Material properties for the flexoelectric and passive layers are given in terms of Young's modulus (c), Poisson's ratio (v), flexoelectric coefficients (f 1 , f 2 ), and electric permittivity ( 1 ) as follows: • Flexoelectric material [64]: The converse flexoelectric constants (h 1 , h 2 ) and direct flexoelectric constants (f 1 , f 2 ) are considered to be equal in the present analysis [65]. The value of mechanical length scale parameters (l 0 = l 1 = l 2 ) are taken as 2 µm [64]. Ren and Sun [19] observed an enhanced electric permittivity of the material at microscale (similar to mechanical stiffness). The electric length scale parameters (l e 0 , l e 1 ) are thus defined to incorporate lower scale effects through a higher-order permittivity term and their values are chosen phenomenologically to accurately predict the flexoelectric material behavior. In this work, the values of electric length scale parameters are considered to be equal to the mechanical length scale parameters, i.e., l e 0 = l e 1 =2 µm [33]. The sensor response is now analyzed in terms of the nondimensional potential Φ defined as where c 0 =1 GPa is the scale for Young's modulus of the material. The flexoelectric beam is analyzed under different loading conditions (SSL, UDL, HSL). The beam is assumed to be simply supported on both ends and the potential difference across the flexoelectric layer is calculated. Potential boundary conditions, defined in Eqs. (92) and (95), are thus imposed. The FEM results are compared with analytical results in Tables 1 for different L/R ratios. The FEM results agree well with analytical results and a maximum error of 0.0221% is noted. The nondimensional potential Φ(L/2, h/2) is maximum for UDL for all L/R ratios. Moreover, Φ increases with increase in the L/R ratio of the flexoelectric curved beam due to curvature effects. Sensor response of a curved flexoelectric beam is now analyzed for different boundary conditions. Specifically, simulations are performed for three sets of boundary conditions: (i) Clamped end (91) with charge boundary condition (94), (ii) Simply supported edge (92) with potential boundary condition (95), and (iii) Free end with higher order potential boundary condition (96). Variation of the nondimensional potential Φ along the span α 1 is presented in Fig. 3 for the three boundary conditions specified above. For a simply supported beam, Φ is maximum at the centre due to symmetric loading and geometry, while it is maximum at the clamped end for other boundary conditions. In addition, the maximum value of Φ is obtained for cantilever boundary conditions among the considered cases. The cantilever flexoelectric beam is thus a suitable geometry for the design of micro sensors.
It may be noted that the potential generated across the flexoelectric layer, due to an applied mechanical load, has contributions from both direct and converse flexoelectric effects. Flexoelectric response due to the converse effect is compared with the response due to the direct effect for different thickness scale values of a cantilever flexoelectric sensor subjected to a uniform distributed load. The results are presented in Fig. 4 for two cases: (i) considering both direct and converse flexoelectric effects (h 1 = h 2 = f 1 = f 2 ), and (ii) considering only direct flexoelectric effect and neglecting converse flexoelectric effect (h 1 = h 2 = 0). It is observed that the converse flexoelectric effect results in a larger response. Thus, it may not be suitable to ignore these effects in all flexoelectric material-based applications.
In what follows, effective piezoelectric coefficient of the flexoelectric material is calculated using the formula (d 31 ) ef f = − 1 E 3 (0, h/2)/(T 11 (0, h/2)) for open circuit conditions at the clamped edge of the beam. The values are compared with the piezoelectric coefficient of Quartz (a commonly used piezoelectric material) for different thickness values in Fig. 5. It is observed that the effective piezoelectric coefficient increases significantly for smaller thickness values. Whereas, the value of d 31 should be independent of its structural thickness [1]. It has been observed by Abdollahi et al. [1] that at smaller thicknesses (of the order of a few micron), atomic force microscopy measurements in piezoelectric materials resulted in ZAMP A gradient electromechanical theory Page 19 of 33 178 Table 1. Now, the effect of electrical length scale parameters (l e 0 , l e 1 ) on the sensor response of the cantilever flexoelectric beam is analyzed. It is subjected to different loading conditions (UDL, SSL, HSL) and the nondimensional potential Φ is evaluated at the clamped end. The size effects on the permittivity of the material is analyzed in terms of the l/h(l = l e 0 = l e 1 ) ratio and the results are presented in Fig. 6. It is observed that an increase in the l/h ratio results in a lower Φ. This can be attributed to an increase in the effective permittivity of the material at smaller scales, which was also observed in recent studies by Ren and Sun [19].  Finally, the effect of flexoelectric layer thickness h f on the output potential Φ(0, h/2) is analyzed. Variation of nondimensional potential Φ with h f /h is presented in Fig. 7 for a cantilever flexoelectric beam subjected to UDL for different L/R ratios. The maximum potential is observed for h f /h = 1. Hence, a flexoelectric beam with no substrate layer is an optimum choice or sensor applications. In case Y. S. Joshan and S. Santapuri ZAMP  of piezoelectric materials, in which the electric field is coupled with bending stresses, maximum potential is obtained for a piezoelectric layer thickness ratio in the range of 0.7-0.8 [66]. In case of flexoelectric beams, the strain gradient G 113 is constant across the thickness of the beam. Increase in thickness of the flexoelectric layer results in larger strain gradients, thus increasing the electrostatic potential generated across the flexoelectric layer.

Bending analysis of flexoelectric curved beams under applied electrical load: actuator mode
We now analyze the curved flexoelectric beam in actuator mode wherein the deformation of the beam is analyzed under a uniform electrical potential (ΔΦ = 100 V) applied across the thickness of the flexoelectric layer. Variation of nondimensional deflection w/h along the span length α 1 /L of the beam is presented in Fig. 8 for different L/R ratios. The FEM results are compared with analytical solution and a good agreement is observed between them. It is observed that increase in the L/R ratio of the beam leads to an increase in the actuator deflection. Further, the cantilever beam actuator is analyzed under different electromechanical loads. Variation of nondimensional deflection w/h along the span of the flexoelectric beam is presented in Fig. 9 for two cases: (i) varying electrical load for a constant mechanical load, and (ii) varying mechanical load for a constant electrical load. As expected, mechanical load induced deformation is dominant compared to electric potential induced deformation.
Next, actuator response of the flexoelectric beam is analyzed for different thickness values. The different values of beam thickness h are taken as 2 µm, 4 µm, 20 µm, and 80 µm. Variation of nondimensional deflection w/h along the span of the flexoelectric beam is presented for classical theory (l 0 = l 1 = l 2 = 0, l e 0 = l e 1 = 2 µm), modified couple stress theory (l 0 = l 1 = 0, l 2 = 2 µm, l e 0 = l e 1 = 2 µm) and strain gradient theory (l 0 = l 1 = l 2 = 2 µm, l e 0 = l e 1 = 2 µm) in Fig. 10. A considerable difference is noted among the results of these theories for lower thickness values i.e., h = 2 µm and 4 µm. There is a negligible difference between the results of the three theories for h = 80 µm, as size effects become negligible with increasing thickness.
We now study the effect of the mechanical length scale parameters on the actuator response of a cantilever flexoelectric beam. The flexoelectric beam is subjected to uniform electrostatic potential load (ΔΦ=100 V). The maximum deflection at the free end of the cantilever flexoelectric beam is calculated. The effect on the actuator response is analyzed in terms of length scale to thickness ratio l/h and the individual effects of three mechanical length scale parameters (l 0 , l 1 and l 2 ) are presented in Fig. 11. The increase in the l/h ratio increases the effective stiffness of the flexoelectric beam actuator, which leads to a decrease in the nondimensional deflection w/h. The size effects are prominently influenced by Y. S. Joshan and S. Santapuri length scale parameter l 0 as compared to l 1 and l 2 . In this study, the values of length scale parameters are considered to be equal. However, different values of these parameters can be chosen to better fit the experimental results. The strain gradient theory gives the choice of three length scale parameters (l 0 , l 1 , l 2 ), which makes it more versatile compared to modified couple stress theory, which considers only one length scale parameter l 2 . As the last example, the effect of flexoelectric layer thickness h f is analyzed on the actuator response of a cantilever flexoelectric beam actuator in terms of thickness ratio (h f /h). In this study, two different types of passive substrates are considered. The material properties of substrate 1 (Aluminium) are taken as c p = 70 GPa, v p =0.3, while for substrate 2 (Steel), the material properties are taken as c p = 220 GPa, v p =0.29. A uniform electrostatic potential difference (ΔΦ=100 V) is applied across the flexoelectric layer. The variation of the nondimensional deflection w/h with h f /h for L/R=0 and 0.5 is presented in Fig. 12. It is noted that the flexoelectric beam with substrate 1 shows a decrease in the nondimensional deflection with increase in the layer thickness ratio. However, the nondimensional deflection increases for the flexoelectric beam with substrate 2 with increase in the thickness ratio. This is due to the fact that the stiffness of substrate 1 made of Aluminium is lesser than that of flexoelectric material. With increase in the thickness ratio, the flexural rigidity of the beam increases which leads to the decrease in the deflection. On the contrary, the flexural rigidity of the beam with substrate 2 decreases with an increase in the thickness ratio, which increase the actuator deflection.

Conclusions
A generalized gradient electromechanical theory is developed for flexoelectric curved beams considering both direct and converse flexoelectric effects. The model is specialized to isotropic materials and the mechanical and electrical size effects are considered using phenomenological length scale parameters. A novel C 2 continuous finite element framework is developed to solve the beam governing equations for different loading and boundary conditions. An analytical solution is also derived to verify the accuracy of the finite element results. Specifically, the sensor and actuator response of flexoelectric curved beams is analyzed for different scales by varying the beam thickness. Our results show that the electromechanical coupling coefficient for flexoelectric materials is more than standard piezoelectric materials for lower values of thickness. For instance, coupling coefficient calculated using d 31 = − 1 E 3 /T 11 , is larger in a non-piezoelectric material exhibiting flexoelectricity than that of Quartz at small thicknesses (h ≤ 8 µm).
It is further observed that the converse flexoelectric effect gives rise to a large flexoelectric response at microscale. Finally, the effect of mechanical and electric length scale parameters on the response of flexoelectric curved beams is analyzed. The length scale parameters account for an increase in the material stiffness and dielectric permittivity at micro/nanoscales. The actuator response of flexoelectric composite beam depends on the stiffness of the passive substrate. A substrate with less stiffness can be preferred Y. S. Joshan and S. Santapuri ZAMP . Substrate 1 is taken as Aluminum and substrate 2 is taken as Steel. The variation of nondimensional deflection with thickness ratio is dependent on the substrate stiffness to obtain higher actuator deflection. In case of the flexoelectric sensors, the maximum voltage output is obtained for flexoelectric layer thickness ratio h f /h = 1 (i.e., no substrate layer). Our results also agree with the trends observed in recent experimental works.
Here, the superscript (l) has been dropped for notational convenience.
[E ca Y. S. Joshan and S. Santapuri ZAMP field gradient components (61)- (63). Further, the explicit form of the governing equations are solved for simply supported conditions on the two ends of the beam. The Fourier series solution is used to express the primary variables in terms of series solution as follows: where p = πm/L. The coefficients U m , W m , Φ m and γ m are unknowns to be solved using the Fourier series solution technique. In addition, transverse load q(α 1 ) and surface charge density ρ s (α 1 ) are expressed in terms of Fourier series as where q m = 2 L L 0 q 3 (α 1 ) sin(pα 1 ) dα 1 , Π m = 2 L L 0 ρ L (α 1 ) sin(pα 1 ) dα 1 .
The series solutions (132)-(137) are substituted into the explicit form of the governing equations, resulting in algebraic equations of the form where U g is the vector of unknowns to be solved, which is given by F g is the resultant force matrix, and R g is the 4 × 4 resultant stiffness matrix. Equation (140) is solved to obtain the displacements and electrostatic potential in terms of Fourier series coefficients.