Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Computational Modeling of Flexoelectricity—A Review

Computational Modeling of Flexoelectricity—A Review energies Review 1,2 3 3 Xiaoying Zhuang , Binh Huy Nguyen , Subbiah Srivilliputtur Nanthakumar , 3 4 4, Thai Quoc Tran , Naif Alajlan and Timon Rabczuk * Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City 758307, Vietnam; xiaoying.zhuang@tdtu.edu.vn Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City 758307, Vietnam Institute for Continuum Mechanics, Leibniz Universität Hannover, Appelstr. 11, 30167 Hannover, Germany; binh.nguyen@ikm.uni-hannover.de (B.H.N.); nanthakumar@ikm.uni-hannover.de (S.S.N.); tran@ikm.uni-hannover.de (T.Q.T.) Department of Computer Engineering, College of Computer and Information Sciences, King Saud University, Riyadh 11543, Saudi Arabia; timon.rabczuk@uni-weimar.de * Correspondence: timon.rabczuk@uni-weimar.de Received: 21 January 2020; Accepted: 26 February 2020; Published: 12 March 2020 Abstract: 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. Keywords: flexoelectricity; numerical methods; modeling 1. 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–5], actuators [6,7] and nano-energy harvesters [8–13]. It can also explain physical phenomena such as strain gradient-driven polarity control [14–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–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 Energies 2020, 13, 1326; doi:10.3390/en13061326 www.mdpi.com/journal/energies Energies 2020, 13, 1326 2 of 29 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–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–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 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–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 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. 2. Ground Work Formulation 2.1. Physical Variational Principle Hu and Shen [31] extended the variational principle from Toupin [48] by the physical variational principle proposed by Kuang [32–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 e , strain gradient e , polarization P ij ij,k i and polarization gradient P , i,j 1 1 1 1 U e , e , P , P = a P P + b P P + c # # + d P # + f # P + h P # + g # # , (1) ij ij,k i i,j ij i j ijkl i,j k,l ijkl ij kl ijk i jk ijkl ij k,l ijkl i jk,l ijklmn ij,k lm,n 2 2 2 2 Energies 2020, 13, 1326 3 of 29 from which the constitutive relations are obtained as ¶U s = = c # + d P + f P (2a) ij ijkl kl ijk k ijkl k,l ¶# ij ¶U t = = h P + g # (2b) ijk ijkl l ijklmn lm,n ¶# ij,k ¶U E = = a P + d # + h # (2c) i ij j ijk jk ijkl jk,l ¶P ¶U V = = b P + f # , (2d) ij ijkl k,l ijkl kl ¶P i,j where s , E , t and V are the (local) stress, electric field, higher-order stress and higher-order ij i ijk ij electric field, respectively, while a , b , c , d , f , and h are the material tensors [32]. In other ij ijkl ijkl ijk ijkl ijkl words, the variation of internal energy can be written as dU = s d# + t d# + E dP + V dP . (3) ij ij ijk ij,k i i ij i,j The energy density is the sum of the internal energy and the internal of the Maxwell self-field [24] W = U + # f f , (4) 0 ,i ,i in which the electric Maxwell self-field can be defined as the gradient of the electric potential f MS E = f . (5) ,i Next, we can define the electric enthalpy as H = W E D = U # f f + f P (6) i i 0 ,i ,i ,i i We consider a dielectric solid occupying a volume v bounded by a surface a that separates from the 0  0 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 d H dv = dW, (7) ext f where W is the external work done by body force f , external electric field E and free charge r on the solid body and traction t , double-traction t , surface charge s and higher-order electric field v ¯ on i i i the surface boundary. Hence, its variation dW is defined as R  R R R R R ext n  f dW = f du + E dP dv + t du da + t ¯ D (du ) da + v ¯ dP da s df da r df dv , (8) s t V f i i i i i i i i i v i a a a a v in which the operator D is explained in the following section. Within physical variational principle approach, variation of the electric enthalpy can be further expressed as [32–35]. Z Z Z d H dv = dH dv + H du dv , (9) k,k v v v e 1 where H = E P + V P is the energy deducted from the total energy minus the pure deformation i i ij i,j energy. Moreover, the migratory variation is the essential feature of physical variational principle, Energies 2020, 13, 1326 4 of 29 in the sense that the virtual displacement not only produces the strain variation but also causes electric potential and polarization variations, such as [31] df = d f + d f = d f + f du , (10a) f u f ,i i df = d f + d f = d f + f du , (10b) ,i f ,i u ,i f ,i ,ij j dP = d P + d P = d P + P du , (10c) i P i i P i i,j j dP = d P + d P = d P + P du = d P + P du , (10d) i,j P i,j u i,j P i,j i,jk k P i,j i,kj k where d () and d () are the variations produced by the (local) variation of the electric potential and polarization, respectively, whereas d  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] Z Z h i ES d H dv = s + t s du dv ij,j ijm,mj i ij,j v v nh i h io ES t t ˜ ˜ + s t + s n + D (n )n n t D n t du da m m ij ijm,m j l j ijm ijm i ij l j + t n n D (du ) da ijm m j i Z Z ext + E V + f + E dP dv + V n dP da i ij,j ,i i ij j i v a Z Z (# f + P ) df dv + (# f + P )n df da , (11) 0 ,i i 0 ,i i i ,i v a t n ˜ ˜ where n is the normal vector. D () = (d n n ) ¶ () and D = n ¶ () are the tangential and i kl k l ,l l ,l normal surface differential operators, respectively, which emerges from the treatment of the variation of displacement gradient du on the boundary (see also [31,49]). Notably, in the above equation, as a i,j ES direct result of the variational process, the generalized electrostatic stress s appears as ij 1 1 ES s = V P + (E P + V P )d D f + # f f + f P d kj k,i k k kl k,l ij j ,i ,k ,k ,k k ij ij 2 2 M M = s + t , (12) ij ijm,m M M with s and t the Maxwell stress and the electrostatic stress corresponding to the strain and ij ijm,m polarization gradients, respectively, defined as 1 1 s = E P d D f + f + # f f + f P d (13a) k k ij j ,j ,i ,k ,k ,k k ij ij 2 2 t = V P + V P . (13b) kj k,j kl k,l ijm,m Upon substituting Equations (11) and (8) into Equation (7) and utilizing the arbitrariness of du , df and dP , 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 Energies 2020, 13, 1326 5 of 29 ES s t + s + f = 0, in v (14a) ij ijm,m i ij ,j ES s = 0, in v’ (14b) ij,j L ext V + E f + E = 0, in v (14c) ij,j ,i i i # f + P = r , in v (14d) 0 ,ii i,i f = 0, in v’ (14e) ,ii Dirichlet boundary conditions u = u ¯ , on a (15a) i i f = f, on a (15b) n u,n D (u ) = u n = u ¯ n , on a (15c) i i,l l i,l l P = P , on a (15d) i i Neumann boundary conditions ES t t s ˜ ˜ s t + [[s ]] n + D (n )n n t D n t = t , on a (16a) m m ij ijm,m j l j ijm ijm i ij l j t n n = t , on a (16b) ijm j i n (# [[f ]] + P ) = s , on a (16c) i 0 ,i i V n = v ¯ , on a (16d) ij j i 2.2. Alternative Formulation The above boundary value problem is formulated from the total energy density U whose strain # , strain gradient # , polarization P and polarization gradient P are independent variables. ij ij,k i i,j 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 # , strain gradient # , electric ij ij,k field E and electric field gradient E [38,40,50]: i i,j 1 1 1 1 G # , # , E , E = k E E b E E + c # # e E # m (E # # E ) (17) ij ij,k i i,j ij i j ijkl i,j k,l ijkl ij kl ijk i jk ijkl k ij,l ij k,l 2 2 2 2 Consequently, the constitutive relations are given as ¶G 1 s = = c # + m E (18a) ij ijkl kl ijkl k,l ¶# 2 ij ¶G 1 t = = m E (18b) ijk ijkl l ¶# 2 ij,k ¶G 1 D = = k E + m # (18c) i ij j ijkl k,l ¶E 2 ¶G 1 Q = = b E m # , (18d) ij ijkl k,l ijkl kl ¶E 2 i,j where s , D , t and V are the (local) stress, electric displacement, higher-order stress and ij i ijk ij higher-order electric displacement, respectively, while k , b , c , e and m are material ij ijkl ijkl ijk ijkl tensors [38]. These relations also imply that the variation of electric Gibbs free energy can be expressed as dG = s d# + t d# D dE + Q dE , (19) ij ij ijk ij,k i i ij i,j Energies 2020, 13, 1326 6 of 29 which then can be used to determine the equilibriums from variational principle d G dv = dW, (20) 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 charges q ¯ are prescribed on the boundary, so that its variation dW reads R  R R R R R ext n n  f ˜ ˜ dW = f du + E dP dv + t du da + t ¯ D (du ) da + q ¯D (df) da s df da r df dv . (21) s t Q f i i i i i i i v i a a a a v By analogy with Equation (10), migratory variation of electric field and electric field gradient are obtained as dE = d E + d E = d f + E du = d f + E du , (22a) i f i u i f ,i i,p p f ,i p,i p dE = d E + d E = d E + E du = d E + E du (22b) f u f p f p i,j i,j i,j i,j i,jp i,j i,pj 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 Z Z M M d G dv = s + t s t du dv ij,j i ijk,kj ij,j ijm,mj v v M M t t ˜ ˜ + s n t n + s n + t n + D (n )n n t D n t du da m m ij j ijk,k j j j l j ijm ijm i ij ijm,m l j Z Z + t n n D (du ) da + D + D df dv ijk k j i i,i ij,ji a v Z Z t t n ˜ ˜ ˜ + D n D n + D (n )n n Q D (Q n ) df da + Q n n D (df) da , (23) i i ij,j i l i j ij i ij j ij i j a a M M where s and t are the Maxwell stress and the higher-order electrostatic stress that are direct ij ijm,m results of the variational process and expressed as s = E D E D d (24a) i j k k ij ij t = E Q E Q d . (24b) ij ijm,m i,k kj k,l kl 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 ES s t + s + f = 0, in v (25a) ij ijm,m i ij ,j ES s = 0, in v’ (25b) ij,j D Q = r , in v (25c) i ij,j ,i D Q = f = 0, in v’ (25d) i ij,j ,ii ,i Dirichlet boundary conditions u = u ¯ , on a (26a) i i f = f, on a (26b) n u,n D (u ) = u n = u ¯ n , on a (26c) i,l l i,l l n f,n ˜ ¯ D (f) = f n = f n , on a (26d) ,l l ,l l Energies 2020, 13, 1326 7 of 29 Neumann boundary conditions ES t t s ˜ ˜ s t + [[s ]] n + D (n )n n t D n t = t , on a (27a) m m ij ijm,m j l j ijm ijm i ij l j t n n = t ¯ , on a (27b) ijk k j i t t  s ˜ ˜ [[D Q ]]n + D (n )n n Q D (Q n ) = s , on a (27c) i ij,j i l i j ij ij j l i Q n n = q ¯, on a (27d) ij i j 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 D Q is degenerated to f in i ij,j ,i vacuum, which also causes the jump [[D Q ]] on the boundary a . i ij,j 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 re results its conjugate variable, double stress t , to be cast in the ijk 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 , to appear ij 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 enters the intramolecular force balance equation and higher-order ij 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 s , whereas the higher-order Maxwell stress are given in their corresponding electric independent ij ES variables. It is worth noting that the electrostatic stress tensor s in vacuum is reduced to ij ES M s = s = # E E E E d , (28) 0 i j k k ij ij ij which is the usual Maxwell stress tensor in electromagnetism. Furthermore, we remark that different formulations from the literature can be deduced from ES the generalized ones presented here. For instance, in the absence of electrostatic stress s and the ij polarization gradient P (in P-formulation) or the electric field gradient E (in E-formulation), we can i,j i,j 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 = f ,i or the usual relation between electric field and electric potential E = f can be retrieved and the i ,i 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 R R R R R s du + t du + (# E + P )df dv = f du dv + t du da r df dv  s df da (29) 0 s s ij i,j ijm i,jm i i ,i i i i i v v a v a for P-formulation and Z Z Z Z Z s du + t du + D df dv = f du dv + t du da r df dv s df da (30) ij i,j ijm i,jm i ,i i i i i s s v v a v a 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 1 1 1 by utilizing the constitutive relation in Equation (2c), i.e., P = a f a d # a h # . i ,j jkl kl jklm kl,m ij ij ij Energies 2020, 13, 1326 8 of 29 On the other hand, as the electric field E = f is an independent variable in E-formulation, it i ,i 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. t t t ˜ ˜ ˜ Table 1. Comparison of strong form. L = D (n )n n t D n t , M = D (n )n n Q i l m j ijm m ijm l i j ij l j l D (Q n ). ij j U(e, P,re,rP) G(e, E,re,rE) ES ES Momentum in v s t + s + f = 0 s t + s + f = 0 ij ijm,m i ij i ijk,k ij ij ,j ,j 0 ES ES Maxwell in v s = 0 s = 0 ij,j ij,j Balance f f Gauss law in v # f + P = r D Q = r 0 ,ii i,i i ij,j ,i L ext Intramolecular force [24] V + E f + E = 0 ij,j ,i i i Gauss law in v f = 0 f = 0 ,ii ,ii Displacement on a u = u ¯ u = u ¯ i i i i ¯ ¯ Electric potential on a f = f f = f u,n Dirichlet BCs Displacement normal derivative on a u n = u ¯ n u n = u ¯ n i,l l i,l l i,l l i,l l f,n Potential normal derivative on a f n = f n ,l l ,l l Surface polarization on a P = P i i ES ES Surface traction on a s t + [[s ]] n + L = t s t + [[s ]] n + L = t ij ijm,m j i i ij ijk,k j i i ij ij Surface charge on a n # [[f ]] + P = s [[D D ]]n + M = s i 0 ,i i i ij,j i Neumann BCs Higher-order traction on a t n n = t ¯ t n n = t ¯ ijm m j i ijm m j i Higher-order surface charge on a Q n n = q ¯ ij i j Higher-order electric field on a V n = v ij j i t t ˜ ˜ Table 2. Comparison of reduced strong form. L = D (n )n n t D n t , M = m m i l j ijm ijm l j t t ˜ ˜ D (n )n n Q D (Q n ). l i j ij ij j l i U(e, P,re) G(e, E,re) References [30,39] References [38,40,43] ES ES Momentum in v s t + s + f = 0 s t + s + f = 0 ij ijm,m i ij ijk,k i ij ij ,j ,j f f Gauss law in v # f + P = r D Q = r 0 ,ii i,i i ij,j Balance ,i L ext Intramolecular force [24] V + E f + E = 0 ij,j ,i i i Gauss law in v f = 0 f = 0 ,ii ,ii Displacement on a u = u ¯ u = u ¯ i i i i ¯ ¯ Dirichlet BCs Electric potential on a f = f f = f u,n Displacement normal derivative on a u n = u ¯ n u n = u ¯ n i,l l i,l l i,l l i,l l s ES ES ¯ ¯ Surface traction on a s t + [[s ]] n + L = t s t + [[s ]] n + L = t ij ijm,m j i i ij j i i ijk,k ij ij Neumann BCs s Surface charge on a n # [[f ]] + P = s [[D D ]]n + M = s i 0 ,i i i ij,j i Higher-order traction on a t n n = t ¯ t n n = t ¯ m m ijm j i ijm j i 3. 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 -continuity of the second-order derivatives of the displacement field by means of higher-order (at least C -continuity) shape functions. In the following, we summarize the numerical methods that have been employed in flexoelectricity, including meshless, mixed FE and IGA. Energies 2020, 13, 1326 9 of 29 3.1. 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 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. 3.2. Mixed Finite Element Method The mixed finite element method offers the advantage of capturing C continuity in the domain by utilizing C 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 -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. 3.3. 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], 0 1 in which the geometrical approximation across the interface is C -continuity while C -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 continuity is imposed across the void–solid phase. A more general treatment for modeling discontinuity across the Energies 2020, 13, 1326 10 of 29 internal boundary can be adopted from the works in [60–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 -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. 1 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ξ ξ (a) (b) Figure 1. (a) C continuity quadratic NURBS basis functions from knot vector X = f0, 0, 0, 1/3, 1/3, 2/3, 2/3, 1, 1, 1g; and (b) C continuity quadratic NURBS basis functions from knot vector X = 0, 0, 0, 1/3, 2/3, 1, 1, 1 f g Perfect interfaces result in continuous displacement and traction fields along the interface, i.e., [[u ]] = 0 (31a) [[t ]] = 0, (31b) 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 n ]] = 0 (32a) i,l l [[t n n ]] = 0, (32b) ijm m j 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 R(ξ) R(ξ) Energies 2020, 13, 1326 11 of 29 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. 3.4. 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 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 continuity displacement gradients, y are introduced in the formulation and the kinematic constraint, y = ru, is imposed using Lagrange multipliers. The weak form of mechanical and electrostatic equilibrium u 1 u,n 1 is obtained by finding the u 2 u = u ¯ on a , u 2 H (v) , y 2 y = y on a , y 2 H (v) , f 2 f 1 2 u 1 f = f on a , f 2 H (v) , l 2 L (v), and similarly ffiu 2 ffiu = 0 on a , ffiu 2 H (v) , dy 2 u,n 1 f 1 2 dy = 0 on a , dy 2 H (v) , df 2 df = 0 on a , df 2 H (v) and dl 2 L (v), such that, Z Z dP = (s : #(du) + t ˆ .h ˆ(dy) D E(df)) dW + l : (dyrdu)dW W W Z Z Z (33) + (yru) : dldW du t dG du b dW = 0 W G W Z Z Z #(du) : C : #(u) dW #(du) : e E(f) dW h ˆ(dy) . h E(f) dW W W W Z Z Z . . . . . . + h ˆ(dy) . g . h ˆ(y)dW E(df) e : #(u) dW E(df) h . h ˆ(y) dW W W W (34) Z Z Z Z E(df) k E(f) dW + dy : l dW rdu : l dW dl : ru dW W W W W Z Z Z + dl : y dW = du t dG + du b dW W G W ¶U where t ˆ = and h ˆ is the third-order tensor related to relaxed displacement gradient y as h ˆ = ¶h ˆ (y + y ). The constraint y = ru is imposed in a weighted residual manner by including jk,i ik,j d (yru) : ldW in Equation (33). W Energies 2020, 13, 1326 12 of 29 Table 3. Comparison of different numerical methods. (a) and (b) indicate quadratic elements in 2D and 3D, respectively. Meshfree IGA Mixed FEM Continuity Depends on kernel functions Knot manipulation Nodal displacement gradient D.o.f Depends on domain of influence 27 (a) or 108 (b) 54 (Q54, [56]); 37 (T37), 47 (Q47), 45 (T45), 59 (Q59) [40]; 233 (3D) [57] Higher-order Dirichlet BCs - Lagrange-multiplier Direct imposition C Interface modeling Lagrange-multiplier, Enriched functions Multi-patch - Energies 2020, 13, 1326 13 of 29 A nine-noded isoparametric element was proposed by Nanthakumar et al. [56]. The degrees of freedoms are u and u at all nodes, relaxed displacement gradients y , y , y and y at the four 1 2 11 12 21 22 corner nodes and electric potential f at the four corner nodes. In addition to these DOFs, the element has four Lagrange multipliers l , l , l and l at the four corner nodes. The flexoelectric element 11 12 21 22 is shown in Figure 2. The biquadratic shape functions, N , are used to interpolate the displacement DOFs u and u . Bilinear shape functions, N , are used to interpolate the relaxed displacement 1 2 gradients, Lagrange multipliers and electric potential, f. The finite element approximation is as shown below, h h l u (x) = N u ; y (x) = N y ; i i å å i i=1 i=1,3,5,7 (35) h l h l f (x) = N f ; l (x) = N l å i å i i i i=1,3,5,7 i=1,3,5,7 h h l du (x) = N du ; dy (x) = N dy ; å i å i i i i=1 i=1,3,5,7 (36) h l h l df (x) = N df ; dl (x) = N dl å i å i i i i=1,3,5,7 i=1,3,5,7 7 6 8 4 1 2 Figure 2. A nine-noded flexoelectric element. DOFs at circles are u and u ; at squares are y , y , y , 1 2 11 12 21 y and f; and at triangles are l , l , l , and l . 22 11 12 21 22 Substituting Equations (35) and (36), into Equation (34) yields the following expression in terms of nodal DOFs: 0 1 Z Z Z T T T T T T T T @ A du B CB dW u + du B e B dW f + dy H h B dW f u f f u u W W 0 1 Z Z Z T T T T T T @ A + dy H g HdW y + df B eB dW u + df B h H dW y f f W W (37) Z Z Z T T T T l T l l df B kB dW f du B N dW l + dy N N dW l f yu W W W Z Z Z Z T T T l T l l T q T q dl N B dW u + dl N N dW y = du N t dG + du N b dW yu W W G W N Energies 2020, 13, 1326 14 of 29 2 3 ¶N 0 0 0 ¶x 6 7 ¶N 6 7 0 0 0 2 3 6 ¶x 7 ¶N " # l l 0 6 7 1 ¶N 1 ¶N ¶x ¶N 6 0 0 7 6 ¶N 7 2 ¶x 2 ¶x ¶x 6 7 where B = , B = , H = , B = u 4 5 f l yu ¶y ¶N 6 7 ¶N 0 0 0 q q 6 ¶y 7 ¶N ¶N ¶y 6 7 ¶y ¶x ¶N 6 7 0 0 0 ¶y 4 5 l l 1 ¶N 1 ¶N 0 0 2 ¶y 2 ¶y 2 3 ¶N ¶x 6 l 7 ¶N 6 0 7 ¶x 6 7 ¶N 6 7 4 5 ¶y ¶N ¶y The algebraic equations to be solved to obtain displacement and electric potential in a flexoelectric structure are, 2 32 3 2 3 K 0 K K u F uu uf u ul 6 76 7 6 7 0 K K K y 0 6 76 7 6 7 yy yf yl 6 76 7 = 6 7 (38) 4 54 5 4 5 K K K 0 f 0 fu fy ff K K 0 0 l 0 lu ly K = B CB dW ; uu u T T T K = B e B dW = K ; uf f u fu T T K = B N dW = K ul yu l lu T T K = H hB dW = K yf f fy K = B gB dW yy y l T K = N N dW = K yl l ly K = B kB dW ; ff f Z Z qT q F = N tdG + N b dW G W 3.5. IGA Approach Within the Galerkin method, both the trial and test functions are approximated by NURBS basis functions (e) (e) u = N u ˜, du = N du ˜, (39a) u u (e) (e) ˜ ˜ f = N f, df = N df, (39b) f f in which N and N are the matrices containing the local NURBS basis functions of element e defined as u f 2 3 N 0 0 i h i 6 7 1 2 n N = 0 N 0 , N = N N . . . N , (40a) 4 5 i u u u u u 0 0 N h i N = N N . . . N , (40b) 1 2 n Energies 2020, 13, 1326 15 of 29 where N is the NURBS basis function at control point ith and n is the number of control points of an element. In the above equation, u ˜ , f, du ˜ and du ˜ 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 h i 1 1 1 n n n u ˜ = u ˜ u ˜ u ˜ . . . u ˜ u ˜ u ˜ , (41) x y z x y z h i 1 1 1 n n n du ˜ = du ˜ du ˜ du ˜ . . . du ˜ du ˜ du ˜ , (42) x y z x y z h i 1 2 n ˜ ˜ ˜ ˜ f = f f . . . f , (43) h i 1 2 n ˜ ˜ ˜ ˜ df = df df . . . df (44) T T The approximation of strain field # = [# # # 2# 2# 2# ] , electric field E = [E E E ] , strain 11 22 33 23 13 12 1 2 3 h i h i T T ¶# ¶# ¶# ¶E ¶E ¶E gradient r# = and electric field gradient rE = can be obtained by using ¶x ¶y ¶z ¶x ¶y ¶z differential operators B ,B ,H and H on Equation (39) such that u f u f (e) (e) # = B u ˜, d# = B du ˜, (45a) u u (e) (e) ˜ ˜ E = B f, dE = B df, (45b) f f (e) (e) r# = H u ˜, dr# = H du ˜, (45c) u u (e) (e) ˜ ˜ rE = H f, drE = H df, (45d) f f in which B , B , H and H contain the first- and second-order partial derivative of the NURBS basis u f u f functions, respectively, defined as 2 3 ¶N 0 0 ¶x 6 7 ¶N i 2 3 6 0 0 7 1 2 n ¶B ¶B ¶B ¶y u u u 6 7 . . . ¶N 6 7 h i ¶x ¶x ¶x 6 7 0 0 1 2 n 6 7 i ¶B ¶B ¶B ¶z 1 2 n 6 u u u 7 B = , B = , H = (46) 6 7 B B . . . B . . . u u u ¶N ¶N u u u i i 4 5 ¶y ¶y ¶y 6 0 7 ¶z ¶y 1 2 n 6 7 ¶B ¶B ¶B u u u 6 ¶N ¶N 7 . . . i i ¶z ¶z ¶z 4 5 ¶z ¶x ¶N ¶N i i ¶y ¶x 2 3 1 2 n 2 3 ¶B ¶B ¶B f f f ¶N . . . ¶x ¶x ¶x ¶x h i 6 7 1 2 n 6 ¶N 7 6 ¶B ¶B ¶B 7 i 1 2 n i f f f B = , B = B B . . . B , H = (47) 4 5 6 7 f f . . . f f f f ¶y ¶y ¶y ¶y 4 5 ¶N 1 2 n ¶B ¶B ¶B f f f ¶z . . . ¶z ¶z ¶z Note that matrices H and H contain second-order derivative with respect to the physical coordinate u f 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 h i R  R R R R n T T e T T T T T f ¯ ˜ ˜ de s + dre t dE D dv = du f dv + du ˜ t da + f r dv df s ¯ da . (48) (e) s,(e) f,(e) e=1 v v a v a 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 " #" # " # K K u ˜ uu uf = , (49) K K f fu ff f Energies 2020, 13, 1326 16 of 29 where the stiffness matrices and force vectors are computed as n Z K = B cB dv , (50a) uu u (e) e=1 n Z T T K = B eB H mB dv , (50b) uf f f u u (e) e=1 T T T T K = B e B B m H dv , (50c) fu u u f f (e) e=1 n Z K = B kB dv , (50d) ff f (e) e=1 n Z Z T T f = N f dv + N t da , (50e) u u u s,(e) v a e=1 Z Z T f T f = N r dv + N t da . (50f) u f f,(e) v a e=1 Note that the flexoelectric tensor m for cubic symmetry crystal k is given as 2 3 m m m 0 0 0 0 0 0 0 0 m 0 0 0 0 m 0 11 12 12 44 44 6 7 m = (51) 4 0 0 0 0 0 m m m m 0 0 0 0 0 0 m 0 05 44 12 11 12 44 0 0 0 0 m 0 0 0 0 m 0 0 m m m 0 0 0 44 44 12 12 11 The symmetry of the flexoelectric tensor can be found in [71,72]. For cubic crystal, the non-zero components of flexoelectric tensor are m , m and m while other coefficients with cycling 1111 1221 1122 indices are equal (e.g., m = m ) and with odd number of indices are zeros (e.g., m = 0) [73]. 1221 1331 1112 In matrix notation, the coefficients m , m and m are denoted as m , m and m , respectively [38]. iiii jiij iijj 11 12 44 4. Numerical Examples 4.1. 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 = 0 and x = 0. 1 2 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 1 1 1 ˆ ˆ U = l# # + m# # + l l# # + 2m# # + f # P + 2 f # P + aP P (52) jj kk jk jk jj,i kk,i jk,i jk,i 12 jj,i i 44 ij,i j i i 2 2 2 and 1 1 1 G = l# # + m# # + l l# # + 2m# # m # E 2m # E kE E , (53) jj kk jk jk jj,i kk,i jk,i jk,i 12 jj,i i 44 ij,i j i i 2 2 2 where l and m are the two Lamé parameters, l is the length scale, m and m are two non-zero 12 44 ˆ ˆ flexoelectric coefficients, and f and f are two flexoelectric coupling coefficients that related to the 1 2 ˆ ˆ flexoelectric coefficients by the susceptibility c = k/# 1 such that m = c f , m = c f , a and k are 0 12 12 44 44 the reciprocal susceptibility and dielectric permittivity, also related to one another by the susceptibility. These material parameters can be found in Table 4. Furthermore, due to the geometry and boundary Energies 2020, 13, 1326 17 of 29 condition, the displacement as well as the electric potential depend only on the radial r of the hollow cylinder such as 1 r r u (r) = c r + c + c I + c K , (54a) r 1 2 3 1 4 1 r l l 0 0 m du u (r) r r f(r) = c ln r + c + + , (54b) ( ) 5 6 k dr r 2 2 where l = l + , I and K are the modified Bessel functions of the first and second kind, 1 1 (l+m)k respectively. The coefficients c , c , c , c , c , c are determined from the boundary conditions 1 2 3 4 5 6 u = u , u = u , (55a) r r r r r=r r=r i o f = f , f = f , (55b) r r r r r=r r=r i o t = t = 0 (55c) rrr rrr r=r r=r i o 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 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 = 0, (56a) x =0 u = 0, (56b) 1,2 x =0 u = 0, (56c) 2,1 x =0 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 direction u ) of the two bottom layer adjacent-bottom bottom of control points are identical, i.e., u = u . This can be implemented by utilizing x x 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. Y n l m m k 12 44 9 6 6 6 9 139 10 Pa 0.3 2 10 m 1 10 C/m 1 10 C/m 1 10 F/m Energies 2020, 13, 1326 18 of 29 −5 x 10 Element edge Control Points 1.8 1.6 1.4 1.2 0.8 0.6 0.4 0.2 0 0.5 1 1.5 2 1 −5 x 10 (b) IGA discretization (a) Schematic of a hollow cylinder. Figure 3. Schematic of an infinitely long hollow cylinder, where the inner and outer radii are r = 10 m and r = 20 m, respectively. The inner surface is imposed by radial displacement u = 0.045  m 0 r and grounded, i.e., f = 0 V, whereas the radial displacement u = 0.05 m and electric potential r r i o f = 1 V are applied on the outer surface. Symmetry boundary conditions are enforced, such that u = u = u = 0 and u = u = u = 0. 1 1,2 2,1 1,2 2,1 x =0 x =0 x =0 x =0 x =0 x =0 1 1 1 2 2 2 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 , electric potential f, radial strain # and circumferential strain # along the radial rr R qq direction with angle q = 45 are presented and compared with the analytical solutions in Figure 5. Excellent agreement between the numerical and analytical results can be seen. 4.2. 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 f 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 ) 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 given as the ratio between the electrostatic energy E and eff elec the elastic energy E elas E k E dv elec 2 v k = = R , (57) eff e C e dv elas 2 v Consequently, the normalized piezoelectric constant is defined as eff e ¯ = , (58) piezo where k is obtained as k without considering the flexoelectric effect. Physically, the normalized piezo eff piezoelectric constant indicates the enhancement of flexoelectric effect on the energy conversion. 2 Energies 2020, 13, 1326 19 of 29 Figure 7 demonstrates the normalized piezoelectric constant e ¯ as the normalized thickness h = e h/m varies under Type-1 configuration. The numerical results in the simplified case, where only 31 12 k and m are non-zeros to mimic the 1D problem, are compared to the analytical solution where 33 12 e +12(m /h) k = and excellent agreement is shown, which illustrates the great enhancement of eff k Y 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]. u r Phi 5.000e-08 1.000e+00 0.66667 4.8333e-8 0.33333 4.6667e-8 4.500e-08 0.000e+00 (a) IGA(E): Distribution of u (b) IGA(E): Distribution of f -8 Radial displacement x10 [m] Electric potential [V] 4.4674 4.5739 4.6804 4.7869 4.8935 - 0.2971 - 0.0375 0.2220 0.4815 0.7410 1.000 5.000 (d) IGA(P): Distribution of f (c) IGA(P): Distribution of u 0.8 0.6 0.4 0.2 46 0 -0.2 (e) Mixed FEM: Distribution of u (f) Mixed FEM: Distribution of f Figure 4. Distribution of radial displacement and electric potential. Energies 2020, 13, 1326 20 of 29 −8 x 10 5.1 Analytical Analytical IGA−E formulation 0.8 IGA−E formulation IGA− P formulation IGA− P formulation Mixed FEM 4.9 Mixed FEM 0.6 4.8 0.4 4.7 0.2 4.6 4.5 −0.2 4.4 −0.4 1 1.2 1.4 1.6 1.8 2 2.2 1 1.2 1.4 1.6 1.8 2 2.2 r −5 r −5 x 10 x 10 −4 −3 x 10 x 10 Analytical IGA−E formulation 4.5 IGA− P formulation Meshfree 3.5 2 Analytical IGA−E formulation IGA− P formulation Meshfree 2.5 −2 −4 2 1 1.2 1.4 1.6 1.8 2 2.2 1 1.2 1.4 1.6 1.8 2 2.2 r r −5 −5 x 10 x 10 Figure 5. Comparison of the exact and numerical results for radial displacement u , electric potential f, radial strain # and circumferential strain # as functions of radius r. rr ff Element edge Control Points -6 0 0.5 1 1.5 -5 (b) x V (a) Figure 6. (a) Schematic of a hollow cylinder. Two different types of electrical boundary condition are imposed. The upper figure describes Type-1 boundary where zero electric potential is imposed on the right edge. Type-2 boundary is demonstrated in the lower figure where zero electric potential is applied on the top edge and an unknown equipotential (due to the electrode) on the bottom edge. (b) IGA discretization (top) and mixed FEM discretization (bottom). Table 5. Material parameters for the problem depicted in Figure 6. Y n e m k k F 31 12 11 33 9 2 6 9 9 6 100 10 Pa 0.37 -4.4 C/m 1 10 C/m 1110 F/m 12.4810 F/m 10010 N ǫ u rr r φφ φ Energies 2020, 13, 1326 21 of 29 Analytical 1D (Flexo+Piezo) IGA−1D (Flexo+Piezo) Analytical 1D (Flexo) IGA−1D (Flexo) IGA− 2D (Flexo) Mixed FEM (Flexo+Piezo) Mixed FEM (Flexo) 0 1 2 3 4 5 6 7 8 Normalized thickness Figure 7. Effective normalized piezoelectric constant obtained with Type-1 boundary condition. We further investigate the difference in the normalized effective piezoelectric constant e 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 (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 mV. 2.5 1.5 IGA-Open circuit IGA-Close circuit 1 Mixed FEM-Open circuit Mixed FEM-Close circuit 0.5 0 1 2 3 4 5 6 7 8 Figure 8. Effective piezoelectric constant under different electrical boundary conditions. e¯ Normalized effective piezoelectric constant Energies 2020, 13, 1326 22 of 29 (a) IGA-Type-1 boundary-f (b) IGA-Type-2 boundary-f (c) Mixed FEM-Type-1 boundary-f (d) Mixed FEM-Type-2 boundary-f Figure 9. Electric potential distribution under Type-1 and Type-2 electrical boundary conditions. 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. phi phi 5.113e+00 5.113e+00 2.7161 2.7161 0.3191 0.3191 -2.078e+00 -2.078e+00 (a) Electric potential (b) Electric potential of the middle surface. Figure 10. Electric potential of a 3D bending flexoelectric beam. 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 jEj = 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. Energies 2020, 13, 1326 23 of 29 x 10 8.2 Experiment LME meshfree IGA 7.8 7.6 7.4 7.2 -5 -10 -8 -6 -4 -2 10 10 10 10 10 6.8 0 0.5 1 1.5 2 2.5 −6 y(m) x 10 (b) Curvature due to inverse flexoelectric effect (a) Electric field across the thickness Figure 11. Beam under electrical load. 4.3. 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 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 -continuity NURBS basis functions, as shown in Figure 12b as well as by mixed FEM as in Figure 12c. −4 x 10 Element edge Control Points 0 0.5 1 1.5 2 1 −3 x 10 (b) IGA discretization (a) Schematic of a truncated pyramid. (c) Mixed FEM discretization Figure 12. Schematic and IGA discretization of a truncated pyramid. Numerical results of electric potential f and compression strain # for the rigid boundary yy 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 E(V/m) 2 Energies 2020, 13, 1326 24 of 29 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  a 1 1 and a  a , respectively. The flexible boundary condition is considered, where the applied force 2 2 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. strain yy Phi 3.416e-06 2.728e-03 0.0018187 -1.7708e-5 0.00090934 -3.5415e-5 -4.971e-05 -2.330e-10 (a) IGA(E) - f (b) IGA(E) - # yy (c) Mixed FEM - f (d) Mixed FEM - # yy Figure 13. Rigid boundary condition. Phi strain Y 5.730e-03 3.416e-06 0.0038201 -2.3075e-5 0.00191 -4.6151e-5 0 -1.263e-10 -6.581e-05 (a) IGA(E) - f (b) IGA(E) - # yy Figure 14. Flexible boundary condition. Energies 2020, 13, 1326 25 of 29 phi Element edge 6.918e+00 Control Points 4.6121 2.306 1.043e-08 (b) Electric potential f (a) IGA mesh phi 6.918e+00 4.6121 2.306 1.043e-08 (c) Electric potential f Figure 15. 3D pyramid. 5. 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 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. Author Contributions: Conceptualization, X.Z., B.H.N., S.S.N., T.Q.T., N.A., and T.R.; methodology, X.Z., B.H.N., S.S.N., T.Q.T., N.A., and T.R.; software, B.H.N., S.S.N., and T.Q.T.; validation, B.H.N., S.S.N., and T.Q.T.; formal analysis, B.H.N., S.S.N., and T.Q.T.; investigation, B.H.N., S.S.N., and T.Q.T.; writing, X.Z., B.H.N., S.S.N., T.Q.T., N.A., and T.R.; visualization, B.H.N., S.S.N., and T.Q.T.; and supervision, X.Z. and T.R. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Acknowledgments: The authors extend their appreciation to the Distinguished Scientist Fellowship Program (DSFP) at King Saud University for funding this work. Conflicts of Interest: The authors declare no conflict of interest. Energies 2020, 13, 1326 26 of 29 Appendix A. Second-Order Derivative of NURBS Basis Functions The computation of second-order derivative with respect to the physical coordinate of the basis functions can be achieved by first noticing that the first-order derivative of the basis functions with respect to the parametric coordinate is computed from the Jacobian matrix 2 3 2 3 2 3 ¶N ¶y ¶x ¶z ¶N ¶x ¶x ¶x ¶x ¶x 6 7 6 7 6 7 ¶N ¶y ¶N i ¶x ¶z i 6 7 6 7 = (A1) 4 5 ¶h ¶h ¶h ¶h ¶y 4 5 4 5 ¶N ¶y ¶N i ¶x ¶z i ¶z ¶z ¶z ¶z ¶z Next, by taking the second-order derivative with respect to the parametric coordinate and re-arranging, the second-order derivative with respect to the physical coordinate is determined from 2 3 2 3 2 3 2 3 2 2 2 x, y, z, 2y z 2x z 2x y N N x y z x x x ,x ,x ,x ,x ,x ,x i,xx i,xx ,xx ,xx ,xx 6 7 2 2 2 6 7 6 7 6 7 2 3 6 x, y, z, 2y z 2x z 2x y 7 N N x y z h h ,h ,h ,h ,h ,h ,h i,yy i,hh ,hh ,hh ,hh x 6 7 6 7 6 7 6 7 6 7 6 7 6 7 i,x 2 2 2 6 7 N N x y z x, y, z, 2y z 2x z 2x y 6 i,zz7 6 i,zz7 6 ,zz ,zz ,zz76 7 6 z z z ,z ,z ,z ,z ,z ,z 7 6 7 = 6 7 6 74N 5 . (A2) i,y 6 7 6 7 6 7 6 7 N N x y z x x y y z z y z + y z x z + x z x y + x y i,yz i,hz ,hz ,hz ,hz 6 ,h ,h ,h ,h ,h ,h ,h ,h ,h7 ,z ,z ,z ,z ,z ,z ,z ,z ,z 6 7 6 7 6 7 6 7 i,z 4 5 4 5 4 5 N N x y z 4 x x y y z z y z + y z x z + x z x y + x y 5 i,xz i,xz ,zx ,zx ,zx ,z ,x ,z ,x ,z ,x ,z ,x ,x ,z ,z ,x ,x ,z ,z ,x ,x ,z N N x y z x x y y z z y z + y z x z + x z x y + x y i,xy i,xh ,xh ,xh ,xh ,h ,h ,h ,h ,h ,h ,h ,h ,h ,x ,x ,x ,x ,x ,x ,x ,x ,x References 1. Tolpygo, K. Long wavelength oscillations of diamond-type crystals including long range forces. Sov. Phys.-Solid State 1963, 4, 1297–1305. 2. Kogan, S.M. Piezoelectric effect during inhomogeneous deformation and acoustic scattering of carriers in crystals. Sov. Phys.-Solid State 1964, 5, 2069–2070. 3. Huang, W.; Yan, X.; Kwon, S.R.; Zhang, S.; Yuan, F.G.; Jiang, X. Flexoelectric strain gradient detection using Ba0. 64Sr0. 36TiO3 for sensing. Appl. Phys. Lett. 2012, 101, 252903. [CrossRef] 4. Kwon, S.R.; Huang, W.; Zhang, S.; Yuan, F.G.; Jiang, X. Study on a flexoelectric microphone using barium strontium titanate. J. Micromech. Microeng. 2016, 26, 045001. [CrossRef] 5. Merupo, V.I.; Guiffard, B.; Seveno, R.; Tabellout, M.; Kassiba, A. Flexoelectric response in soft polyurethane films and their use for large curvature sensing. J. Appl. Phys. 2017, 122, 144101. [CrossRef] 6. Bhaskar, U.K.; Banerjee, N.; Abdollahi, A.; Wang, Z.; Schlom, D.G.; Rijnders, G.; Catalan, G. A flexoelectric microelectromechanical system on silicon. Nat. Nanotechnol. 2016, 11, 263. [CrossRef] 7. Zhang, S.; Liu, K.; Xu, M.; Shen, S. A curved resonant flexoelectric actuator. Appl. Phys. Lett. 2017, 111, 082904. [CrossRef] 8. Rey, A.D.; Servio, P.; Herrera-Valencia, E. Bioinspired model of mechanical energy harvesting based on flexoelectric membranes. Phys. Rev. E 2013, 87, 022505. [CrossRef] 9. Jiang, X.; Huang, W.; Zhang, S. Flexoelectric nano-generator: Materials, structures and devices. Nano Energy 2013, 2, 1079–1092. [CrossRef] 10. Deng, Q.; Kammoun, M.; Erturk, A.; Sharma, P. Nanoscale flexoelectric energy harvesting. Int. J. Solids Struct. 2014, 51, 3218–3225. [CrossRef] 11. Choi, S.B.; Kim, G.W. Measurement of flexoelectric response in polyvinylidene fluoride films for piezoelectric vibration energy harvesters. J. Phys. D Appl. Phys. 2017, 50, 075502. [CrossRef] 12. Liang, X.; Hu, S.; Shen, S. Nanoscale mechanical energy harvesting using piezoelectricity and flexoelectricity. Smart Mater. Struct. 2017, 26, 035050. [CrossRef] 13. Zhu, R.; Wang, Z.; Ma, H.; Yuan, G.; Wang, F.; Cheng, Z.; Kimura, H. Poling-free energy harvesters based on robust self-poled ferroelectric fibers. Nano Energy 2018, 50, 97–105. [CrossRef] 14. Lu, H.; Bark, C.W.; De Los Ojos, D.E.; Alcala, J.; Eom, C.B.; Catalan, G.; Gruverman, A. Mechanical writing of ferroelectric polarization. Science 2012, 336, 59–61. [CrossRef] [PubMed] 15. Lu, H.; Liu, S.; Ye, Z.; Yasui, S.; Funakubo, H.; Rappe, A.M.; Gruverman, A. Asymmetry in mechanical polarization switching. Appl. Phys. Lett. 2017, 110, 222903. [CrossRef] Energies 2020, 13, 1326 27 of 29 16. Guo, R.; Shen, L.; Wang, H.; Lim, Z.; Lu, W.; Yang, P.; Gruverman, A.; Venkatesan, T.; Feng, Y.P.; Chen, J. Tailoring Self-Polarization of BaTiO3 Thin Films by Interface Engineering and Flexoelectric Effect. Adv. Mater. Interfaces 2016, 3, 1600737. [CrossRef] 17. Yang, M.M.; Kim, D.J.; Alexe, M. Flexo-photovoltaic effect. Science 2018, 360, 904–907. [CrossRef] 18. Liu, Y.; Chen, J.; Deng, H.; Hu, G.; Zhu, D.; Dai, N. Anomalous thermoelectricity in strained Bi 2 Te 3 films. Sci. Rep. 2016, 6, 32661. [CrossRef] 19. Nguyen, T.D.; Mao, S.; Yeh, Y.W.; Purohit, P.K.; McAlpine, M.C. Nanoscale flexoelectricity. Adv. Mater. 2013, 25, 946–974. [CrossRef] 20. Zubko, P.; Catalan, G.; Tagantsev, A.K. Flexoelectric effect in solids. Annu. Rev. Mater. Res. 2013, 43. [CrossRef] 21. Yudin, P.; Tagantsev, A. Fundamentals of flexoelectricity in solids. Nanotechnology 2013, 24, 432001. [CrossRef] [PubMed] 22. Krichen, S.; Sharma, P. Flexoelectricity: A perspective on an unusual electromechanical coupling. J. Appl. Mech. 2016, 83, 030801. [CrossRef] 23. Maranganti, R.; Sharma, N.; Sharma, P. Electromechanical coupling in nonpiezoelectric materials due to nanoscale nonlocal size effects: Green’s function solutions and embedded inclusions. Phys. Rev. B 2006, 74, 014110. [CrossRef] 24. Mindlin, R.D. Polarization gradient in elastic dielectrics. Int. J. Solids Struct. 1968, 4, 637–642. [CrossRef] 25. Sahin, E.; Dost, S. A strain-gradients theory of elastic dielectrics with spatial dispersion. Int. J. Eng. Sci. 1988, 26, 1231–1245. [CrossRef] 26. Majdoub, M.; Sharma, P.; Çagin, ˘ T. Dramatic enhancement in energy harvesting for a narrow range of dimensions in piezoelectric nanostructures. Phys. Rev. B 2008, 78, 121407. [CrossRef] 27. Majdoub, M.; Sharma, P.; Çagin, ˘ T. Erratum: Dramatic enhancement in energy harvesting for a narrow range of dimensions in piezoelectric nanostructures [Phys. Rev. B 78, 121407 (R)(2008)]. Phys. Rev. B 2009, 79, 159901. [CrossRef] 28. Majdoub, M.; Sharma, P.; Cagin, T. Enhanced size-dependent piezoelectricity and elasticity in nanostructures due to the flexoelectric effect. Phys. Rev. B 2008, 77, 125424. [CrossRef] 29. Majdoub, M.S.; Sharma, P.; Çagin, ˘ T. Erratum: Enhanced size-dependent piezoelectricity and elasticity in nanostructures due to the flexoelectric effect [Phys. Rev. B77, 125424 (2008)]. Phys. Rev. B 2009, 79. [CrossRef] 30. Sharma, N.; Landis, C.; Sharma, P. Piezoelectric thin-film superlattices without using piezoelectric materials. J. Appl. Phys. 2010, 108, 024304. [CrossRef] 31. Hu, S.; Shen, S. Variational principles and governing equations in nano-dielectrics with the flexoelectric effect. Sci. China Physics, Mech. Astron. 2010, 53, 1497–1504. [CrossRef] 32. Kuang, Z.B. Some problems in electrostrictive and magnetostrictive materials. Acta Mech. Solida Sin. 2007, 20, 219–227. [CrossRef] 33. Kuang, Z. Some variational principles in electroelastic media under finite deformation. Sci. China Ser. Physics, Mech. Astron. 2008, 51, 1390–1402. [CrossRef] 34. Kuang, Z.b. Internal energy variational principles and governing equations in electroelastic analysis. Int. J. Solids Struct. 2009, 46, 902–911. [CrossRef] 35. Kuang, Z.B. Theory Electroelasticity; Springer: Berlin/Heidelberg, Germany, 2014. 36. Liu, L. On energy formulations of electrostatics for continuum media. J. Mech. Phys. Solids 2013, 61, 968–990. [CrossRef] 37. Liu, L. An energy formulation of continuum magneto-electro-elasticity with applications. J. Mech. Phys. Solids 2014, 63, 451–480. [CrossRef] 38. Abdollahi, A.; Peco, C.; Millan, D.; Arroyo, M.; Arias, I. Computational evaluation of the flexoelectric effect in dielectric solids. J. Appl. Phys. 2014, 116, 093502. [CrossRef] 39. Mao, S.; Purohit, P.K.; Aravas, N. Mixed finite-element formulations in piezoelectricity and flexoelectricity. Proc. R. Soc. A 2016, 472, 20150879. [CrossRef] 40. Deng, F.; Deng, Q.; Yu, W.; Shen, S. Mixed Finite Elements for Flexoelectric Solids. J. Appl. Mech. 2017, 84, 081004. [CrossRef] 41. Ghasemi, H.; Park, H.S.; Rabczuk, T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Comput. Methods Appl. Mech. Eng. 2017, 313, 239–258. [CrossRef] Energies 2020, 13, 1326 28 of 29 42. Thai, T.Q.; Rabczuk, T.; Zhuang, X. A large deformation isogeometric approach for flexoelectricity and soft materials. Comput. Methods Appl. Mech. Eng. 2018, 341, 718–739. 43. Nguyen, B.; Zhuang, X.; Rabczuk, T. Numerical model for the characterization of Maxwell-Wagner relaxation in piezoelectric and flexoelectric composite material. Comput. Struct. 2018, 208, 75–91. [CrossRef] 44. Nguyen, B.; Zhuang, X.; Rabczuk, T. NURBS-based formulation for nonlinear electro-gradient elasticity in semiconductors. Comput. Methods Appl. Mech. Eng. 2019, 346, 1074–1095. 45. Yvonnet, J.; Liu, L. A numerical framework for modeling flexoelectricity and Maxwell stress in soft dielectrics at finite strains. Comput. Methods Appl. Mech. Eng. 2017, 313, 450–482. [CrossRef] 46. Codony, D.; Marco, O.; Fernández-Méndez, S.; Arias, I. An Immersed Boundary Hierarchical B-spline method for flexoelectricity. arXiv 2019, arXiv:1902.02567. 47. Roy, P.; Roy, D. Peridynamics model for flexoelectricity and damage. Appl. Math. Model. 2019, 68, 82–112. [CrossRef] 48. Toupin, R.A. The elastic dielectric. J. Ration. Mech. Anal. 1956, 5, 849–915. [CrossRef] 49. Mindlin, R.D.; Eshel, N. On first strain-gradient theories in linear elasticity. Int. J. Solids Struct. 1968, 4, 109–124. [CrossRef] 50. Abdollahi, A.; Millán, D.; Peco, C.; Arroyo, M.; Arias, I. Revisiting pyramid compression to quantify flexoelectricity: A three-dimensional simulation study. Phys. Rev. B 2015, 91, 104103. [CrossRef] 51. Abdollahi, A.; Domingo, N.; Arias, I.; Catalan, G. Converse flexoelectricity yields large piezoresponse force microscopy signals in non-piezoelectric materials. Nat. Commun. 2019, 10, 1266. [CrossRef] 52. He, B.; Javvaji, B.; Zhuang, X. Characterizing Flexoelectricity in Composite Material Using the Element-Free Galerkin Method. Energies 2019, 12, 271. [CrossRef] 53. Liu, G.R. Meshfree Methods: Moving Beyond the Finite Element Method; CRC Press: Boca Raton, FL, USA, 2009. 54. Cordes, L.; Moran, B. Treatment of material discontinuity in the element-free Galerkin method. Comput. Methods Appl. Mech. Eng. 1996, 139, 75–89. [CrossRef] 55. Krongauz, Y.; Belytschko, T. EFG approximation with discontinuous derivatives. Int. J. Numer. Methods Eng. 1998, 41, 1215–1233. [CrossRef] 56. Nanthakumar, S.; Zhuang, X.; Park, H.S.; Rabczuk, T. Topology optimization of flexoelectric structures. J. Mech. Phys. Solids 2017, 105, 217–234. [CrossRef] 57. Deng, F.; Deng, Q.; Shen, S. A three-dimensional mixed finite element for flexoelectricity. J. Appl. Mech. 2018, 85, 031009. [CrossRef] 58. Les Piegl, W.T. The NURBS Book; Springer: Berlin/Heidelberg, Germany, 1997. 59. Cottrell, J.; Hughes, T.; Bazilevs, Y. Isogeometric Analysis: Toward Integration of CAD and FEA; John Wiley & Sons Inc.: Hoboken, NJ, USA, 2009. 60. Camacho, G.T.; Ortiz, M. Computational modelling of impact damage in brittle materials. Int. J. Solids Struct. 1996, 33, 2899–2938. [CrossRef] 61. Schrefler, B.A.; Secchi, S.; Simoni, L. On adaptive refinement techniques in multi-field problems including cohesive fracture. Comput. Methods Appl. Mech. Eng. 2006, 195, 444–461. [CrossRef] 62. de Borst, R. Fluid flow in fractured and fracturing porous media: A unified view. Mech. Res. Commun. 2017, 80, 47–57. [CrossRef] 63. Sukumar, N.; Chopp, D.L.; Moës, N.; Belytschko, T. Modeling holes and inclusions by level sets in the extended finite-element method. Comput. Methods Appl. Mech. Eng. 2001, 190, 6183–6200. [CrossRef] 64. Sukumar, N.; Huang, Z.; Prévost, J.H.; Suo, Z. Partition of unity enrichment for bimaterial interface cracks. Int. J. Numer. Methods Eng. 2004, 59, 1075–1102. [CrossRef] 65. Shu, J.Y.; King, W.E.; Fleck, N.A. Finite elements for materials with strain gradient effects. Int. J. Numer. Methods Eng. 1999, 44, 373–391. [CrossRef] 66. Soh, A.K.; Wanji, C. Finite element formulations of strain gradient theory for microstructures and the C0–1 patch test. Int. J. Numer. Methods Eng. 2004, 61, 433–454. [CrossRef] 67. Gu, S.; He, Q.C. Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces. J. Mech. Phys. Solids 2011, 59, 1413–1426. [CrossRef] 68. Liu, J.T.; He, B.; Gu, S.T.; He, Q.C. Implementation of a physics-based general elastic imperfect interface model in the XFEM and LSM context. Int. J. Numer. Methods Eng. 2018, 115, 1499–1525. [CrossRef] 69. Nguyen, V.P.; Anitescu, C.; Bordas, S.P.; Rabczuk, T. Isogeometric analysis: an overview and computer implementation aspects. Math. Comput. Simul. 2015, 117, 89–116. [CrossRef] Energies 2020, 13, 1326 29 of 29 70. Fischer, P.; Klassen, M.; Mergheim, J.; Steinmann, P.; Müller, R. Isogeometric analysis of 2D gradient elasticity. Comput. Mech. 2011, 47, 325–334. [CrossRef] 71. Shu, L.; Wei, X.; Pang, T.; Yao, X.; Wang, C. Symmetry of flexoelectric coefficients in crystalline medium. J. Appl. Phys. 2011, 110, 104106. [CrossRef] 72. Le Quang, H.; He, Q.C. The number and types of all possible rotational symmetries for flexoelectric tensors. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 2011, 467, 2369–2386. 73. Hong, J.; Vanderbilt, D. First-principles theory and calculation of flexoelectricity. Phys. Rev. B 2013, 88, 174107. [CrossRef] 74. Gao, X.L.; Park, S. Variational formulation of a simplified strain gradient elasticity theory and its application to a pressurized thick-walled cylinder problem. Int. J. Solids Struct. 2007, 44, 7486–7499. [CrossRef] 75. Aravas, N. Plane-strain problems for a class of gradient elasticity models—A stress function approach. J. Elast. 2011, 104, 45–70. [CrossRef] 76. Bursian, E.; OI, Z. Changes in curvature of a ferroelectric film due to polarization. Sov. Phys. Solid State, USSR 1968, 10, 1121–1124. 77. Lumentut, M.; Shu, Y. A unified electromechanical finite element dynamic analysis of multiple segmented smart plate energy harvesters: circuit connection patterns. Acta Mech. 2018, 229, 4575–4604. [CrossRef] 78. Lumentut, M.; Howard, I. Intrinsic electromechanical dynamic equations for piezoelectric power harvesters. Acta Mech. 2017, 228, 631–650. [CrossRef] Sample Availability: Detailed data and codes are available from the authors by request. c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Energies Multidisciplinary Digital Publishing Institute

Computational Modeling of Flexoelectricity—A Review

Loading next page...
 
/lp/multidisciplinary-digital-publishing-institute/computational-modeling-of-flexoelectricity-a-review-EYT0hIuDev

References (79)

Publisher
Multidisciplinary Digital Publishing Institute
Copyright
© 1996-2020 MDPI (Basel, Switzerland) unless otherwise stated Disclaimer The statements, opinions and data contained in the journals are solely those of the individual authors and contributors and not of the publisher and the editor(s). Terms and Conditions Privacy Policy
ISSN
1996-1073
DOI
10.3390/en13061326
Publisher site
See Article on Publisher Site

Abstract

energies Review 1,2 3 3 Xiaoying Zhuang , Binh Huy Nguyen , Subbiah Srivilliputtur Nanthakumar , 3 4 4, Thai Quoc Tran , Naif Alajlan and Timon Rabczuk * Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City 758307, Vietnam; xiaoying.zhuang@tdtu.edu.vn Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City 758307, Vietnam Institute for Continuum Mechanics, Leibniz Universität Hannover, Appelstr. 11, 30167 Hannover, Germany; binh.nguyen@ikm.uni-hannover.de (B.H.N.); nanthakumar@ikm.uni-hannover.de (S.S.N.); tran@ikm.uni-hannover.de (T.Q.T.) Department of Computer Engineering, College of Computer and Information Sciences, King Saud University, Riyadh 11543, Saudi Arabia; timon.rabczuk@uni-weimar.de * Correspondence: timon.rabczuk@uni-weimar.de Received: 21 January 2020; Accepted: 26 February 2020; Published: 12 March 2020 Abstract: 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. Keywords: flexoelectricity; numerical methods; modeling 1. 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–5], actuators [6,7] and nano-energy harvesters [8–13]. It can also explain physical phenomena such as strain gradient-driven polarity control [14–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–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 Energies 2020, 13, 1326; doi:10.3390/en13061326 www.mdpi.com/journal/energies Energies 2020, 13, 1326 2 of 29 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–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–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 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–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 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. 2. Ground Work Formulation 2.1. Physical Variational Principle Hu and Shen [31] extended the variational principle from Toupin [48] by the physical variational principle proposed by Kuang [32–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 e , strain gradient e , polarization P ij ij,k i and polarization gradient P , i,j 1 1 1 1 U e , e , P , P = a P P + b P P + c # # + d P # + f # P + h P # + g # # , (1) ij ij,k i i,j ij i j ijkl i,j k,l ijkl ij kl ijk i jk ijkl ij k,l ijkl i jk,l ijklmn ij,k lm,n 2 2 2 2 Energies 2020, 13, 1326 3 of 29 from which the constitutive relations are obtained as ¶U s = = c # + d P + f P (2a) ij ijkl kl ijk k ijkl k,l ¶# ij ¶U t = = h P + g # (2b) ijk ijkl l ijklmn lm,n ¶# ij,k ¶U E = = a P + d # + h # (2c) i ij j ijk jk ijkl jk,l ¶P ¶U V = = b P + f # , (2d) ij ijkl k,l ijkl kl ¶P i,j where s , E , t and V are the (local) stress, electric field, higher-order stress and higher-order ij i ijk ij electric field, respectively, while a , b , c , d , f , and h are the material tensors [32]. In other ij ijkl ijkl ijk ijkl ijkl words, the variation of internal energy can be written as dU = s d# + t d# + E dP + V dP . (3) ij ij ijk ij,k i i ij i,j The energy density is the sum of the internal energy and the internal of the Maxwell self-field [24] W = U + # f f , (4) 0 ,i ,i in which the electric Maxwell self-field can be defined as the gradient of the electric potential f MS E = f . (5) ,i Next, we can define the electric enthalpy as H = W E D = U # f f + f P (6) i i 0 ,i ,i ,i i We consider a dielectric solid occupying a volume v bounded by a surface a that separates from the 0  0 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 d H dv = dW, (7) ext f where W is the external work done by body force f , external electric field E and free charge r on the solid body and traction t , double-traction t , surface charge s and higher-order electric field v ¯ on i i i the surface boundary. Hence, its variation dW is defined as R  R R R R R ext n  f dW = f du + E dP dv + t du da + t ¯ D (du ) da + v ¯ dP da s df da r df dv , (8) s t V f i i i i i i i i i v i a a a a v in which the operator D is explained in the following section. Within physical variational principle approach, variation of the electric enthalpy can be further expressed as [32–35]. Z Z Z d H dv = dH dv + H du dv , (9) k,k v v v e 1 where H = E P + V P is the energy deducted from the total energy minus the pure deformation i i ij i,j energy. Moreover, the migratory variation is the essential feature of physical variational principle, Energies 2020, 13, 1326 4 of 29 in the sense that the virtual displacement not only produces the strain variation but also causes electric potential and polarization variations, such as [31] df = d f + d f = d f + f du , (10a) f u f ,i i df = d f + d f = d f + f du , (10b) ,i f ,i u ,i f ,i ,ij j dP = d P + d P = d P + P du , (10c) i P i i P i i,j j dP = d P + d P = d P + P du = d P + P du , (10d) i,j P i,j u i,j P i,j i,jk k P i,j i,kj k where d () and d () are the variations produced by the (local) variation of the electric potential and polarization, respectively, whereas d  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] Z Z h i ES d H dv = s + t s du dv ij,j ijm,mj i ij,j v v nh i h io ES t t ˜ ˜ + s t + s n + D (n )n n t D n t du da m m ij ijm,m j l j ijm ijm i ij l j + t n n D (du ) da ijm m j i Z Z ext + E V + f + E dP dv + V n dP da i ij,j ,i i ij j i v a Z Z (# f + P ) df dv + (# f + P )n df da , (11) 0 ,i i 0 ,i i i ,i v a t n ˜ ˜ where n is the normal vector. D () = (d n n ) ¶ () and D = n ¶ () are the tangential and i kl k l ,l l ,l normal surface differential operators, respectively, which emerges from the treatment of the variation of displacement gradient du on the boundary (see also [31,49]). Notably, in the above equation, as a i,j ES direct result of the variational process, the generalized electrostatic stress s appears as ij 1 1 ES s = V P + (E P + V P )d D f + # f f + f P d kj k,i k k kl k,l ij j ,i ,k ,k ,k k ij ij 2 2 M M = s + t , (12) ij ijm,m M M with s and t the Maxwell stress and the electrostatic stress corresponding to the strain and ij ijm,m polarization gradients, respectively, defined as 1 1 s = E P d D f + f + # f f + f P d (13a) k k ij j ,j ,i ,k ,k ,k k ij ij 2 2 t = V P + V P . (13b) kj k,j kl k,l ijm,m Upon substituting Equations (11) and (8) into Equation (7) and utilizing the arbitrariness of du , df and dP , 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 Energies 2020, 13, 1326 5 of 29 ES s t + s + f = 0, in v (14a) ij ijm,m i ij ,j ES s = 0, in v’ (14b) ij,j L ext V + E f + E = 0, in v (14c) ij,j ,i i i # f + P = r , in v (14d) 0 ,ii i,i f = 0, in v’ (14e) ,ii Dirichlet boundary conditions u = u ¯ , on a (15a) i i f = f, on a (15b) n u,n D (u ) = u n = u ¯ n , on a (15c) i i,l l i,l l P = P , on a (15d) i i Neumann boundary conditions ES t t s ˜ ˜ s t + [[s ]] n + D (n )n n t D n t = t , on a (16a) m m ij ijm,m j l j ijm ijm i ij l j t n n = t , on a (16b) ijm j i n (# [[f ]] + P ) = s , on a (16c) i 0 ,i i V n = v ¯ , on a (16d) ij j i 2.2. Alternative Formulation The above boundary value problem is formulated from the total energy density U whose strain # , strain gradient # , polarization P and polarization gradient P are independent variables. ij ij,k i i,j 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 # , strain gradient # , electric ij ij,k field E and electric field gradient E [38,40,50]: i i,j 1 1 1 1 G # , # , E , E = k E E b E E + c # # e E # m (E # # E ) (17) ij ij,k i i,j ij i j ijkl i,j k,l ijkl ij kl ijk i jk ijkl k ij,l ij k,l 2 2 2 2 Consequently, the constitutive relations are given as ¶G 1 s = = c # + m E (18a) ij ijkl kl ijkl k,l ¶# 2 ij ¶G 1 t = = m E (18b) ijk ijkl l ¶# 2 ij,k ¶G 1 D = = k E + m # (18c) i ij j ijkl k,l ¶E 2 ¶G 1 Q = = b E m # , (18d) ij ijkl k,l ijkl kl ¶E 2 i,j where s , D , t and V are the (local) stress, electric displacement, higher-order stress and ij i ijk ij higher-order electric displacement, respectively, while k , b , c , e and m are material ij ijkl ijkl ijk ijkl tensors [38]. These relations also imply that the variation of electric Gibbs free energy can be expressed as dG = s d# + t d# D dE + Q dE , (19) ij ij ijk ij,k i i ij i,j Energies 2020, 13, 1326 6 of 29 which then can be used to determine the equilibriums from variational principle d G dv = dW, (20) 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 charges q ¯ are prescribed on the boundary, so that its variation dW reads R  R R R R R ext n n  f ˜ ˜ dW = f du + E dP dv + t du da + t ¯ D (du ) da + q ¯D (df) da s df da r df dv . (21) s t Q f i i i i i i i v i a a a a v By analogy with Equation (10), migratory variation of electric field and electric field gradient are obtained as dE = d E + d E = d f + E du = d f + E du , (22a) i f i u i f ,i i,p p f ,i p,i p dE = d E + d E = d E + E du = d E + E du (22b) f u f p f p i,j i,j i,j i,j i,jp i,j i,pj 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 Z Z M M d G dv = s + t s t du dv ij,j i ijk,kj ij,j ijm,mj v v M M t t ˜ ˜ + s n t n + s n + t n + D (n )n n t D n t du da m m ij j ijk,k j j j l j ijm ijm i ij ijm,m l j Z Z + t n n D (du ) da + D + D df dv ijk k j i i,i ij,ji a v Z Z t t n ˜ ˜ ˜ + D n D n + D (n )n n Q D (Q n ) df da + Q n n D (df) da , (23) i i ij,j i l i j ij i ij j ij i j a a M M where s and t are the Maxwell stress and the higher-order electrostatic stress that are direct ij ijm,m results of the variational process and expressed as s = E D E D d (24a) i j k k ij ij t = E Q E Q d . (24b) ij ijm,m i,k kj k,l kl 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 ES s t + s + f = 0, in v (25a) ij ijm,m i ij ,j ES s = 0, in v’ (25b) ij,j D Q = r , in v (25c) i ij,j ,i D Q = f = 0, in v’ (25d) i ij,j ,ii ,i Dirichlet boundary conditions u = u ¯ , on a (26a) i i f = f, on a (26b) n u,n D (u ) = u n = u ¯ n , on a (26c) i,l l i,l l n f,n ˜ ¯ D (f) = f n = f n , on a (26d) ,l l ,l l Energies 2020, 13, 1326 7 of 29 Neumann boundary conditions ES t t s ˜ ˜ s t + [[s ]] n + D (n )n n t D n t = t , on a (27a) m m ij ijm,m j l j ijm ijm i ij l j t n n = t ¯ , on a (27b) ijk k j i t t  s ˜ ˜ [[D Q ]]n + D (n )n n Q D (Q n ) = s , on a (27c) i ij,j i l i j ij ij j l i Q n n = q ¯, on a (27d) ij i j 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 D Q is degenerated to f in i ij,j ,i vacuum, which also causes the jump [[D Q ]] on the boundary a . i ij,j 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 re results its conjugate variable, double stress t , to be cast in the ijk 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 , to appear ij 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 enters the intramolecular force balance equation and higher-order ij 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 s , whereas the higher-order Maxwell stress are given in their corresponding electric independent ij ES variables. It is worth noting that the electrostatic stress tensor s in vacuum is reduced to ij ES M s = s = # E E E E d , (28) 0 i j k k ij ij ij which is the usual Maxwell stress tensor in electromagnetism. Furthermore, we remark that different formulations from the literature can be deduced from ES the generalized ones presented here. For instance, in the absence of electrostatic stress s and the ij polarization gradient P (in P-formulation) or the electric field gradient E (in E-formulation), we can i,j i,j 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 = f ,i or the usual relation between electric field and electric potential E = f can be retrieved and the i ,i 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 R R R R R s du + t du + (# E + P )df dv = f du dv + t du da r df dv  s df da (29) 0 s s ij i,j ijm i,jm i i ,i i i i i v v a v a for P-formulation and Z Z Z Z Z s du + t du + D df dv = f du dv + t du da r df dv s df da (30) ij i,j ijm i,jm i ,i i i i i s s v v a v a 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 1 1 1 by utilizing the constitutive relation in Equation (2c), i.e., P = a f a d # a h # . i ,j jkl kl jklm kl,m ij ij ij Energies 2020, 13, 1326 8 of 29 On the other hand, as the electric field E = f is an independent variable in E-formulation, it i ,i 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. t t t ˜ ˜ ˜ Table 1. Comparison of strong form. L = D (n )n n t D n t , M = D (n )n n Q i l m j ijm m ijm l i j ij l j l D (Q n ). ij j U(e, P,re,rP) G(e, E,re,rE) ES ES Momentum in v s t + s + f = 0 s t + s + f = 0 ij ijm,m i ij i ijk,k ij ij ,j ,j 0 ES ES Maxwell in v s = 0 s = 0 ij,j ij,j Balance f f Gauss law in v # f + P = r D Q = r 0 ,ii i,i i ij,j ,i L ext Intramolecular force [24] V + E f + E = 0 ij,j ,i i i Gauss law in v f = 0 f = 0 ,ii ,ii Displacement on a u = u ¯ u = u ¯ i i i i ¯ ¯ Electric potential on a f = f f = f u,n Dirichlet BCs Displacement normal derivative on a u n = u ¯ n u n = u ¯ n i,l l i,l l i,l l i,l l f,n Potential normal derivative on a f n = f n ,l l ,l l Surface polarization on a P = P i i ES ES Surface traction on a s t + [[s ]] n + L = t s t + [[s ]] n + L = t ij ijm,m j i i ij ijk,k j i i ij ij Surface charge on a n # [[f ]] + P = s [[D D ]]n + M = s i 0 ,i i i ij,j i Neumann BCs Higher-order traction on a t n n = t ¯ t n n = t ¯ ijm m j i ijm m j i Higher-order surface charge on a Q n n = q ¯ ij i j Higher-order electric field on a V n = v ij j i t t ˜ ˜ Table 2. Comparison of reduced strong form. L = D (n )n n t D n t , M = m m i l j ijm ijm l j t t ˜ ˜ D (n )n n Q D (Q n ). l i j ij ij j l i U(e, P,re) G(e, E,re) References [30,39] References [38,40,43] ES ES Momentum in v s t + s + f = 0 s t + s + f = 0 ij ijm,m i ij ijk,k i ij ij ,j ,j f f Gauss law in v # f + P = r D Q = r 0 ,ii i,i i ij,j Balance ,i L ext Intramolecular force [24] V + E f + E = 0 ij,j ,i i i Gauss law in v f = 0 f = 0 ,ii ,ii Displacement on a u = u ¯ u = u ¯ i i i i ¯ ¯ Dirichlet BCs Electric potential on a f = f f = f u,n Displacement normal derivative on a u n = u ¯ n u n = u ¯ n i,l l i,l l i,l l i,l l s ES ES ¯ ¯ Surface traction on a s t + [[s ]] n + L = t s t + [[s ]] n + L = t ij ijm,m j i i ij j i i ijk,k ij ij Neumann BCs s Surface charge on a n # [[f ]] + P = s [[D D ]]n + M = s i 0 ,i i i ij,j i Higher-order traction on a t n n = t ¯ t n n = t ¯ m m ijm j i ijm j i 3. 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 -continuity of the second-order derivatives of the displacement field by means of higher-order (at least C -continuity) shape functions. In the following, we summarize the numerical methods that have been employed in flexoelectricity, including meshless, mixed FE and IGA. Energies 2020, 13, 1326 9 of 29 3.1. 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 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. 3.2. Mixed Finite Element Method The mixed finite element method offers the advantage of capturing C continuity in the domain by utilizing C 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 -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. 3.3. 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], 0 1 in which the geometrical approximation across the interface is C -continuity while C -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 continuity is imposed across the void–solid phase. A more general treatment for modeling discontinuity across the Energies 2020, 13, 1326 10 of 29 internal boundary can be adopted from the works in [60–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 -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. 1 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ξ ξ (a) (b) Figure 1. (a) C continuity quadratic NURBS basis functions from knot vector X = f0, 0, 0, 1/3, 1/3, 2/3, 2/3, 1, 1, 1g; and (b) C continuity quadratic NURBS basis functions from knot vector X = 0, 0, 0, 1/3, 2/3, 1, 1, 1 f g Perfect interfaces result in continuous displacement and traction fields along the interface, i.e., [[u ]] = 0 (31a) [[t ]] = 0, (31b) 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 n ]] = 0 (32a) i,l l [[t n n ]] = 0, (32b) ijm m j 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 R(ξ) R(ξ) Energies 2020, 13, 1326 11 of 29 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. 3.4. 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 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 continuity displacement gradients, y are introduced in the formulation and the kinematic constraint, y = ru, is imposed using Lagrange multipliers. The weak form of mechanical and electrostatic equilibrium u 1 u,n 1 is obtained by finding the u 2 u = u ¯ on a , u 2 H (v) , y 2 y = y on a , y 2 H (v) , f 2 f 1 2 u 1 f = f on a , f 2 H (v) , l 2 L (v), and similarly ffiu 2 ffiu = 0 on a , ffiu 2 H (v) , dy 2 u,n 1 f 1 2 dy = 0 on a , dy 2 H (v) , df 2 df = 0 on a , df 2 H (v) and dl 2 L (v), such that, Z Z dP = (s : #(du) + t ˆ .h ˆ(dy) D E(df)) dW + l : (dyrdu)dW W W Z Z Z (33) + (yru) : dldW du t dG du b dW = 0 W G W Z Z Z #(du) : C : #(u) dW #(du) : e E(f) dW h ˆ(dy) . h E(f) dW W W W Z Z Z . . . . . . + h ˆ(dy) . g . h ˆ(y)dW E(df) e : #(u) dW E(df) h . h ˆ(y) dW W W W (34) Z Z Z Z E(df) k E(f) dW + dy : l dW rdu : l dW dl : ru dW W W W W Z Z Z + dl : y dW = du t dG + du b dW W G W ¶U where t ˆ = and h ˆ is the third-order tensor related to relaxed displacement gradient y as h ˆ = ¶h ˆ (y + y ). The constraint y = ru is imposed in a weighted residual manner by including jk,i ik,j d (yru) : ldW in Equation (33). W Energies 2020, 13, 1326 12 of 29 Table 3. Comparison of different numerical methods. (a) and (b) indicate quadratic elements in 2D and 3D, respectively. Meshfree IGA Mixed FEM Continuity Depends on kernel functions Knot manipulation Nodal displacement gradient D.o.f Depends on domain of influence 27 (a) or 108 (b) 54 (Q54, [56]); 37 (T37), 47 (Q47), 45 (T45), 59 (Q59) [40]; 233 (3D) [57] Higher-order Dirichlet BCs - Lagrange-multiplier Direct imposition C Interface modeling Lagrange-multiplier, Enriched functions Multi-patch - Energies 2020, 13, 1326 13 of 29 A nine-noded isoparametric element was proposed by Nanthakumar et al. [56]. The degrees of freedoms are u and u at all nodes, relaxed displacement gradients y , y , y and y at the four 1 2 11 12 21 22 corner nodes and electric potential f at the four corner nodes. In addition to these DOFs, the element has four Lagrange multipliers l , l , l and l at the four corner nodes. The flexoelectric element 11 12 21 22 is shown in Figure 2. The biquadratic shape functions, N , are used to interpolate the displacement DOFs u and u . Bilinear shape functions, N , are used to interpolate the relaxed displacement 1 2 gradients, Lagrange multipliers and electric potential, f. The finite element approximation is as shown below, h h l u (x) = N u ; y (x) = N y ; i i å å i i=1 i=1,3,5,7 (35) h l h l f (x) = N f ; l (x) = N l å i å i i i i=1,3,5,7 i=1,3,5,7 h h l du (x) = N du ; dy (x) = N dy ; å i å i i i i=1 i=1,3,5,7 (36) h l h l df (x) = N df ; dl (x) = N dl å i å i i i i=1,3,5,7 i=1,3,5,7 7 6 8 4 1 2 Figure 2. A nine-noded flexoelectric element. DOFs at circles are u and u ; at squares are y , y , y , 1 2 11 12 21 y and f; and at triangles are l , l , l , and l . 22 11 12 21 22 Substituting Equations (35) and (36), into Equation (34) yields the following expression in terms of nodal DOFs: 0 1 Z Z Z T T T T T T T T @ A du B CB dW u + du B e B dW f + dy H h B dW f u f f u u W W 0 1 Z Z Z T T T T T T @ A + dy H g HdW y + df B eB dW u + df B h H dW y f f W W (37) Z Z Z T T T T l T l l df B kB dW f du B N dW l + dy N N dW l f yu W W W Z Z Z Z T T T l T l l T q T q dl N B dW u + dl N N dW y = du N t dG + du N b dW yu W W G W N Energies 2020, 13, 1326 14 of 29 2 3 ¶N 0 0 0 ¶x 6 7 ¶N 6 7 0 0 0 2 3 6 ¶x 7 ¶N " # l l 0 6 7 1 ¶N 1 ¶N ¶x ¶N 6 0 0 7 6 ¶N 7 2 ¶x 2 ¶x ¶x 6 7 where B = , B = , H = , B = u 4 5 f l yu ¶y ¶N 6 7 ¶N 0 0 0 q q 6 ¶y 7 ¶N ¶N ¶y 6 7 ¶y ¶x ¶N 6 7 0 0 0 ¶y 4 5 l l 1 ¶N 1 ¶N 0 0 2 ¶y 2 ¶y 2 3 ¶N ¶x 6 l 7 ¶N 6 0 7 ¶x 6 7 ¶N 6 7 4 5 ¶y ¶N ¶y The algebraic equations to be solved to obtain displacement and electric potential in a flexoelectric structure are, 2 32 3 2 3 K 0 K K u F uu uf u ul 6 76 7 6 7 0 K K K y 0 6 76 7 6 7 yy yf yl 6 76 7 = 6 7 (38) 4 54 5 4 5 K K K 0 f 0 fu fy ff K K 0 0 l 0 lu ly K = B CB dW ; uu u T T T K = B e B dW = K ; uf f u fu T T K = B N dW = K ul yu l lu T T K = H hB dW = K yf f fy K = B gB dW yy y l T K = N N dW = K yl l ly K = B kB dW ; ff f Z Z qT q F = N tdG + N b dW G W 3.5. IGA Approach Within the Galerkin method, both the trial and test functions are approximated by NURBS basis functions (e) (e) u = N u ˜, du = N du ˜, (39a) u u (e) (e) ˜ ˜ f = N f, df = N df, (39b) f f in which N and N are the matrices containing the local NURBS basis functions of element e defined as u f 2 3 N 0 0 i h i 6 7 1 2 n N = 0 N 0 , N = N N . . . N , (40a) 4 5 i u u u u u 0 0 N h i N = N N . . . N , (40b) 1 2 n Energies 2020, 13, 1326 15 of 29 where N is the NURBS basis function at control point ith and n is the number of control points of an element. In the above equation, u ˜ , f, du ˜ and du ˜ 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 h i 1 1 1 n n n u ˜ = u ˜ u ˜ u ˜ . . . u ˜ u ˜ u ˜ , (41) x y z x y z h i 1 1 1 n n n du ˜ = du ˜ du ˜ du ˜ . . . du ˜ du ˜ du ˜ , (42) x y z x y z h i 1 2 n ˜ ˜ ˜ ˜ f = f f . . . f , (43) h i 1 2 n ˜ ˜ ˜ ˜ df = df df . . . df (44) T T The approximation of strain field # = [# # # 2# 2# 2# ] , electric field E = [E E E ] , strain 11 22 33 23 13 12 1 2 3 h i h i T T ¶# ¶# ¶# ¶E ¶E ¶E gradient r# = and electric field gradient rE = can be obtained by using ¶x ¶y ¶z ¶x ¶y ¶z differential operators B ,B ,H and H on Equation (39) such that u f u f (e) (e) # = B u ˜, d# = B du ˜, (45a) u u (e) (e) ˜ ˜ E = B f, dE = B df, (45b) f f (e) (e) r# = H u ˜, dr# = H du ˜, (45c) u u (e) (e) ˜ ˜ rE = H f, drE = H df, (45d) f f in which B , B , H and H contain the first- and second-order partial derivative of the NURBS basis u f u f functions, respectively, defined as 2 3 ¶N 0 0 ¶x 6 7 ¶N i 2 3 6 0 0 7 1 2 n ¶B ¶B ¶B ¶y u u u 6 7 . . . ¶N 6 7 h i ¶x ¶x ¶x 6 7 0 0 1 2 n 6 7 i ¶B ¶B ¶B ¶z 1 2 n 6 u u u 7 B = , B = , H = (46) 6 7 B B . . . B . . . u u u ¶N ¶N u u u i i 4 5 ¶y ¶y ¶y 6 0 7 ¶z ¶y 1 2 n 6 7 ¶B ¶B ¶B u u u 6 ¶N ¶N 7 . . . i i ¶z ¶z ¶z 4 5 ¶z ¶x ¶N ¶N i i ¶y ¶x 2 3 1 2 n 2 3 ¶B ¶B ¶B f f f ¶N . . . ¶x ¶x ¶x ¶x h i 6 7 1 2 n 6 ¶N 7 6 ¶B ¶B ¶B 7 i 1 2 n i f f f B = , B = B B . . . B , H = (47) 4 5 6 7 f f . . . f f f f ¶y ¶y ¶y ¶y 4 5 ¶N 1 2 n ¶B ¶B ¶B f f f ¶z . . . ¶z ¶z ¶z Note that matrices H and H contain second-order derivative with respect to the physical coordinate u f 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 h i R  R R R R n T T e T T T T T f ¯ ˜ ˜ de s + dre t dE D dv = du f dv + du ˜ t da + f r dv df s ¯ da . (48) (e) s,(e) f,(e) e=1 v v a v a 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 " #" # " # K K u ˜ uu uf = , (49) K K f fu ff f Energies 2020, 13, 1326 16 of 29 where the stiffness matrices and force vectors are computed as n Z K = B cB dv , (50a) uu u (e) e=1 n Z T T K = B eB H mB dv , (50b) uf f f u u (e) e=1 T T T T K = B e B B m H dv , (50c) fu u u f f (e) e=1 n Z K = B kB dv , (50d) ff f (e) e=1 n Z Z T T f = N f dv + N t da , (50e) u u u s,(e) v a e=1 Z Z T f T f = N r dv + N t da . (50f) u f f,(e) v a e=1 Note that the flexoelectric tensor m for cubic symmetry crystal k is given as 2 3 m m m 0 0 0 0 0 0 0 0 m 0 0 0 0 m 0 11 12 12 44 44 6 7 m = (51) 4 0 0 0 0 0 m m m m 0 0 0 0 0 0 m 0 05 44 12 11 12 44 0 0 0 0 m 0 0 0 0 m 0 0 m m m 0 0 0 44 44 12 12 11 The symmetry of the flexoelectric tensor can be found in [71,72]. For cubic crystal, the non-zero components of flexoelectric tensor are m , m and m while other coefficients with cycling 1111 1221 1122 indices are equal (e.g., m = m ) and with odd number of indices are zeros (e.g., m = 0) [73]. 1221 1331 1112 In matrix notation, the coefficients m , m and m are denoted as m , m and m , respectively [38]. iiii jiij iijj 11 12 44 4. Numerical Examples 4.1. 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 = 0 and x = 0. 1 2 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 1 1 1 ˆ ˆ U = l# # + m# # + l l# # + 2m# # + f # P + 2 f # P + aP P (52) jj kk jk jk jj,i kk,i jk,i jk,i 12 jj,i i 44 ij,i j i i 2 2 2 and 1 1 1 G = l# # + m# # + l l# # + 2m# # m # E 2m # E kE E , (53) jj kk jk jk jj,i kk,i jk,i jk,i 12 jj,i i 44 ij,i j i i 2 2 2 where l and m are the two Lamé parameters, l is the length scale, m and m are two non-zero 12 44 ˆ ˆ flexoelectric coefficients, and f and f are two flexoelectric coupling coefficients that related to the 1 2 ˆ ˆ flexoelectric coefficients by the susceptibility c = k/# 1 such that m = c f , m = c f , a and k are 0 12 12 44 44 the reciprocal susceptibility and dielectric permittivity, also related to one another by the susceptibility. These material parameters can be found in Table 4. Furthermore, due to the geometry and boundary Energies 2020, 13, 1326 17 of 29 condition, the displacement as well as the electric potential depend only on the radial r of the hollow cylinder such as 1 r r u (r) = c r + c + c I + c K , (54a) r 1 2 3 1 4 1 r l l 0 0 m du u (r) r r f(r) = c ln r + c + + , (54b) ( ) 5 6 k dr r 2 2 where l = l + , I and K are the modified Bessel functions of the first and second kind, 1 1 (l+m)k respectively. The coefficients c , c , c , c , c , c are determined from the boundary conditions 1 2 3 4 5 6 u = u , u = u , (55a) r r r r r=r r=r i o f = f , f = f , (55b) r r r r r=r r=r i o t = t = 0 (55c) rrr rrr r=r r=r i o 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 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 = 0, (56a) x =0 u = 0, (56b) 1,2 x =0 u = 0, (56c) 2,1 x =0 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 direction u ) of the two bottom layer adjacent-bottom bottom of control points are identical, i.e., u = u . This can be implemented by utilizing x x 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. Y n l m m k 12 44 9 6 6 6 9 139 10 Pa 0.3 2 10 m 1 10 C/m 1 10 C/m 1 10 F/m Energies 2020, 13, 1326 18 of 29 −5 x 10 Element edge Control Points 1.8 1.6 1.4 1.2 0.8 0.6 0.4 0.2 0 0.5 1 1.5 2 1 −5 x 10 (b) IGA discretization (a) Schematic of a hollow cylinder. Figure 3. Schematic of an infinitely long hollow cylinder, where the inner and outer radii are r = 10 m and r = 20 m, respectively. The inner surface is imposed by radial displacement u = 0.045  m 0 r and grounded, i.e., f = 0 V, whereas the radial displacement u = 0.05 m and electric potential r r i o f = 1 V are applied on the outer surface. Symmetry boundary conditions are enforced, such that u = u = u = 0 and u = u = u = 0. 1 1,2 2,1 1,2 2,1 x =0 x =0 x =0 x =0 x =0 x =0 1 1 1 2 2 2 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 , electric potential f, radial strain # and circumferential strain # along the radial rr R qq direction with angle q = 45 are presented and compared with the analytical solutions in Figure 5. Excellent agreement between the numerical and analytical results can be seen. 4.2. 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 f 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 ) 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 given as the ratio between the electrostatic energy E and eff elec the elastic energy E elas E k E dv elec 2 v k = = R , (57) eff e C e dv elas 2 v Consequently, the normalized piezoelectric constant is defined as eff e ¯ = , (58) piezo where k is obtained as k without considering the flexoelectric effect. Physically, the normalized piezo eff piezoelectric constant indicates the enhancement of flexoelectric effect on the energy conversion. 2 Energies 2020, 13, 1326 19 of 29 Figure 7 demonstrates the normalized piezoelectric constant e ¯ as the normalized thickness h = e h/m varies under Type-1 configuration. The numerical results in the simplified case, where only 31 12 k and m are non-zeros to mimic the 1D problem, are compared to the analytical solution where 33 12 e +12(m /h) k = and excellent agreement is shown, which illustrates the great enhancement of eff k Y 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]. u r Phi 5.000e-08 1.000e+00 0.66667 4.8333e-8 0.33333 4.6667e-8 4.500e-08 0.000e+00 (a) IGA(E): Distribution of u (b) IGA(E): Distribution of f -8 Radial displacement x10 [m] Electric potential [V] 4.4674 4.5739 4.6804 4.7869 4.8935 - 0.2971 - 0.0375 0.2220 0.4815 0.7410 1.000 5.000 (d) IGA(P): Distribution of f (c) IGA(P): Distribution of u 0.8 0.6 0.4 0.2 46 0 -0.2 (e) Mixed FEM: Distribution of u (f) Mixed FEM: Distribution of f Figure 4. Distribution of radial displacement and electric potential. Energies 2020, 13, 1326 20 of 29 −8 x 10 5.1 Analytical Analytical IGA−E formulation 0.8 IGA−E formulation IGA− P formulation IGA− P formulation Mixed FEM 4.9 Mixed FEM 0.6 4.8 0.4 4.7 0.2 4.6 4.5 −0.2 4.4 −0.4 1 1.2 1.4 1.6 1.8 2 2.2 1 1.2 1.4 1.6 1.8 2 2.2 r −5 r −5 x 10 x 10 −4 −3 x 10 x 10 Analytical IGA−E formulation 4.5 IGA− P formulation Meshfree 3.5 2 Analytical IGA−E formulation IGA− P formulation Meshfree 2.5 −2 −4 2 1 1.2 1.4 1.6 1.8 2 2.2 1 1.2 1.4 1.6 1.8 2 2.2 r r −5 −5 x 10 x 10 Figure 5. Comparison of the exact and numerical results for radial displacement u , electric potential f, radial strain # and circumferential strain # as functions of radius r. rr ff Element edge Control Points -6 0 0.5 1 1.5 -5 (b) x V (a) Figure 6. (a) Schematic of a hollow cylinder. Two different types of electrical boundary condition are imposed. The upper figure describes Type-1 boundary where zero electric potential is imposed on the right edge. Type-2 boundary is demonstrated in the lower figure where zero electric potential is applied on the top edge and an unknown equipotential (due to the electrode) on the bottom edge. (b) IGA discretization (top) and mixed FEM discretization (bottom). Table 5. Material parameters for the problem depicted in Figure 6. Y n e m k k F 31 12 11 33 9 2 6 9 9 6 100 10 Pa 0.37 -4.4 C/m 1 10 C/m 1110 F/m 12.4810 F/m 10010 N ǫ u rr r φφ φ Energies 2020, 13, 1326 21 of 29 Analytical 1D (Flexo+Piezo) IGA−1D (Flexo+Piezo) Analytical 1D (Flexo) IGA−1D (Flexo) IGA− 2D (Flexo) Mixed FEM (Flexo+Piezo) Mixed FEM (Flexo) 0 1 2 3 4 5 6 7 8 Normalized thickness Figure 7. Effective normalized piezoelectric constant obtained with Type-1 boundary condition. We further investigate the difference in the normalized effective piezoelectric constant e 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 (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 mV. 2.5 1.5 IGA-Open circuit IGA-Close circuit 1 Mixed FEM-Open circuit Mixed FEM-Close circuit 0.5 0 1 2 3 4 5 6 7 8 Figure 8. Effective piezoelectric constant under different electrical boundary conditions. e¯ Normalized effective piezoelectric constant Energies 2020, 13, 1326 22 of 29 (a) IGA-Type-1 boundary-f (b) IGA-Type-2 boundary-f (c) Mixed FEM-Type-1 boundary-f (d) Mixed FEM-Type-2 boundary-f Figure 9. Electric potential distribution under Type-1 and Type-2 electrical boundary conditions. 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. phi phi 5.113e+00 5.113e+00 2.7161 2.7161 0.3191 0.3191 -2.078e+00 -2.078e+00 (a) Electric potential (b) Electric potential of the middle surface. Figure 10. Electric potential of a 3D bending flexoelectric beam. 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 jEj = 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. Energies 2020, 13, 1326 23 of 29 x 10 8.2 Experiment LME meshfree IGA 7.8 7.6 7.4 7.2 -5 -10 -8 -6 -4 -2 10 10 10 10 10 6.8 0 0.5 1 1.5 2 2.5 −6 y(m) x 10 (b) Curvature due to inverse flexoelectric effect (a) Electric field across the thickness Figure 11. Beam under electrical load. 4.3. 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 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 -continuity NURBS basis functions, as shown in Figure 12b as well as by mixed FEM as in Figure 12c. −4 x 10 Element edge Control Points 0 0.5 1 1.5 2 1 −3 x 10 (b) IGA discretization (a) Schematic of a truncated pyramid. (c) Mixed FEM discretization Figure 12. Schematic and IGA discretization of a truncated pyramid. Numerical results of electric potential f and compression strain # for the rigid boundary yy 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 E(V/m) 2 Energies 2020, 13, 1326 24 of 29 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  a 1 1 and a  a , respectively. The flexible boundary condition is considered, where the applied force 2 2 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. strain yy Phi 3.416e-06 2.728e-03 0.0018187 -1.7708e-5 0.00090934 -3.5415e-5 -4.971e-05 -2.330e-10 (a) IGA(E) - f (b) IGA(E) - # yy (c) Mixed FEM - f (d) Mixed FEM - # yy Figure 13. Rigid boundary condition. Phi strain Y 5.730e-03 3.416e-06 0.0038201 -2.3075e-5 0.00191 -4.6151e-5 0 -1.263e-10 -6.581e-05 (a) IGA(E) - f (b) IGA(E) - # yy Figure 14. Flexible boundary condition. Energies 2020, 13, 1326 25 of 29 phi Element edge 6.918e+00 Control Points 4.6121 2.306 1.043e-08 (b) Electric potential f (a) IGA mesh phi 6.918e+00 4.6121 2.306 1.043e-08 (c) Electric potential f Figure 15. 3D pyramid. 5. 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 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. Author Contributions: Conceptualization, X.Z., B.H.N., S.S.N., T.Q.T., N.A., and T.R.; methodology, X.Z., B.H.N., S.S.N., T.Q.T., N.A., and T.R.; software, B.H.N., S.S.N., and T.Q.T.; validation, B.H.N., S.S.N., and T.Q.T.; formal analysis, B.H.N., S.S.N., and T.Q.T.; investigation, B.H.N., S.S.N., and T.Q.T.; writing, X.Z., B.H.N., S.S.N., T.Q.T., N.A., and T.R.; visualization, B.H.N., S.S.N., and T.Q.T.; and supervision, X.Z. and T.R. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Acknowledgments: The authors extend their appreciation to the Distinguished Scientist Fellowship Program (DSFP) at King Saud University for funding this work. Conflicts of Interest: The authors declare no conflict of interest. Energies 2020, 13, 1326 26 of 29 Appendix A. Second-Order Derivative of NURBS Basis Functions The computation of second-order derivative with respect to the physical coordinate of the basis functions can be achieved by first noticing that the first-order derivative of the basis functions with respect to the parametric coordinate is computed from the Jacobian matrix 2 3 2 3 2 3 ¶N ¶y ¶x ¶z ¶N ¶x ¶x ¶x ¶x ¶x 6 7 6 7 6 7 ¶N ¶y ¶N i ¶x ¶z i 6 7 6 7 = (A1) 4 5 ¶h ¶h ¶h ¶h ¶y 4 5 4 5 ¶N ¶y ¶N i ¶x ¶z i ¶z ¶z ¶z ¶z ¶z Next, by taking the second-order derivative with respect to the parametric coordinate and re-arranging, the second-order derivative with respect to the physical coordinate is determined from 2 3 2 3 2 3 2 3 2 2 2 x, y, z, 2y z 2x z 2x y N N x y z x x x ,x ,x ,x ,x ,x ,x i,xx i,xx ,xx ,xx ,xx 6 7 2 2 2 6 7 6 7 6 7 2 3 6 x, y, z, 2y z 2x z 2x y 7 N N x y z h h ,h ,h ,h ,h ,h ,h i,yy i,hh ,hh ,hh ,hh x 6 7 6 7 6 7 6 7 6 7 6 7 6 7 i,x 2 2 2 6 7 N N x y z x, y, z, 2y z 2x z 2x y 6 i,zz7 6 i,zz7 6 ,zz ,zz ,zz76 7 6 z z z ,z ,z ,z ,z ,z ,z 7 6 7 = 6 7 6 74N 5 . (A2) i,y 6 7 6 7 6 7 6 7 N N x y z x x y y z z y z + y z x z + x z x y + x y i,yz i,hz ,hz ,hz ,hz 6 ,h ,h ,h ,h ,h ,h ,h ,h ,h7 ,z ,z ,z ,z ,z ,z ,z ,z ,z 6 7 6 7 6 7 6 7 i,z 4 5 4 5 4 5 N N x y z 4 x x y y z z y z + y z x z + x z x y + x y 5 i,xz i,xz ,zx ,zx ,zx ,z ,x ,z ,x ,z ,x ,z ,x ,x ,z ,z ,x ,x ,z ,z ,x ,x ,z N N x y z x x y y z z y z + y z x z + x z x y + x y i,xy i,xh ,xh ,xh ,xh ,h ,h ,h ,h ,h ,h ,h ,h ,h ,x ,x ,x ,x ,x ,x ,x ,x ,x References 1. Tolpygo, K. Long wavelength oscillations of diamond-type crystals including long range forces. Sov. Phys.-Solid State 1963, 4, 1297–1305. 2. Kogan, S.M. Piezoelectric effect during inhomogeneous deformation and acoustic scattering of carriers in crystals. Sov. Phys.-Solid State 1964, 5, 2069–2070. 3. Huang, W.; Yan, X.; Kwon, S.R.; Zhang, S.; Yuan, F.G.; Jiang, X. Flexoelectric strain gradient detection using Ba0. 64Sr0. 36TiO3 for sensing. Appl. Phys. Lett. 2012, 101, 252903. [CrossRef] 4. Kwon, S.R.; Huang, W.; Zhang, S.; Yuan, F.G.; Jiang, X. Study on a flexoelectric microphone using barium strontium titanate. J. Micromech. Microeng. 2016, 26, 045001. [CrossRef] 5. Merupo, V.I.; Guiffard, B.; Seveno, R.; Tabellout, M.; Kassiba, A. Flexoelectric response in soft polyurethane films and their use for large curvature sensing. J. Appl. Phys. 2017, 122, 144101. [CrossRef] 6. Bhaskar, U.K.; Banerjee, N.; Abdollahi, A.; Wang, Z.; Schlom, D.G.; Rijnders, G.; Catalan, G. A flexoelectric microelectromechanical system on silicon. Nat. Nanotechnol. 2016, 11, 263. [CrossRef] 7. Zhang, S.; Liu, K.; Xu, M.; Shen, S. A curved resonant flexoelectric actuator. Appl. Phys. Lett. 2017, 111, 082904. [CrossRef] 8. Rey, A.D.; Servio, P.; Herrera-Valencia, E. Bioinspired model of mechanical energy harvesting based on flexoelectric membranes. Phys. Rev. E 2013, 87, 022505. [CrossRef] 9. Jiang, X.; Huang, W.; Zhang, S. Flexoelectric nano-generator: Materials, structures and devices. Nano Energy 2013, 2, 1079–1092. [CrossRef] 10. Deng, Q.; Kammoun, M.; Erturk, A.; Sharma, P. Nanoscale flexoelectric energy harvesting. Int. J. Solids Struct. 2014, 51, 3218–3225. [CrossRef] 11. Choi, S.B.; Kim, G.W. Measurement of flexoelectric response in polyvinylidene fluoride films for piezoelectric vibration energy harvesters. J. Phys. D Appl. Phys. 2017, 50, 075502. [CrossRef] 12. Liang, X.; Hu, S.; Shen, S. Nanoscale mechanical energy harvesting using piezoelectricity and flexoelectricity. Smart Mater. Struct. 2017, 26, 035050. [CrossRef] 13. Zhu, R.; Wang, Z.; Ma, H.; Yuan, G.; Wang, F.; Cheng, Z.; Kimura, H. Poling-free energy harvesters based on robust self-poled ferroelectric fibers. Nano Energy 2018, 50, 97–105. [CrossRef] 14. Lu, H.; Bark, C.W.; De Los Ojos, D.E.; Alcala, J.; Eom, C.B.; Catalan, G.; Gruverman, A. Mechanical writing of ferroelectric polarization. Science 2012, 336, 59–61. [CrossRef] [PubMed] 15. Lu, H.; Liu, S.; Ye, Z.; Yasui, S.; Funakubo, H.; Rappe, A.M.; Gruverman, A. Asymmetry in mechanical polarization switching. Appl. Phys. Lett. 2017, 110, 222903. [CrossRef] Energies 2020, 13, 1326 27 of 29 16. Guo, R.; Shen, L.; Wang, H.; Lim, Z.; Lu, W.; Yang, P.; Gruverman, A.; Venkatesan, T.; Feng, Y.P.; Chen, J. Tailoring Self-Polarization of BaTiO3 Thin Films by Interface Engineering and Flexoelectric Effect. Adv. Mater. Interfaces 2016, 3, 1600737. [CrossRef] 17. Yang, M.M.; Kim, D.J.; Alexe, M. Flexo-photovoltaic effect. Science 2018, 360, 904–907. [CrossRef] 18. Liu, Y.; Chen, J.; Deng, H.; Hu, G.; Zhu, D.; Dai, N. Anomalous thermoelectricity in strained Bi 2 Te 3 films. Sci. Rep. 2016, 6, 32661. [CrossRef] 19. Nguyen, T.D.; Mao, S.; Yeh, Y.W.; Purohit, P.K.; McAlpine, M.C. Nanoscale flexoelectricity. Adv. Mater. 2013, 25, 946–974. [CrossRef] 20. Zubko, P.; Catalan, G.; Tagantsev, A.K. Flexoelectric effect in solids. Annu. Rev. Mater. Res. 2013, 43. [CrossRef] 21. Yudin, P.; Tagantsev, A. Fundamentals of flexoelectricity in solids. Nanotechnology 2013, 24, 432001. [CrossRef] [PubMed] 22. Krichen, S.; Sharma, P. Flexoelectricity: A perspective on an unusual electromechanical coupling. J. Appl. Mech. 2016, 83, 030801. [CrossRef] 23. Maranganti, R.; Sharma, N.; Sharma, P. Electromechanical coupling in nonpiezoelectric materials due to nanoscale nonlocal size effects: Green’s function solutions and embedded inclusions. Phys. Rev. B 2006, 74, 014110. [CrossRef] 24. Mindlin, R.D. Polarization gradient in elastic dielectrics. Int. J. Solids Struct. 1968, 4, 637–642. [CrossRef] 25. Sahin, E.; Dost, S. A strain-gradients theory of elastic dielectrics with spatial dispersion. Int. J. Eng. Sci. 1988, 26, 1231–1245. [CrossRef] 26. Majdoub, M.; Sharma, P.; Çagin, ˘ T. Dramatic enhancement in energy harvesting for a narrow range of dimensions in piezoelectric nanostructures. Phys. Rev. B 2008, 78, 121407. [CrossRef] 27. Majdoub, M.; Sharma, P.; Çagin, ˘ T. Erratum: Dramatic enhancement in energy harvesting for a narrow range of dimensions in piezoelectric nanostructures [Phys. Rev. B 78, 121407 (R)(2008)]. Phys. Rev. B 2009, 79, 159901. [CrossRef] 28. Majdoub, M.; Sharma, P.; Cagin, T. Enhanced size-dependent piezoelectricity and elasticity in nanostructures due to the flexoelectric effect. Phys. Rev. B 2008, 77, 125424. [CrossRef] 29. Majdoub, M.S.; Sharma, P.; Çagin, ˘ T. Erratum: Enhanced size-dependent piezoelectricity and elasticity in nanostructures due to the flexoelectric effect [Phys. Rev. B77, 125424 (2008)]. Phys. Rev. B 2009, 79. [CrossRef] 30. Sharma, N.; Landis, C.; Sharma, P. Piezoelectric thin-film superlattices without using piezoelectric materials. J. Appl. Phys. 2010, 108, 024304. [CrossRef] 31. Hu, S.; Shen, S. Variational principles and governing equations in nano-dielectrics with the flexoelectric effect. Sci. China Physics, Mech. Astron. 2010, 53, 1497–1504. [CrossRef] 32. Kuang, Z.B. Some problems in electrostrictive and magnetostrictive materials. Acta Mech. Solida Sin. 2007, 20, 219–227. [CrossRef] 33. Kuang, Z. Some variational principles in electroelastic media under finite deformation. Sci. China Ser. Physics, Mech. Astron. 2008, 51, 1390–1402. [CrossRef] 34. Kuang, Z.b. Internal energy variational principles and governing equations in electroelastic analysis. Int. J. Solids Struct. 2009, 46, 902–911. [CrossRef] 35. Kuang, Z.B. Theory Electroelasticity; Springer: Berlin/Heidelberg, Germany, 2014. 36. Liu, L. On energy formulations of electrostatics for continuum media. J. Mech. Phys. Solids 2013, 61, 968–990. [CrossRef] 37. Liu, L. An energy formulation of continuum magneto-electro-elasticity with applications. J. Mech. Phys. Solids 2014, 63, 451–480. [CrossRef] 38. Abdollahi, A.; Peco, C.; Millan, D.; Arroyo, M.; Arias, I. Computational evaluation of the flexoelectric effect in dielectric solids. J. Appl. Phys. 2014, 116, 093502. [CrossRef] 39. Mao, S.; Purohit, P.K.; Aravas, N. Mixed finite-element formulations in piezoelectricity and flexoelectricity. Proc. R. Soc. A 2016, 472, 20150879. [CrossRef] 40. Deng, F.; Deng, Q.; Yu, W.; Shen, S. Mixed Finite Elements for Flexoelectric Solids. J. Appl. Mech. 2017, 84, 081004. [CrossRef] 41. Ghasemi, H.; Park, H.S.; Rabczuk, T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Comput. Methods Appl. Mech. Eng. 2017, 313, 239–258. [CrossRef] Energies 2020, 13, 1326 28 of 29 42. Thai, T.Q.; Rabczuk, T.; Zhuang, X. A large deformation isogeometric approach for flexoelectricity and soft materials. Comput. Methods Appl. Mech. Eng. 2018, 341, 718–739. 43. Nguyen, B.; Zhuang, X.; Rabczuk, T. Numerical model for the characterization of Maxwell-Wagner relaxation in piezoelectric and flexoelectric composite material. Comput. Struct. 2018, 208, 75–91. [CrossRef] 44. Nguyen, B.; Zhuang, X.; Rabczuk, T. NURBS-based formulation for nonlinear electro-gradient elasticity in semiconductors. Comput. Methods Appl. Mech. Eng. 2019, 346, 1074–1095. 45. Yvonnet, J.; Liu, L. A numerical framework for modeling flexoelectricity and Maxwell stress in soft dielectrics at finite strains. Comput. Methods Appl. Mech. Eng. 2017, 313, 450–482. [CrossRef] 46. Codony, D.; Marco, O.; Fernández-Méndez, S.; Arias, I. An Immersed Boundary Hierarchical B-spline method for flexoelectricity. arXiv 2019, arXiv:1902.02567. 47. Roy, P.; Roy, D. Peridynamics model for flexoelectricity and damage. Appl. Math. Model. 2019, 68, 82–112. [CrossRef] 48. Toupin, R.A. The elastic dielectric. J. Ration. Mech. Anal. 1956, 5, 849–915. [CrossRef] 49. Mindlin, R.D.; Eshel, N. On first strain-gradient theories in linear elasticity. Int. J. Solids Struct. 1968, 4, 109–124. [CrossRef] 50. Abdollahi, A.; Millán, D.; Peco, C.; Arroyo, M.; Arias, I. Revisiting pyramid compression to quantify flexoelectricity: A three-dimensional simulation study. Phys. Rev. B 2015, 91, 104103. [CrossRef] 51. Abdollahi, A.; Domingo, N.; Arias, I.; Catalan, G. Converse flexoelectricity yields large piezoresponse force microscopy signals in non-piezoelectric materials. Nat. Commun. 2019, 10, 1266. [CrossRef] 52. He, B.; Javvaji, B.; Zhuang, X. Characterizing Flexoelectricity in Composite Material Using the Element-Free Galerkin Method. Energies 2019, 12, 271. [CrossRef] 53. Liu, G.R. Meshfree Methods: Moving Beyond the Finite Element Method; CRC Press: Boca Raton, FL, USA, 2009. 54. Cordes, L.; Moran, B. Treatment of material discontinuity in the element-free Galerkin method. Comput. Methods Appl. Mech. Eng. 1996, 139, 75–89. [CrossRef] 55. Krongauz, Y.; Belytschko, T. EFG approximation with discontinuous derivatives. Int. J. Numer. Methods Eng. 1998, 41, 1215–1233. [CrossRef] 56. Nanthakumar, S.; Zhuang, X.; Park, H.S.; Rabczuk, T. Topology optimization of flexoelectric structures. J. Mech. Phys. Solids 2017, 105, 217–234. [CrossRef] 57. Deng, F.; Deng, Q.; Shen, S. A three-dimensional mixed finite element for flexoelectricity. J. Appl. Mech. 2018, 85, 031009. [CrossRef] 58. Les Piegl, W.T. The NURBS Book; Springer: Berlin/Heidelberg, Germany, 1997. 59. Cottrell, J.; Hughes, T.; Bazilevs, Y. Isogeometric Analysis: Toward Integration of CAD and FEA; John Wiley & Sons Inc.: Hoboken, NJ, USA, 2009. 60. Camacho, G.T.; Ortiz, M. Computational modelling of impact damage in brittle materials. Int. J. Solids Struct. 1996, 33, 2899–2938. [CrossRef] 61. Schrefler, B.A.; Secchi, S.; Simoni, L. On adaptive refinement techniques in multi-field problems including cohesive fracture. Comput. Methods Appl. Mech. Eng. 2006, 195, 444–461. [CrossRef] 62. de Borst, R. Fluid flow in fractured and fracturing porous media: A unified view. Mech. Res. Commun. 2017, 80, 47–57. [CrossRef] 63. Sukumar, N.; Chopp, D.L.; Moës, N.; Belytschko, T. Modeling holes and inclusions by level sets in the extended finite-element method. Comput. Methods Appl. Mech. Eng. 2001, 190, 6183–6200. [CrossRef] 64. Sukumar, N.; Huang, Z.; Prévost, J.H.; Suo, Z. Partition of unity enrichment for bimaterial interface cracks. Int. J. Numer. Methods Eng. 2004, 59, 1075–1102. [CrossRef] 65. Shu, J.Y.; King, W.E.; Fleck, N.A. Finite elements for materials with strain gradient effects. Int. J. Numer. Methods Eng. 1999, 44, 373–391. [CrossRef] 66. Soh, A.K.; Wanji, C. Finite element formulations of strain gradient theory for microstructures and the C0–1 patch test. Int. J. Numer. Methods Eng. 2004, 61, 433–454. [CrossRef] 67. Gu, S.; He, Q.C. Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces. J. Mech. Phys. Solids 2011, 59, 1413–1426. [CrossRef] 68. Liu, J.T.; He, B.; Gu, S.T.; He, Q.C. Implementation of a physics-based general elastic imperfect interface model in the XFEM and LSM context. Int. J. Numer. Methods Eng. 2018, 115, 1499–1525. [CrossRef] 69. Nguyen, V.P.; Anitescu, C.; Bordas, S.P.; Rabczuk, T. Isogeometric analysis: an overview and computer implementation aspects. Math. Comput. Simul. 2015, 117, 89–116. [CrossRef] Energies 2020, 13, 1326 29 of 29 70. Fischer, P.; Klassen, M.; Mergheim, J.; Steinmann, P.; Müller, R. Isogeometric analysis of 2D gradient elasticity. Comput. Mech. 2011, 47, 325–334. [CrossRef] 71. Shu, L.; Wei, X.; Pang, T.; Yao, X.; Wang, C. Symmetry of flexoelectric coefficients in crystalline medium. J. Appl. Phys. 2011, 110, 104106. [CrossRef] 72. Le Quang, H.; He, Q.C. The number and types of all possible rotational symmetries for flexoelectric tensors. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 2011, 467, 2369–2386. 73. Hong, J.; Vanderbilt, D. First-principles theory and calculation of flexoelectricity. Phys. Rev. B 2013, 88, 174107. [CrossRef] 74. Gao, X.L.; Park, S. Variational formulation of a simplified strain gradient elasticity theory and its application to a pressurized thick-walled cylinder problem. Int. J. Solids Struct. 2007, 44, 7486–7499. [CrossRef] 75. Aravas, N. Plane-strain problems for a class of gradient elasticity models—A stress function approach. J. Elast. 2011, 104, 45–70. [CrossRef] 76. Bursian, E.; OI, Z. Changes in curvature of a ferroelectric film due to polarization. Sov. Phys. Solid State, USSR 1968, 10, 1121–1124. 77. Lumentut, M.; Shu, Y. A unified electromechanical finite element dynamic analysis of multiple segmented smart plate energy harvesters: circuit connection patterns. Acta Mech. 2018, 229, 4575–4604. [CrossRef] 78. Lumentut, M.; Howard, I. Intrinsic electromechanical dynamic equations for piezoelectric power harvesters. Acta Mech. 2017, 228, 631–650. [CrossRef] Sample Availability: Detailed data and codes are available from the authors by request. c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Journal

EnergiesMultidisciplinary Digital Publishing Institute

Published: Mar 12, 2020

There are no references for this article.