Computational Modeling of Flexoelectricity—A Review

Electromechanical coupling devices have been playing an indispensable role in modern engineering. Particularly, flexoelectricity, an electromechanical coupling effect that involves strain gradients, has shown promising potential for future miniaturized electromechanical coupling devices. Therefore, simulation of flexoelectricity is necessary and inevitable. In this paper, we provide an overview of numerical procedures on modeling flexoelectricity. Specifically, we summarize a generalized formulation including the electrostatic stress tensor, which can be simplified to retrieve other formulations from the literature. We further show the weak and discretization forms of the boundary value problem for different numerical methods, including isogeometric analysis and mixed FEM. Several benchmark problems are presented to demonstrate the numerical implementation. The source code for the implementation can be utilized to analyze and develop more complex flexoelectric nano-devices.


Introduction
In recent years, huge attention in various research areas, such as material science, mechanical engineering and chemistry, has been paid to flexoelectricity [1,2], an electromechanical coupling phenomenon in all dielectric materials (regardless the material symmetry). It describes the induced electric polarization due to strain gradients. Despite the small coefficients as compared to its counterpart, piezoelectricity, flexoelectric effect has gained more attention and is appreciable thanks to the miniaturized trend, where many electromechanical systems are operated at submicron or nano length scale, leading a tremendous increase in the strain gradient. Due to size effects from strain gradients, flexoelectricity holds promise in micro/nano-electromechanical systems such as sensors [3][4][5], actuators [6,7] and nano-energy harvesters [8][9][10][11][12][13]. It can also explain physical phenomena such as strain gradient-driven polarity control [14][15][16], flexo-photovoltaic effects [17] and flexo-caloric effects [18]. Interested readers are referred to other excellent review articles about flexoelectricity for more details [19][20][21][22].
To our best knowledge, the first continuum theory for flexoelectricity in dielectric solid was proposed by Maranganti et al. [23], who considered both direct and converse flexoelectric effect in the internal energy that takes strain, strain gradient, polarization and polarization gradient as independent variables. As noted in their work, the formulation can be considered as a complement for the polarization-gradient theory proposed by Mindlin [24] and later on by Sahin and Host [25]. Based on their formulation, Sharma's group further studied the enhancement of flexoelectricity in nano-structures [26][27][28][29][30]. It should be remarked that, in the above-mentioned formulation, the electric force, which is the divergence of the Maxwell stress, was not taken into account. The consideration of Maxwell stress may be important in nano-dielectric solid, as pointed out by Hu and Shen [31], who utilized the physical variation principle proposed by Kuang [32][33][34][35] so that the electrostatic stress appears naturally from the minimization of the electric enthalpy. The significance of the Maxwell stress was also discussed in the unified energy formulation of continuum magneto-electro-elasticity in the work of Liu [36,37].
Nonetheless, by virtue of the non-locality, i.e., strain gradient or polarization gradient, flexoelectricity in a solid dielectric is governed by a fourth-order partial differential equation. Therefore, modeling flexoelectricity is mostly simplified to a 1D model or problems with simple geometries with axisymmetry (such as a cylinder or a disk) or plate. As in strain gradient elasticity, the difficulty in fully multi-dimensional modeling for flexoelectricity arises from the C 1 continuity in the approximation of the displacement field. This issue was solved firstly by Abdollahi et al. [38] using the local maximum-entropy (LME) meshfree method. Analogous to strain gradient elasticity, mixed FEM was also adopted to solve flexoelectric problems by Mao et al. [39,40]. Isogeometric analysis (IGA) is also a robust approach to flexoelectricity, as in [41][42][43][44], thanks to its straightforward handling of higher-order continuity. Alternatively, the Argyris triangular element can also be used as a remedy to the C 1 requirement, as in the work of Yvonnet and Liu [45]. More recently, in a similar spirit to IGA approach, Codony et al. [46] proposed the immersed boundary method on hierarchical B-Spline for flexoelectricity problems with higher efficiency. Besides, flexoelectric effect and the associated damage behavior were also approached with peridynamics method by Roy and Roy [47].
As the research field of flexoelectricity is rapidly developing, the need of numerical simulations is apparently inevitable. However, there has not been a review on the computational aspects of flexoelectricity yet. Therefore, in this paper, we aim to provide an overview as well as detailed procedure on modeling flexoelectricity. Different numerical methods, including isogeometric analysis, mixed FEM and meshfree, are employed to solve the boundary value problem. The remainder of the paper is structured as follows. In Section 2, motivated by the existence of the Maxwell stress tensor in dielectric solid, we summarize the physical variational principle in which the electrostatic stress tensor appears naturally as the minimization of the electric enthalpy from the work of Hu and Shen [31]. Furthermore, we present an alternative formulation, derived from the electric Gibbs free energy (but with different independent variables), which also results in the electrostatic stress tensor. Noticeably, the two formulations can be considered as a generalized formulation, for which, when the electrostatic stress is omitted, the strong forms from previous literature can be retrieved. Consequently, we obtain the weak forms from Galerkin method for different numerical approaches. Details on the numerical implementation of the discretization forms can be found in Section 3 together with source code wherever possible. Several numerical examples to validate as well as illustrate the numerical implementation are presented in Section 4.

Physical Variational Principle
Hu and Shen [31] extended the variational principle from Toupin [48] by the physical variational principle proposed by Kuang [32][33][34][35] to derive the governing equations for flexoelectricity, in which the Maxwell's stress is taken into account. We briefly summarize the theory and formulation in this section. Assuming the internal energy as a function of strain ij , strain gradient ij,k , polarization P i and polarization gradient P i,j , U ij , ij,k , P i , P i,j = 1 2 a ij P i P j + 1 2 b ijkl P i,j P k,l + 1 2 c ijkl ε ij ε kl + d ijk P i ε jk + f ijkl ε ij P k,l + h ijkl P i ε jk,l + 1 2 g ijklmn ε ij,k ε lm,n , (1) from which the constitutive relations are obtained as where σ ij , E i , τ ijk and V ij are the (local) stress, electric field, higher-order stress and higher-order electric field, respectively, while a ij , b ijkl , c ijkl , d ijk , f ijkl , and h ijkl are the material tensors [32]. In other words, the variation of internal energy can be written as The energy density is the sum of the internal energy and the internal of the Maxwell self-field [24] in which the electric Maxwell self-field can be defined as the gradient of the electric potential φ Next, we can define the electric enthalpy as We consider a dielectric solid occupying a volume v bounded by a surface a that separates from the surrounding environment v (considered to be air in this study). Taking v * = v + v , the equilibrium of such system in the the isothermal process can be defined by the variational principle where W is the external work done by body force f i , external electric field E ext i and free charge ρ f on the solid body and tractiont i , double-tractionτ i , surface charge σ * and higher-order electric fieldv i on the surface boundary. Hence, its variation δW is defined as (8) in which the operatorD n is explained in the following section.
where H e = 1 2 E i P i + V ij P i,j is the energy deducted from the total energy minus the pure deformation energy. Moreover, the migratory variation is the essential feature of physical variational principle, in the sense that the virtual displacement not only produces the strain variation but also causes electric potential and polarization variations, such as [31] where δ φ (•) and δ P (•) are the variations produced by the (local) variation of the electric potential and polarization, respectively, whereas δ u (•) is the migratory variation caused by the virtual displacement. Subsequently, by substituting Equation (10) into Equation (9) with the use of Equations (6) and (3) and employing Gauss divergence theorem, through a lengthy but straightforward manipulation, the variation of the electric enthalpy can be re-expressed as [31] δ where n i is the normal vector.D t k (•) = (δ kl − n k n l ) ∂ ,l (•) andD n = n l ∂ ,l (•) are the tangential and normal surface differential operators, respectively, which emerges from the treatment of the variation of displacement gradient δu i,j on the boundary (see also [31,49]). Notably, in the above equation, as a direct result of the variational process, the generalized electrostatic stress σ ES ij appears as with σ M ij and τ M ijm,m the Maxwell stress and the electrostatic stress corresponding to the strain and polarization gradients, respectively, defined as Upon substituting Equations (11) and (8) into Equation (7) and utilizing the arbitrariness of δu i , δφ and δP i , we can obtain the governing equations and the boundary conditions with the remark that mechanical displacement and polarization do not exist in the surrounding air environment Dirichlet boundary conditions D n (u i ) = u i,l n l =ū i,l n l , on a u,n (15c) Neumann boundary conditions

Alternative Formulation
The above boundary value problem is formulated from the total energy density U whose strain ε ij , strain gradient ε ij,k , polarization P i and polarization gradient P i,j are independent variables. Alternatively, we can also choose the electric Gibbs free energy (which is identical to the electric enthalpy in isothermal system), whose independent variables are strain ε ij , strain gradient ε ij,k , electric field E i and electric field gradient E i,j [38,40,50]: Consequently, the constitutive relations are given as where σ ij , D i , τ ijk and V ij are the (local) stress, electric displacement, higher-order stress and higher-order electric displacement, respectively, while κ ij , b ijkl , c ijkl , e ijk and µ ijkl are material tensors [38]. These relations also imply that the variation of electric Gibbs free energy can be expressed as which then can be used to determine the equilibriums from variational principle where W is again the external work, but instead of having the contribution of the external electric field and higher-order electric field as in the previous formulation, higher-order surface chargesq are prescribed on the boundary, so that its variation δW reads By analogy with Equation (10), migratory variation of electric field and electric field gradient are obtained as Now, substituting Equation (22) into Equation (20) with the use of Equation (19) and employing Gauss divergence theorem, the physical variational the electric Gibbs free energy can be obtained where σ M ij and τ M ijm,m are the Maxwell stress and the higher-order electrostatic stress that are direct results of the variational process and expressed as Note that Equations (24a) and (13a) are identical. Finally, by substituting Equations (23) and (21) into Equation (20), we can obtain the governing equations and boundary condition as Dirichlet boundary conditions D n (u i ) = u i,l n l =ū i,l n l , on a u,n (26c) D n (φ) = φ ,l n l =φ ,l n l , on a φ,n (26d) Neumann boundary conditions In the above equations, although it is not trivial, we need to bear in mind that in vacuum displacement and polarization do not exist; hence, the quantity For the sake of clarity, we compare the two formulations in Table 1. Let us denote the strong form derived from the internal energy U and the electric Gibbs energy G in Sections 2.1 and 2.2 as P-formulation and E-formulation, respectively. As compared to the classical (local) theory, in both formulations, strain gradient ∇ results its conjugate variable, double stress τ ijk , to be cast in the balance of linear momentum as well as surface traction and surface higher-order traction. In addition to the classical Dirichlet boundary condition, the displacement gradient in the normal direction is also specified on the boundary. In the same manner, the electric field gradient dependence in electric Gibbs energy of E-formulation causes its conjugate variable, higher-order electric displacement Q ij , to appear in the modification of Gauss law and surface charge. The electric potential gradient in normal direction is also specified on the boundary. On the other hand, in P-formulation, when polarization gradient is considered, its conjugate variable V ij enters the intramolecular force balance equation and higher-order surface polarization, while the Gauss law and surface charge condition are unchanged as compared to classical theory. Moreover, as mentioned above, both formulations give identical Maxwell stress tensor σ M ij , whereas the higher-order Maxwell stress are given in their corresponding electric independent variables. It is worth noting that the electrostatic stress tensor σ ES ij in vacuum is reduced to which is the usual Maxwell stress tensor in electromagnetism. Furthermore, we remark that different formulations from the literature can be deduced from the generalized ones presented here. For instance, in the absence of electrostatic stress σ ES ij and the polarization gradient P i,j (in P-formulation) or the electric field gradient E i,j (in E-formulation), we can retrieve the strong form given by Sharma [30], Mao [39], Abdollahi [38], Deng [40] and Nguyen [43]. We summarize the two reduced formulations in Table 2. It should also be noted that, in P-formulation, when the external electric field is also omitted, the intramolecular force balance is reduced to E L i = φ ,i or the usual relation between electric field and electric potential E i = −φ ,i can be retrieved and the strong form of two formulations are identical. To this end, which formulation to be used for numerical computation is simply a matter of choice, as the weak form of both formulations are identical and can be written as for E-formulation. Note that Equation (29) will be solved with respect to displacement and electric potential fields. Hence, the electric polarization needs to be expressed in terms of the electric potential by utilizing the constitutive relation in Equation (2c), i.e., On the other hand, as the electric field E i = −φ ,i is an independent variable in E-formulation, it is straightforward to solve for the displacement and electric potential fields from Equation (30). Therefore, in the following sections, we employ different numerical methods to solve the weak form in Equation (30) of simplified E-formulation. When full formulations (Table 1) are chosen, the electric Gibbs free energy is still more favorable. In this case, the expression of the electric polarization in terms of electric potential is done in both domain and surface integrals. In other words, the application of the surface polarization in the P-formulation is equivalent to that of the potential normal derivative in E-formulation.
Displacement normal derivative on a u,n u i,l n l =ū i,l n l u i,l n l =ū i,l n l Potential normal derivative on a φ,n φ ,l n l =φ ,l n l Surface polarization on a P P i =P i

Neumann BCs
Surface traction on a σ Higher-order traction on a τ τ ijm n m n j =τ i τ ijm n m n j =τ i Higher-order surface charge on a Q Q ij n i n j =q Higher-order electric field on a V V ij n j =v i Table 2. Comparison of reduced strong form. [30,39] References [38,40,43] Balance Displacement normal derivative on a u,n u i,l n l =ū i,l n l u i,l n l =ū i,l n l

Neumann BCs
Surface traction on a σ Higher-order traction on a τ τ ijm n m n j =τ i τ ijm n m n j =τ i

Computational Modeling
The modeling of flexoelectricity is to solve the boundary value problem described in Table 2 or to minimize the weak forms (Equation (29) or Equation (30)). To do this numerically, the existence of strain gradient entails the chosen numerical methods must ensure at least C 0 -continuity of the second-order derivatives of the displacement field by means of higher-order (at least C 1 -continuity) shape functions. In the following, we summarize the numerical methods that have been employed in flexoelectricity, including meshless, mixed FE and IGA.

Meshfree Method
Some of the initial prominent works on computational evaluation of flexoelectric structures is by adopting a meshfree method using local maximum entropy (LME) approximants, as was done by Abdollahi et al. [38]. The meshfree formulation is validated by comparing the analytical and numerical energy conversion factor in a one-dimensional problem. Inspired by the experimental works of Ma and Cross, the LME based meshfree method is adopted to analyze a two-dimensional pyramid structure in [38] and a three-dimensional pyramid structure in [50]. The converse flexoelectricity that leads to strain due to indentation of atomic force microscopy tip onto a SrTiO 3 sample is experimentally studied and simulated using LME-based meshfree method in [51]. Besides the works on LME, Bo He et al. [52] presented analysis of 2D flexoelectric structures by element free Galerkin method which uses moving least square (MLS) approximants. Nevertheless, more works on flexoelectricity using meshfree methods are still required, as both MLS and LME shape functions lack Kronecker delta property, which leads to additional treatment on imposing Dirichlet boundary conditions. Moreover, in the multi-material problem that possesses material interfaces, some special techniques such as domain partitioning [53], Lagrange multiplier [54] and Jump function [55] are utilized to capture material discontinuities. On the other side, the higher-order Dirichlet boundary condition has not been mentioned in the meshfree approach in flexoelectricity.

Mixed Finite Element Method
The mixed finite element method offers the advantage of capturing C 1 continuity in the domain by utilizing C 0 finite elements. The first work on utilizing mixed finite element was by Mao et al. [39] for analyzing flexoelectric plate and crack in periodic flexoelectric structure, in which they proposed a mixed FE element denoted as I9-87. The formulation for flexoelectricity adopted in [39] is based on P-formulation and so the electrical degrees of freedom include polarization in addition to electric potential. The mechanical DOFs comprise displacements and relaxed displacement gradients and kinematic constraints are imposed by Lagrange multipliers. Mixed FE was also implemented for E-formulation by the authors of [40,56,57]. Nanthakumar et al. [56] adopted a nine-noded quadrilateral element with 54 DOFs to analyze and optimize two-dimensional flexoelectric structures. Deng et al. [40] presented two triangular elements, T37 and T45, and two quadrilateral elements, Q47 and Q59. A hexahedral element was presented by Deng et al. [57], who utilized the mixed FE to analyze a micro-pyramid structure.
The main advantage of mixed FE is its simplicity in imposing higher-order Dirichlet boundary conditions due to the displacement gradient degrees of freedom; hence, it can be directly implemented in FE commercial software. Moreover, due to its C 0 -approximation, mixed FEM can be used to alleviate the modeling of interface and composite flexoelectric structure. However, the method suffers one major drawback: expensive computational cost due to the many degrees of freedom per element, especially in 3D.

IGA Approach
Due to the flexibility of manipulating continuity of the basis functions, IGA has been used in solving flexoelectricity problems for topology optimization [41], soft dielectric material [42] and semiconductor [44]. In those works, NURBS basis functions are chosen to approximate both geometry and field variables (displacement and electric potential fields). By using knot insertion technique [58,59], the desired continuity order of NURBS basis functions can be obtained, as shown in Figure 1. Moreover, by controlling the continuity, IGA approach can be used to model interface [42], in which the geometrical approximation across the interface is C 0 -continuity while C 1 -continuity is required in the domains as usual. Besides, in topology optimization works [41], the voids are assumed to be solid material whose density is much smaller than the host material; however, C 1 continuity is imposed across the void-solid phase. A more general treatment for modeling discontinuity across the internal boundary can be adopted from the works in [60][61][62]. The approach has been applied for porous media, composites, multi-field and multi-material problems. In these works, the interface element has been applied to construct C 0 -continuous interpolation along both sides of the redefined interface. The element enables the possibility in enforcing the boundary condition via adding mechanical traction defined by a constitutive relation at the internal discontinuity into the equilibrium equation.
where [[]] indicates the jump. Such conditions of weakly-discontinuity displacement can be implemented with suitable enrichment functions as in XFEM [63,64]. Moreover, in strain gradient elasticity problems, additional continuity conditions associated with higher-order terms such as displacement normal derivative and double traction [65,66] [[u i,l n l ]] = 0 (32a) should be satisfied along the interface. Furthermore, the interfaces in microstructure problems are often considered to be imperfect (displacement and/or traction fields are not continuous along the interface) [67,68]. However, these imperfect interface models are limited to local elasticity and do not consider strain gradients. Although the interface element was applied for mechanical problems, the same concept can be inherited and applied for multi-field problems such as piezoelectric and flexoelectric composites.
It is worth discussing the imposition of the boundary conditions in IGA approach. While the Dirichlet boundary condition can be imposed by means of least-square method [69] due to the lack of Kronecker-delta property of NURBS basis functions, the displacement gradients conditions can be applied by imposing constraint on displacement degrees of freedom with Lagrange multipliers. This idea is often used in solving strain-gradient elasticity [70] but has not been discussed in flexoelectricity modeling works. Numerical implementation of the displacement gradients conditions is shown in the next section.
In the reduced formulations without consider electric field or polarization gradient, the electric potential gradient or surface polarization needs not to be imposed. However, for modeling flexoelectric transducers, it is necessary to model the electrode on the surface of the structure by imposing the equipotential condition. This can be done by utilizing the Lagrange multiplier, as we show in the next section. Note that the equipotential has not been used in previous works on flexoelectricity using IGA. It should be remarked that recently Codony et al. [46] proposed the immersed boundary method on hierarchical B-spline for flexoelectricity problems with higher efficiency in terms of complex domain shapes smf boundary conditions handling. Moreover, the non-local conditions on non-smooth boundary (which we have neglected) was also treated.
In Table 3, we present a comparison of different numerical methods on some special features in flexoelectricity modeling, including continuity handling, computational cost (via degrees of freedom per element), the ability to enforce higher-order Dirichlet boundary conditions and material discontinuity approximation. Each numerical method comes with pros and cons and should be considered based on the desired application. For instance, mixed FE is more suitable for a simplified 2D flexoelectric problems, whereas flexoelectric domains consisting of voids or inclusions can be modeled with meshless methods.

Implementations
In this section, we present the computation procedure of IGA and mixed FE approaches to solve problems governed by flexoelectricity, which can serve as a numerical tool for simulating the nano-energy harvesters, sensors or actuators. For the reasons mentioned above, we consider the simplified E-formulation in Table 2, which was also the choice of the authors of [38,40], hence reasonably suitable for comparison and validation.

Mixed FEM
The requirement of C 1 continuous basis functions due to the fourth-order PDEs of flexoelectricity excludes (or at least complicates) the use of classical Lagrange polynomials. too ensure C 1 continuity displacement gradients, ψ are introduced in the formulation and the kinematic constraint, ψ = ∇u, is imposed using Lagrange multipliers. The weak form of mechanical and electrostatic equilibrium is obtained by finding the u ∈ u =ū on a u , whereτ = ∂U ∂η andη is the third-order tensor related to relaxed displacement gradient ψ asη = 1 2 (ψ jk,i + ψ ik,j ). The constraint ψ = ∇u is imposed in a weighted residual manner by including δ Ω (ψ − ∇u) : λdΩ in Equation (33). A nine-noded isoparametric element was proposed by Nanthakumar et al. [56]. The degrees of freedoms are u 1 and u 2 at all nodes, relaxed displacement gradients ψ 11 , ψ 12 , ψ 21 and ψ 22 at the four corner nodes and electric potential φ at the four corner nodes. In addition to these DOFs, the element has four Lagrange multipliers λ 11 , λ 12 , λ 21 and λ 22 at the four corner nodes. The flexoelectric element is shown in Figure 2. The biquadratic shape functions, N q , are used to interpolate the displacement DOFs u 1 and u 2 . Bilinear shape functions, N l , are used to interpolate the relaxed displacement gradients, Lagrange multipliers and electric potential, φ. The finite element approximation is as shown below, Substituting Equations (35) and (36), into Equation (34) yields the following expression in terms of nodal DOFs: where The algebraic equations to be solved to obtain displacement and electric potential in a flexoelectric structure are, 

IGA Approach
Within the Galerkin method, both the trial and test functions are approximated by NURBS basis functions in which N u and N φ are the matrices containing the local NURBS basis functions of element e defined as where N i is the NURBS basis function at control point ith and n is the number of control points of an element. In the above equation,ũ,φ, δũ and δũ are the nodal control variables associated with displacement and electric potential of the trial and test functions, respectively, which can be written in matrix form asũ in which B u , B φ , H u and H φ contain the first-and second-order partial derivative of the NURBS basis functions, respectively, defined as Note that matrices H u and H φ contain second-order derivative with respect to the physical coordinate of the basis functions. The details are presented in Appendix A. Now, based on the discretized finite elements, the weak form in Equation (30) can be re-expressed as (e) δφ Tσ * da . (48) Finally, by substituting Equations (39) and (45) into Equation (48) and utilizing the arbitrariness of the test functions, a system of linear equations can be obtained as where the stiffness matrices and force vectors are computed as Note that the flexoelectric tensor µ for cubic symmetry crystal κ is given as The symmetry of the flexoelectric tensor can be found in [71,72]. For cubic crystal, the non-zero components of flexoelectric tensor are µ 1111 , µ 1221 and µ 1122 while other coefficients with cycling indices are equal (e.g., µ 1221 = µ 1331 ) and with odd number of indices are zeros (e.g., µ 1112 = 0) [73]. In matrix notation, the coefficients µ iiii , µ jiij and µ iijj are denoted as µ 11 , µ 12 and µ 44 , respectively [38].

Hollow Cylinder
In this first numerical example, we consider an infinitely long hollow cylinder subject to mechanical displacement and electric potential at the inner and outer surface, as depicted schematically in Figure 3a. Taking advantage of the axisymmetry, a quarter of the hollow cylinder is modeled as shown in Figure 3b, in which symmetry boundary conditions are enforced at x 1 = 0 and x 2 = 0. This problem is an oft-used benchmark in gradient elasticity [74,75], and has been extended to flexoelectricity [39,40]. In our work, we adapt both internal energy U from [39] and electric Gibbs energy G from [40] as follows and where λ and µ are the two Lamé parameters, l is the length scale, µ 12 and µ 44 are two non-zero flexoelectric coefficients, andf 1 andf 2 are two flexoelectric coupling coefficients that related to the flexoelectric coefficients by the susceptibility χ = κ/ε 0 − 1 such that µ 12 = χf 12 , µ 44 = χf 44 , α and κ are the reciprocal susceptibility and dielectric permittivity, also related to one another by the susceptibility. Furthermore, due to the geometry and boundary condition, the displacement as well as the electric potential depend only on the radial r of the hollow cylinder such as where l 2 0 = l 2 + µ 2

11
(λ+µ)κ , I 1 and K 1 are the modified Bessel functions of the first and second kind, respectively. The coefficients c 1 , c 2 , c 3 , c 4 , c 5 , c 6 are determined from the boundary conditions For numerical modeling, both formulations in Table 2 are solved with IGA, henceforth denoted as IGA(P) and IGA(E) for the formulation obtained from electric enthalpy and Gibbs energy, respectively. In addition, we also demonstrate results from second formulation with mixed FEM. For the IGA approach, we employ 25 × 25 second-order NURBS basis functions that have C 1 continuity to model a quarter of the hollow cylinder, as shown in Figure 3b, and impose additional symmetric boundary condition on the left and bottom edges. While the homogeneous and inhomogeneous Dirichlet boundary conditions are trivial to impose [69], the boundary conditions for the displacement gradients require constraining the degrees of freedom associated with the boundary control points and their neighbors. For instance, the symmetry boundary conditions on the bottom edge in Figure 3b read u 1,2 where Equation (56c) is automatically fulfilled as a result of Equation (56a). To impose Equation (56b), the nodal parameters (associated with the displacement in x 1 direction u x 1 ) of the two bottom layer of control points are identical, i.e., u bottom . This can be implemented by utilizing Lagrange multipliers. On the other hand, the displacement gradient boundary conditions can be imposed effortlessly for the mixed FEM with the nodal displacement gradient. Table 4. Material parameter for the problem depicted in Figure 3a.  and r 0 = 20µm, respectively. The inner surface is imposed by radial displacement u r i = 0.045µ m and grounded, i.e., φ r i = 0V, whereas the radial displacement u r o = 0.05µm and electric potential φ r o = 1V are applied on the outer surface. Symmetry boundary conditions are enforced, such that u 1 x 1 =0 = u 1,2 x 1 =0 = u 2,1 x 1 =0 = 0 and u 2 x 2 =0 = u 1,2 x 2 =0 = u 2,1 x 2 =0 = 0. Figure 4 presents numerical results of displacement and electric potential obtained from three approaches, showing good agreement between them. For a closer look, numerical results of radial displacement u R , electric potential φ, radial strain ε rr and circumferential strain ε θθ along the radial direction with angle θ = 45 o are presented and compared with the analytical solutions in Figure 5. Excellent agreement between the numerical and analytical results can be seen.

Cantilever Beam
In this example, we consider a 2D cantilever beam subjected to a point load on the right-end and having two different electrode configurations, which serve as a demonstration for a nanogenerator [38]. As depicted in Figure 6, the beam is constrained at the left-end, whereas two different electrode configuration are set up for two electrical boundary conditions. In Type-1 boundary condition, the electrode at the right-end of the beam is grounded, i.e., the electrical potential φ is prescribed as zero. In Type-2 condition, an electric potential difference is applied between the top and bottom electrode, e.g., the electric potential on the top is fixed to zero, whereas the bottom is induced to have an unknown equipotential V, which is solved by utilizing Lagrange multiplier. The material parameter in this study is single barium titanate (BaTiO 3 ) crystal adopted from [38] and presented in Table 5. The mechanical deformation can be transform into electrical energy via the electromechanical coupling effect in the nanogenerator. The performance of this energy conversion can be evaluated by the electromechanical coupling factor k eff given as the ratio between the electrostatic energy E elec and the elastic energy E elas Consequently, the normalized piezoelectric constant is defined as where k piezo is obtained as k eff without considering the flexoelectric effect. Physically, the normalized piezoelectric constant indicates the enhancement of flexoelectric effect on the energy conversion. and excellent agreement is shown, which illustrates the great enhancement of flexoelectricity with the decrease in beam thickness. In addition, this enhancement can also be seen in the full 2D model, which is also in good agreement with the work of Abdollahi et al. [38].  Table 5. Material parameters for the problem depicted in Figure 6. We further investigate the difference in the normalized effective piezoelectric constantē with two different electrical boundary conditions and present the numerical results in Figure 8. In both boundary condition setups, as the beam thickness decreases, the normalized effective piezoelectric constant increases, which is expected due to the size effect nature of flexoelectricity. Notably, Type-1 boundary condition seems to be more effective than Type-2. Furthermore, both IGA and mixed FEM approach give the same prediction on the performance of two boundary condition types. For a closer look, we present the electric potential distribution for two electrical boundary conditions computed from IGA and mixed FEM in Figure 9 with the height and length of the beam being h = 1.2468 × 10 −7 (m) and L = 10h, respectively. The results from the two methods are in agreement. For Type-1 boundary condition, the electric potential difference is highest at the left end of the beam, where strain gradient is largest. As for the Type-2 boundary condition, the induced electric potential at the top electrode is around −9 µV.  Additionally, we also carried out a three-dimensional simulation of a cantilever beam under Type-2 electrical boundary condition by using IGA. The 3D beam is assumed to have square cross-section, where the width is equal to the height h. The length of the beam as well as the applied force on the right surface are unchanged as in the two-dimensional case. Figure 10 shows the induced electric potential. A similar response has been observed as in 2D case, where largest electric potential is at the left end. Nevertheless, the electric potential at the top electrode has the value of 3.21 V.  We also study the converse flexoelectric effect from the electrical boundary condition. While the mechanical point load is omitted, a uniform electric field with magnitude |E| = 8 MV/m is applied by setting the electric potential to be V = −8h (MV) on the bottom electrode. In this case, the beam behaves as an actuator and induces bending curvature, which is shown in Figure 11b. The calculated curvature of the bending beam due to converse flexoelectric effect is in good agreement with both numerical results from [38] and the experimental ones from Bursian and OI [76]. Further investigation on the distribution of electric field in the thickness direction of the beam is shown in Figure 11a, which is also in agreement with the result from [38]. Non-uniform distribution of electric field can be observed, especially at the region near the top and bottom edge of the beam, which yields high electric field gradient and consequently induces deformation. Note that the considered flexoelectric beams do not account for the resistive load [10] or the rectifying circuit [77,78]. In [77], a piezoelectric structure is analyzed along with its circuit connections. In addition, in [78], comparisons between the charge type Hamiltonian and voltage type Hamiltonian are performed to identify the output power and voltage of the piezoelectric energy harvester. Development of a fully combined flexoelectric generator and interface circuit might be of significance in the future.

Truncated Pyramid
This section presents the numerical results of a two-dimensional truncated pyramid under compression with different type of boundary conditions, as adopted from [38]. Figure 12 depicts the schematic of a truncated pyramid whose top surface is grounded with zero electric potential, whereas the bottom surface is attached with an electrode that results in the equipotential condition on it. Two mechanical boundary conditions are considered, namely the rigid boundary condition where the bottom surface is constrained in x 2 direction and the top surface is subjected to uniform load F, and the flexible boundary condition where both the bottom and top surfaces are subjected to uniform compression load F. For numerical modeling, the problem is discretized by quadratic C 1 -continuity NURBS basis functions, as shown in Figure 12b as well as by mixed FEM as in Figure 12c. (a) Schematic of a truncated pyramid. Numerical results of electric potential φ and compression strain ε yy for the rigid boundary condition from IGA and Mixed FEM are presented in Figure 13 and for flexible boundary condition from IGA are shown in Figure 14. Good agreement can be observed between the two methods as well as with the reference results [38]. In rigid model, our IGA and mixed FEM approaches predict the equipotential of bottom electrode to be 2.7 and 2.5 mV, respectively, while the reference result is 2.6 mV [38]. Similarly, in flexible model, the bottom equipotential is 5.7 mV, agreeing with the value of 5.6 mV from [38].
Additionally, a three-dimensional version of the truncated pyramid is also studied. The geometric parameters are kept as in 2D case, so that the top and bottom surfaces now have dimension of a 1 × a 1 and a 2 × a 2 , respectively. The flexible boundary condition is considered, where the applied force F is identical from the 2D case. The numerical result from IGA is shown in Figure 15. As can be seen in Figure 15c, the electric potential in 3D model is very similar to that of 2D model. However, the induced potential on the top electrode is now 6.918 V, three orders of magnitude greater than the result predicted from 2D model.

Discussion
Since its discovery, flexoelectricity bears significance because of its potential to aid developing novel electromechanical coupling devices. To fully utilize this effect, both experimental and simulation approaches are necessary. In view of simulation, in this paper, we provide a generalized formulation of flexoelectricity as the extension from gradient elasticity. Remarkably, the boundary value problems that involve fourth-order PDE in flexoelectricity, necessitating C 1 continuity, are obtained from electric enthalpy as well as electric Gibbs free energy. To overcome the high-order continuity, we employ isogeometric analysis and mixed finite element, where detailed implementations are reported. We solved benchmark problems using these methods on 2D and 3D flexoelectric problems. The source codes are provided and can be used to study more complex problems.