Modeling flexoelectricity in soft dielectrics at finite deformation
MModeling flexoelectricity in soft dielectrics at finitedeformation
D. Codony , P. Gupta ,O. Marco , I. Arias , , ∗ Laboratori de C`alcul Num`eric (LaC`aN), Universitat Polit`ecnica de Catalunya (UPC),Campus Nord UPC-C2, E-08034 Barcelona, Spain Centre Internacional de M`etodes Num`erics en Enginyeria (CIMNE),08034 Barcelona, Spain ∗ Corresponding author; E-mail: [email protected].
Abstract
This paper develops the equilibrium equations describing the flexoelectric e ff ect in softdielectrics under large deformations. Previous works have developed related theories using aflexoelectric coupling tensor of mixed material-spatial character. Here, we formulate the modelin terms of a flexoelectric tensor completely defined in the material frame, with the same sym-metries of the small-strain flexocoupling tensor and leading naturally to objective flexoelectricpolarization fields. The energy potential and equilibrium equations are first expressed in termsof deformation and polarization, and then rewritten in terms of deformation and electric po-tential, yielding an unconstrained system of fourth order partial di ff erential equations (PDEs).We further develop a theory of geometrically nonlinear extensible flexoelectric rods underopen and closed circuit conditions, with which we examine analytically cantilever bendingand buckling under mechanical and electrical actuation. Besides being a simple and explicitmodel pertinent to slender structures, this rod theory also allows us to test our general theoryand its numerical implementation using B-splines. This numerical implementation is robust asit handles the electromechanical instabilities in soft flexoelectric materials. Keywords:
Soft dielectrics , Maxwell equations , Flexoelectricity , Electrostriction , Buck-ling , Special Cosserat Rod
Flexoelectricity is a two-way coupling between electric polarization and strain gradient, presentin any dielectric material. The direct flexoelectric e ff ect is understood as the material polarizationdue to inhomogeneous deformation (e.g. bending, twisting), and the converse flexoelectric e ff ect1 a r X i v : . [ phy s i c s . c o m p - ph ] O c t onsists on the generation of stress due to the presence of an inhomogeneous electric field. Theflexoelectric e ff ect is size dependent due to its intrinsic scaling with strain-gradients, and thereforeit is only relevant at the micro- and nanoscale.Flexoelectric e ff ects have been observed and widely studied in hard materials (Tolpygo, 1963,Kogan, 1964, Hong and Vanderbilt, 2011, Resta, 2010, Maranganti et al., 2006), mainly crystallineceramics such as ferroelectric perovskites (Zubko et al., 2007, Ma and Cross, 2001a, 2002, Fuet al., 2006, Ma and Cross, 2001b, 2003, 2005, 2006). However, they are also present in softmaterials, such as liquid crystals (Meyer, 1969, Petrov, 1975, Prost and Marcerou, 1977, Pikin andindenbom, 1978, Lagerwall and Dahl, 1984, Barbero et al., 1986, ˇCepiˇc et al., 2000, Kuczynski andHo ff mann, 2005, Harden et al., 2006, Trabi et al., 2008), cellular membranes (Petrov et al., 1989,Todorov et al., 1991, 1994, Sun, 1997, Petrov, 2002, Gao et al., 2008, Jewell, 2011) and polymers(Breger et al., 1976, Marvan and Havr´anek, 1998, Baskaran et al., 2011, 2012, Deng et al., 2014c,Zhang et al., 2016b, Zhou et al., 2017).The mechanism of flexoelectricity in hard materials can be intuitively understood by the ioniccrystal model under bending, in which a non-zero net dipole moment arises due to a shift betweenthe centers of gravity of the negative and the positive ions (Zubko et al., 2013). In soft materialssuch as liquid crystals or lipid bilayers, flexoelectricity results from the reorientation of irregu-larly shaped polarized molecules under strain gradients, see e.g. Meyer (1969), Petrov (1999), Rey(2006), Ahmadpoor et al. (2013), Liu and Sharma (2013), Mohammadi et al. (2014), Ahmadpoorand Sharma (2015), Morozovska et al. (2018). In these materials, flexoelectricity has been mech-anistically linked to the arrangement of not only dipolar but also quadrupolar constituents (Prostand Marcerou, 1977, Marcerou and Prost, 1980, Derzhanski et al., 1990, de Gennes and Prost,1993), and theories accounting for thermal fluctuations have been proposed (Osipov and Pikin,1995, Liu and Sharma, 2013). The mechanisms leading to flexoelectricity in polymers, however,are not known (Krichen and Sharma, 2016), although they likely involve rearrangements of glassyand crystalline components (Baskaran et al., 2011, 2012). Of note is the conceptual model by Mar-van and Havr´anek (1998), in which flexoelectric polarization results from strain gradient-inducedasymmetry of the free-volume of a fluctuating dipole. This or other mechanisms, however, havenot been demonstrated. We refer to Yudin and Tagantsev (2013), Nguyen et al. (2013), Zubko et al.(2013), Krichen and Sharma (2016), Wang et al. (2019) for excellent and comprehensive reviewsof flexoelectricity in solids.In recent years, several reasons justify an increasing interest in flexoelectricity in polymer ma-terials. On the one hand, a large flexoelectric response is expected. Experiments suggest that theflexoelectric coe ffi cients of polymers are at least the same order of magnitude as those of hardcrystalline materials (Chu and Salem, 2012, Baskaran et al., 2011, 2012), but being much more de-formable, much larger flexoelectric polarization is possible. On the other hand, electromechanicalactuation of polymers by flexoelectricity overcomes the current limitations of traditional actuationbased on electrostriction, which are: (i) one-way coupling, i.e. mechanical deformation does notproduce an electric field, (ii) very large electric fields are required (which may lead to dielectricbreakdown), and (iii) reversal of electric field does not reverse the direction of the deformation(Pelrine et al., 1998, O’Halloran et al., 2008, Krichen and Sharma, 2016, Rosset and Shea, 2016).2urthermore, only a few polymers exhibit significant piezoelectricity (Bauer and Bauer, 2008).Thus, quantifying flexoelectricity at large deformations may enable the design of e ffi cient elec-tromechanical elastomeric devices, such as sensors, actuators and energy harvesters, based on theflexoelectric e ff ect (Jiang et al., 2013, Huang et al., 2018, Wang et al., 2019).The literature about continuum theories of flexoelectricity in bulk solids ranges from the earlyworks by Mashkevich and Tolpygo (1957), Tolpygo (1963), Kogan (1964), Indenbom et al. (1981a,b),Tagantsev (1985, 1986), Sahin and Dost (1988), Tagantsev (1991) to the more recent developmentsby Maranganti et al. (2006), Shen and Hu (2010), Hu and Shen (2010), Hadjesfandiari (2013), Liu(2014), Anqing et al. (2015), to name a few. However, most of these works assume infinitesimaldeformations, and are therefore suitable to model crystalline ceramics only. E ff orts have been re-cently made to extend the theory to polymers or elastomers undergoing large deformations, butthe literature is still scarce (Liu, 2014, Yvonnet and Liu, 2017, Thai et al., 2018, Poya et al., 2019,McBride et al., 2019, Zhuang et al., 2019, Nguyen et al., 2019). Some of these works model flexo-electricity as a linear coupling between strain gradients and the electric displacement (Poya et al.,2019) or the electric field (McBride et al., 2019, Nguyen et al., 2019, Zhuang et al., 2019) insteadof the electric polarization, which however is the most natural choice (Toupin, 1956, Lifshitz andLandau, 1951, Devonshire, 1949, 1951, 1954, Lines and Glass, 1979). Furthermore, works mod-eling flexoelectricity as a coupling between strain gradients and electric polarization consider acoupling tensor of mixed material-spatial character (Liu, 2014, Yvonnet and Liu, 2017, Thai et al.,2018), leading in general to a lack of objectivity in the resulting polarization as argued in Section2.3.The equations of flexoelectricity can only be solved analytically in very simple settings, suchas simplified Euler-Bernoulli (E-B) (Liang et al., 2014, Deng et al., 2014a) and Timoshenko beam(Zhang et al., 2016a) models. Such models have been extended to large deformations but mod-erate rotations `a la von Karmann (Baroudi and Najar, 2019). Otherwise, it is necessary to resortto computational flexoelectricity (Zhuang et al., 2020). The major challenge is to handle the C continuity of the state variables required by the fourth-order PDE system. To address this, severalnumerical alternatives have been proposed, such as mesh-free approximations (Abdollahi et al.,2014, 2015a,b, Abdollahi and Arias, 2015, Zhuang et al., 2019), isogeometric analysis (Ghasemiet al., 2017, Nanthakumar et al., 2017, Thai et al., 2018, Hamdia et al., 2018, Ghasemi et al.,2018, Nguyen et al., 2019), C Argyris triangular element approximation (Yvonnet and Liu, 2017)and the B-spline-based immersed boundary method (Codony et al., 2019). Another family of nu-merical methods are those circumventing the C continuity requirement by introducing additionalvariables, such as mixed formulations (Mao et al., 2016, Deng et al., 2017, 2018), or those based onmicromorphic theories of continua (Poya et al., 2019, McBride et al., 2019). Recently, a few worksreport the application of these methods to large deformation flexoelectricity (Thai et al., 2018, Poyaet al., 2019, McBride et al., 2019, Yvonnet and Liu, 2017, Zhuang et al., 2019, Nguyen et al., 2019)but the continuum formulation at finite deformation is still open, see previous paragraph, and thereis a need for validation of the computational results.To provide a general tool to assess flexoelectricity under large deformations, we propose aformulation with a fully material flexoelectric coupling between strain gradient and electric po-3arization, leading by construction to objective polarization fields. To facilitate the solution of theassociated boundary value problem, we reformulate the balance equations in terms of displace-ments and electric potential as primal unknowns, yielding an unconstrained system of fourth-orderPDE. We solve this system computationally with open uniform B-spline basis in body-fitted Carte-sian meshes. We further derive large deformation models for geometrically nonlinear extensibleflexoelectric rods under open and closed circuit conditions and derive closed-form solutions forcantilever bending and buckling. We report excellent agreement well into the nonlinear regimebetween numerical and analytical solutions in conditions mimicking the assumptions of the ana-lytical models, which serves as validation. We then explore general flexoelectric problems beyondthe simplifying assumptions of the analytical models and analyze the role of the flexoelectric ma-terial parameters in the electromechanical response of the rod.The paper is organized as follows. In Section 2 the free energy density and correspondingbalance equations of a flexoelectric body are reviewed, the mathematical expression of the flex-oelectric coupling is discussed, and the boundary value problem is stated. The numerical imple-mentation used to solve the boundary value problem is presented in Section 3, and the analyticalsolutions for one-dimensional geometrically nonlinear flexoelectric rods are derived in Section 4.In Section 5 the numerical and analytical results of bending and buckling of rods under open / closedcircuit are shown. The paper is concluded in Section 6. Consider a deformable dielectric body described by Ω in the reference (or undeformed) config-uration, and by Ω in the current (or deformed) configuration. The deformation map χ : Ω → Ω maps every material point X ∈ Ω to the spatial point x = χ ( X ) ∈ Ω . Whenever index notationsare used, uppercase and lowercase indexes refer to quantities in the reference and the current con-figurations, respectively. The deformation gradient F , the Jacobian determinant J , and the rightand left Cauchy-Green deformation tensors C , B are defined as F iI ( X ) (cid:66) ∂χ i ( X ) ∂ X I , J (cid:66) det( F ) , C IJ (cid:66) F kI F kJ , B i j (cid:66) F iK F jK . (1)Standard strain measures in the reference and the current configurations are the Green-Lagrangian E and the Almansi-Eulerian e strain tensors given by E IJ (cid:66)
12 ( C IJ − δ IJ ) , e i j (cid:66) (cid:16) δ i j − B − i j (cid:17) = E IJ F − Ii F − J j . (2)Since the flexoelectricity theory involves high-order derivatives, let us define the gradient of thedeformation gradient (cid:101) F , the gradient of the Cauchy-Green deformation tensor (cid:101) C and the Green-4agrangian strain gradient (cid:101) E as (cid:101) F iJK (cid:66) ∂ F iJ ∂ X K = ∂ x i ∂ X J ∂ X K , (cid:101) C IJK (cid:66) ∂ C IJ ∂ X K = IJ (cid:16) (cid:101) F kIK F kJ (cid:17) , (cid:101) E IJK (cid:66) ∂ E IJ ∂ X K = (cid:101) C IJK ; (3)where symm IJ ( A IJ ) : = ( A IJ + A JI ) /
2. Note that the relation (cid:101) E ( (cid:101) F ) in Eq. (3) is inverted as (cid:101) F iJK = (cid:16)(cid:101) E IJK + (cid:101) E KIJ − (cid:101) E KJI (cid:17) F − Ii , (4)analogously to the relation between second derivative of displacement and strain gradients in thelimit of infinitesimal deformation (Schia ffi no et al., 2019).This body in equilibrium necessarily satisfies mechanical balance laws of linear and angularmomentum, and Maxwell equations. In the absence of a magnetic field, they can be expressed inan Eulerian frame as ∇· σ + b = , (5a) σ = σ T , (5b) ∇ × e = , (5c) ∇· d − q =
0; (5d)where σ is the physical stress, e is the the electric field, d is the electric displacement, and b and q are the body force and electric charge per unit volume. Equation (5c) implies the existence of anelectric potential φ such that e = −∇ φ . The linear constitutive law for d for a dielectric material is d ( p , e ) = (cid:15) e + p or, equivalently, d ( p , φ ) = − (cid:15) ∇ φ + p , (6)where p is the electric polarization, which is work-conjugate to e , and (cid:15) is the electric permittivityof vacuum.To formulate the problem in a material frame, the Lagrangian second Piola-Kirchho ff physicalstress tensor S is defined from the work-conjugacy relation σ i j e i j = J S IJ E IJ , (7)where σ i j e i j is a mechanical work density per unit physical volume and S IJ E IJ a mechanical workdensity per unit reference volume, leading to S IJ = JF − Ii F − J j σ i j , (8)where strictly speaking we should write S IJ ◦ χ − = JF − Ii F − J j σ i j to account for the fact that some ofthese fields are over Ω and others are over Ω . To follow an analogous procedure with the electricdisplacement (Lax and Nelson, 1976, Dorfmann and Ogden, 2005, Vu et al., 2007, Dorfmann andOgden, 2014, 2017, Steinmann and Vu, 2017), we first identify the nominal or material electric5eld. The electric potential can be expressed in the material frame as Φ ( X ) = φ ( χ ( X )), and thenominal electric field E defined as the negative of its material gradient. By the chain rule, we thusfind that E I = − ∂ Φ ∂ X I = − ∂φ∂ x i ∂χ i ∂ X I = e i F iI . (9)Then, from the work-conjugacy relation d i e i = J D I E I , (10)we identify the nominal electric displacement as D I = JF − Ii d i . (11)Since electric displacement and polarization are physically equivalent quantities, we analogouslyfind P I = JF − Ii p i . (12)Using Eq. (1), (8), (9), (11), (12), the balance equations in Eq. (5a)-(5d) and the constitutive lawfor dielectrics in Eq. (6) are written in material form as( F iI S IJ ) , J + B i = i , (13a) S IJ = S JI , (13b) E L + Φ , L = , (13c) D K = (cid:15) JC − KL E L + P K , (13d) D K , K − Q = , (13e)with B = J b and Q = Jq . We define the Lagrangian internal energy density per unit reference volume of the flexoelectricsolid as Ψ Int ( E , (cid:101) E , P ) = Ψ Mech ( E , (cid:101) E ) + Ψ Diele ( E , P ) + Ψ Flexo ( P , (cid:101) E ) . (14)We allow Ψ Mech to depend on Lagrangian strain and strain gradient as required for stability (Liu,2014). The isotropic dielectric energy per unit reference volume follows by transforming thespatial expression per unit physical volume ψ Diele ( p ) = (cid:15) − (cid:15) ) p i p i (Liu, 2014) by recallingEq. (12), resulting in Ψ Diele ( E , P ) = J ( (cid:15) − (cid:15) ) P I C IJ P J , (15)6here (cid:15) denotes the electric permittivity of the material. The flexoelectric coupling linking polar-ization and strain gradient is encoded by Ψ Flexo , which for simplicity we assume to be independenton strain.The spatial expression of the electrostatic energy density ψ Elec ( e ) = (cid:15) e i e i (Liu, 2014) can alsobe expressed in the material frame by recalling Eq. (9), resulting in the energy density per unitreference volume Ψ Elec ( E , E ) = J (cid:15) E I C − IJ E J . (16)To formulate a unified potential self-consistently accounting for the material electromechan-ics and for electrostatics, Ψ Int ( E , (cid:101) E , P ) and Ψ Elec ( E , E ) must be expressed in terms of the samevariables. To accomplish this, we resort to a partial Legendre transform and define the followinginternal dual potential ¯ Ψ Int ( E , (cid:101) E , E ) = min P (cid:16) Ψ Int ( E , (cid:101) E , P ) − P · E (cid:17) . (17)The stationarity condition of the minimization results in E ( E , (cid:101) E , P ) = ∂ Ψ Int ∂ P . (18)In principle, this expression can be inverted to find P ( E , E , (cid:101) E ), which plugged into Ψ Int ( E , (cid:101) E , P ) − P · E results in the dual potential ¯ Ψ Int ( E , (cid:101) E , E ).If we postulate the following flexoelectric coupling Ψ Flexo ( P , (cid:101) E ) = − P L f LIJK (cid:101) E IJK , (19)where f LIJK is a purely Lagrangian tensor as further discussed later, this inversion can be madeexplicit yielding E L = J ( (cid:15) − (cid:15) ) C LM P M − f LIJK (cid:101) E IJK ⇒ (20) P M = J ( (cid:15) − (cid:15) ) C − ML (cid:16) E L + f LIJK (cid:101) E IJK (cid:17) = J ( (cid:15) − (cid:15) ) C − ML (cid:16) E L + E Flexo L (cid:17) , (21)where we have defined E Flexo L = f LIJK (cid:101) E IJK for convenience. Replacing this expression for P inEq. (17) and rearranging terms, we find¯ Ψ Int ( E , (cid:101) E , E ) = Ψ Mech ( E , (cid:101) E ) − J (cid:15) − (cid:15) ) E Flexo I C − IJ E Flexo J − J (cid:15) − (cid:15) ) E I C − IJ E J − J ( (cid:15) − (cid:15) ) E I C − IJ E Flexo J . (22)Now, the total electromechanical enthalpy accounting for electrostatics ¯ Ψ Enth = ¯ Ψ Int − Ψ Elec (Liu,2014, Dorfmann and Ogden, 2014, 2017) can be written as¯ Ψ Enth ( E , (cid:101) E , E ) = ¯ Ψ Mech ( E , (cid:101) E ) + ¯ Ψ Diele ( E , E ) + ¯ Ψ Flexo ( E , (cid:101) E , E ) , (23)7ith ¯ Ψ Diele ( E , E ) = − J (cid:15) E M C − ML E L , (24)¯ Ψ Flexo ( E , (cid:101) E , E ) = − JC − ML E M µ LIJK (cid:101) E IJK ; (25)where µ = ( (cid:15) − (cid:15) ) f is the flexoelectricity tensor (Zubko et al., 2013, Wang et al., 2019), describedin Eq. (A.3). The e ff ective mechanical energy density of the system (Wang et al., 2019) is¯ Ψ Mech ( E , (cid:101) E ) =Ψ Mech ( E , (cid:101) E ) − J (cid:15) − (cid:15) ) E Flexo M C − ML E Flexo L , =Ψ Mech ( E , (cid:101) E ) − (cid:101) E IJK (cid:32) µ AIJK JC − AB µ BLMN (cid:15) − (cid:15) (cid:33) (cid:101) E LMN . (26)The standard mechanical contribution accounting for strain gradient elasticity can be written as Ψ Mech ( E , (cid:101) E ) = Ψ Elast ( E ) + (cid:101) E IJK h IJKLMN (cid:101) E LMN , (27)where Ψ Elast can be any classical hyperelastic potential, e.g. Saint-Venant–Kirchho ff , cf. Eq. (A.1),or Neo-Hookean, cf. Eq. (A.2), constitutive models, and h is the sixth-order strain gradient elastic-ity tensor. Upon inspection, it is clear that the second contribution in Eq. (26), i.e. the flexoelectricity-induced mechanical energy, has the same structure as the strain gradient elasticity potential. Forconvenience, we thus define¯ Ψ Mech ( E , (cid:101) E ) =Ψ Elast ( E ) + (cid:101) E IJK ¯ h IJKLMN (cid:101) E LMN , (28)where ¯ h IJKLMN = h IJKLMN − µ AIJK JC − AB µ BLMN (cid:15) − (cid:15) (29)is the e ff ective strain gradient elasticity tensor as described in Eq. (A.4). To preserve the positivedefiniteness of ¯ Ψ Mech , it is clear from Eq. (28) that ¯ h has to be semidefinite positive and thus astability condition can be derived from Eq. (29) depending on both h and µ (Yudin et al., 2014,2015, Morozovska et al., 2016). The boundary of the reference body, ∂ Ω , is split in several disjoint Dirichlet and Neumann sets asfollows: ∂ Ω = ∂ Ω χ ∪ ∂ Ω T = ∂ Ω V ∪ ∂ Ω R = ∂ Ω Φ ∪ ∂ Ω W . (30)8n the Dirichlet boundaries ∂ Ω χ , ∂ Ω V and ∂ Ω Φ , the deformation map χ , normal derivatives ofthe deformation map ∂ N χ , and electric potential Φ are prescribed, respectively. On the Neumannboundaries ∂ Ω T , ∂ Ω R and ∂ Ω W , their respective work conjugate quantities (per unit reference vol-ume) are prescribed, i.e. the surface traction T ( χ , Φ ) = T , the surface double traction R ( χ , Φ ) = R and the surface charge W ( χ , Φ ) = W . As a result of the strain-gradient elasticity potential (Mindlin,1964, Mindlin and Eshel, 1968), additional loads arise in non-smooth regions of ∂ Ω , i.e. edges C in a three-dimensional domain (Mao and Purohit, 2014, Codony et al., 2019). We also split themin Dirichlet in Neumann sets as C = C χ ∪ C J , (31)depending on whether the deformation map χ or edge forces (per unit reference volume) J ( χ , Φ ) = J are prescribed. For simplicity, dead loads are considered.The enthalpy functional governing the physics of a flexoelectric body is written as Π [ χ , Φ ] = (cid:90) Ω (cid:16) ¯ Ψ Enth ( E , (cid:101) E , −∇ Φ ) − B i χ i + Q Φ (cid:17) d Ω − (cid:90) ∂ Ω T T i χ i d Γ − (cid:90) ∂ Ω R R i ∂ N χ i d Γ − (cid:90) C J J i χ i ds + (cid:90) ∂ Ω W W Φ d Γ , (32)where we have used E = −∇ Φ from Eq. (13c). Equilibrium states { χ ∗ , Φ ∗ } are its saddle pointssatisfying { χ ∗ , Φ ∗ } = arg min χ ∈X max Φ ∈P Π [ χ , Φ ] , (33)where X and P are the functional spaces for χ and Φ with su ffi cient regularity fulfilling Dirichletboundary conditions.A necessary condition for equilibrium is the vanishing of the first variation of Π [ χ , Φ ]0 = δ Π [ χ , Φ ; δ χ , δ Φ ] = (cid:90) Ω (cid:32) ∂ ¯ Ψ Enth ∂ E IJ δ E IJ + ∂ ¯ Ψ Enth ∂ (cid:101) E IJK δ (cid:101) E IJK + ∂ ¯ Ψ Enth ∂ E L δ E L − B i δχ i + Q δ Φ (cid:33) d Ω − (cid:90) ∂ Ω T T i δχ i d Γ − (cid:90) ∂ Ω R R i ∂ N δχ i d Γ − (cid:90) C J J i δχ i dS + (cid:90) ∂ Ω Φ W δ Φ d Γ = (cid:90) Ω (cid:16)(cid:98) S IJ δ E IJ + (cid:101) S MJK δ (cid:101) E MJK − D L δ E L − B i δχ i + Q δ Φ (cid:17) d Ω − (cid:90) ∂ Ω T T i δχ i d Γ − (cid:90) ∂ Ω R R i ∂ N δχ i d Γ − (cid:90) C J J i δχ i dS + (cid:90) ∂ Ω Φ W δ Φ d Γ , (34)9or all admissible variations δ χ and δ Φ , and where δ E L (cid:66) − ∂ ( δ Φ ) ∂ X L , δ F iI (cid:66) ∂ ( δχ i ) ∂ X I , δ (cid:101) F iIJ (cid:66) ∂ ( δχ i ) ∂ X I ∂ X J , (35) δ E IJ = δ C IJ (cid:66) symm IJ ( δ F kI F kJ ) , δ (cid:101) E IJK = δ (cid:101) C IJK (cid:66) symm IJ (cid:16) δ F kI (cid:101) F kJK + F kI δ (cid:101) F kJK (cid:17) . (36)We have introduced the local second Piola-Kirchho ff stress (cid:98) S , the second Piola-Kirchho ff doublestress (cid:101) S and the electric displacement D defined as follows: (cid:98) S IJ ( χ , Φ ) = ∂ ¯ Ψ Enth ∂ E IJ = ∂ Ψ Elast ( C ) ∂ C IJ + J C MLIJ E M (cid:32) (cid:15) E L + µ LABK (cid:101) E ABK (cid:33) , (37) (cid:101) S IJK ( χ , Φ ) = ∂ ¯ Ψ Enth ∂ (cid:101) E IJK = ¯ h IJKLMN (cid:101) E LMN − JC − LM E M µ LIJK , (38) D L ( χ , Φ ) = − ∂ ¯ Ψ Enth ∂ E L = JC − KL (cid:16) (cid:15) E K + µ KIJM (cid:101) E IJM (cid:17) , (39)with C ABCD = J ∂ (cid:16) − JC − AB (cid:17) ∂ C CD = (cid:16) C − AC C − BD + C − BC C − AD − C − AB C − CD (cid:17) . (40)Analogously to the infinitesimal strain theory of flexoelectricity (Mao and Purohit, 2014, Codonyet al., 2019), Eq. (34) can be integrated by parts and, by invoking the divergence and surface di-vergence theorems, the strong form in Eq. (13) is recovered along with the following definitionsof the physical second Piola-Kirchho ff stress S , the surface traction T , the double traction R , thesurface charge density W and the edge forces J : S IJ ( χ , Φ ) (cid:66) (cid:98) S IJ ( χ , Φ ) − (cid:101) S IJK , K ( χ , Φ ) = ∂ Ψ Elast ( C ) ∂ C IJ − ¯ h IJKLMN (cid:101) E LMN , K + J C MLIJ E M (cid:15) E L + JC − LM E M , K µ LIJK in Ω , (41a) T i ( χ , Φ ) (cid:66) F iI (cid:104)(cid:16) S IJ ( χ , Φ ) − (cid:101) S IKJ , N P NK (cid:17) N J + (cid:101) S IJK (cid:101) N JK (cid:105) − (cid:101) F iIN P NK (cid:101) S IKJ N J on ∂ Ω , (41b) R i ( χ , Φ ) (cid:66) F iI (cid:101) S IJK N J N K on ∂ Ω , (41c) W ( χ , Φ ) (cid:66) − D L N L on ∂ Ω , (41d) J i ( χ , Φ ) (cid:66) (cid:20) F iI (cid:101) S IJK M J N K (cid:21) on C ; (41e)where N is the outward unit normal vector on ∂ Ω , M is the outward unit co-normal vector on C , P = I − N × N is the projection operator on ∂ Ω , (cid:101) N = ∇ N : P ( N × N ) − ∇ N · P is thesecond-order geometry tensor on ∂ Ω and (cid:126) (cid:127) is the jump operator defined on C as the sum of its10rgument evaluated at each boundary adjacent to C (we refer to Codony et al. (2019) for a detaileddefinition of the quantities involved here).Upon inspection, the second Piola-Kirchho ff stress tensor S in Eq. (41a) is composed by fourterms. The first two terms correspond to the classical and high-order mechanical stresses, respec-tively. The third one corresponds to the total second Piola-Maxwell stress tensor S Maxwell . Thisbecomes evident by expanding it as S Maxwell IJ (cid:66) J C MLIJ E M (cid:15) E L = JF − Ii F − J j (cid:15) (cid:34)(cid:16) E M F − Mi (cid:17) (cid:16) E L F − L j (cid:17) − (cid:16) E M F − Ma (cid:17) (cid:16) E L F − La (cid:17) δ i j (cid:35) , (42)and obtaining its spatial counterpart by using Eq. (8) and (9) as σ Maxwell (cid:66) (cid:15) (cid:32) e ⊗ e − | e | I (cid:33) . (43)The last term in Eq. (41a) corresponds to the total flexoelectricity-induced stress, and is analogousto the term appearing in the linear theory of flexoelectricity, cf. Eqs. (31-33) in Codony et al.(2019).Equations (20) and (39) show that the Lagrangian flexoelectric polarization in the present the-ory is P M = JC − ML ( (cid:15) − (cid:15) ) E L + JC − ML µ LIJK (cid:101) E IJK , and hence its spatial counterpart is derived withEq. (9) and (12) as p m = ( (cid:15) − (cid:15) ) e m + F − Lm µ LIJK (cid:101) E IJK . (44)In the present formulation, µ LIJK is a purely Lagrangian tensor, and hence it is meaningful to view itas a material constant with the same material symmetries and intrinsic symmetry ( µ LIJK = µ LJIK ) asthe small strain flexoelectric tensor (Majdoub et al., 2008, Zubko et al., 2013, Krichen and Sharma,2016). We note, however, that in previous literature a distinct notion of polarization per unitundeformed volume is introduced as p r = J p , i.e. a volume-normalized spatial polarization relatedto our material or nominal polarization by p r i = F iI P I (Liu, 2014, Dorfmann and Ogden, 2014,2017). The polarization p r is not work-conjugate to the Lagrangian electric field E . Furthermore,when it is used to formulate flexoelectric models it can be problematic. Indeed, the flexoelectriccoupling has been defined in terms of p r (Liu, 2014, Deng et al., 2014b,c, Yvonnet and Liu, 2017,Thai et al., 2018) as Ψ Flexo ( (cid:101) F , p r ) = − p r l (cid:70) liJK (cid:101) F iJK , (45)with (cid:70) a mixed spatial-material flexoelectric tensor, which unlike the infinitesimal flexoelectrictensor is intrinsically symmetric with respect to its last two indices ( (cid:70) liJK = (cid:70) liKJ ). By comparingEq. (45) and (19), using Eq. (3) and (4), the relation p r i = F iI P I and the chain rule, we find therelation between f and (cid:70) as (cid:70) liJK = − ∂ Ψ Flexo ∂ p r l ∂ (cid:101) F iJK = symm JK ( f LIJK ) F iI F − Ll , (46a) f LIJK = − ∂ Ψ Flexo ∂ P L ∂ (cid:101) E IJK = (cid:16) (cid:70) liJK F − Ii + (cid:70) l jIK F − J j − (cid:70) lkIJ F − Kk (cid:17) F lL . (46b)11n the limit of infinitesimal deformation, (cid:70) and f correspond to the so-called type-I ( f I ) andtype-II ( f II ) flexocoupling tensors, respectively, and choosing one or the other is just a matter ofconvenience (Schia ffi no et al., 2019). However this equivalence does not hold anymore in a finitedeformation framework, since f is purely Lagrangian whereas (cid:70) is not.Equation (46b) clearly shows that taking (cid:70) as a material constant, as done in Yvonnet and Liu(2017) and Thai et al. (2018), directly implies a very particular dependence of the Lagrangian flex-oelectric tensor f on deformation. Thus, formulating flexoelectricity as in Eq. (45) leads implicitlyto a material flexoelectric tensor whose magnitude and symmetry depend on deformation in a waythat is unphysical. To illustrate this assertion, consider a particular case in which χ corresponds toa rigid body deformation map, and thus F is a rotation matrix R (cid:44) I . Then, Eq. (46b) leads to f LIJK = (cid:16) (cid:70) liJK R Ii + (cid:70) l jIK R J j − (cid:70) lkIJ R Kk (cid:17) R lL (47)showing that, if Eq. (45) is used to model flexoelectricity, then the Lagrangian flexoelectric materialtensor, and hence the enthalpy functional Π [ χ , Φ ], are not invariant with respect to a superimposedrigid body motion and hence not objective. In this Section, we develop a direct numerical approach to solve the boundary value problem inSection 2.3. We restrict ourselves to 2D rod-like geometries, which can be easily discretized byCartesian grids. The state variables { χ , Φ } are approximated by an open uniform B-spline basis(de Boor, 2001, Rogers, 2001, Piegl and Tiller, 2012) of degree p ≥ ξ Figure 1: Univariate open uniform B-spline basis of degree p =
3. Each basis function is a smooth( C p − ) piece-wise polynomial on a compact ( ≤ p +
1) support. Multivariate B-spline bases areconstructed by means of the tensor product of multiple univariate bases.The discretization of Eq. (34) yields a nonlinear system of equations (for the sake of brevity,we keep the same notation to denote discretized quantities). In order to solve it, we considera modified-step Newton-Raphson algorithm. At the k -th iteration, an increment of the solution12 ∆ χ , ∆Φ } ( k ) is found by vanishing the first order Taylor expansion of the residual R in Eq. (34)around the previous solution { χ , Φ } ( k − : R [ χ ( k ) , Φ ( k ) ; δ χ , δ Φ ] ≈ R [ χ ( k − , Φ ( k − ; δ χ , δ Φ ] + ∂ R [ χ ( k − , Φ ( k − ; δ χ , δ Φ ] ∂ χ ∆ χ ( k ) + ∂ R [ χ ( k − , Φ ( k − ; δ χ , δ Φ ] ∂ Φ ∆Φ ( k ) = , (48)leading to an algebraic system of equations for { ∆ χ , ∆Φ } ( k ) of the form (cid:34) H χχ H χ Φ H Φ χ H ΦΦ (cid:35) ( k − · (cid:34) ∆ χ ∆Φ (cid:35) ( k ) = − (cid:34) R χ R Φ (cid:35) ( k − , (49)given { χ , Φ } ( k − at the previous iteration. The explicit form of the variations of the residual R canbe found in Appendix B.Once { ∆ χ , ∆Φ } ( k ) are found, we compute the modified increments of the solution at the k -thiteration, namely { ∆ χ , ∆Φ } ( k ) , by ensuring that the total increment i) leads to an enthalpy decreasealong χ , ii) leads to an enthalpy increase along Φ , and iii) has a predefined maximum norm γ max ∈ R + . The first two conditions are required in accordance to the variational principle in Eq. (33),whereas the latter is just a numerical requirement to avoid too large increments of the solution ateach iteration. To formulate those conditions mathematically, let us recast the variational principlein Eq. (33) as (cid:98) Φ ( χ ) : = argmax Φ ∈P (cid:16) Π [ χ , Φ ] (cid:17) ; (50a) χ ∗ = argmin χ ∈X (cid:16)(cid:98) Π [ χ ] (cid:17) , with (cid:98) Π [ χ ] : = Π [ χ , (cid:98) Φ ( χ )]; (50b) Φ ∗ = (cid:98) Φ ( χ ∗ ) . (50c)Numerically, Eq. (50) is equivalent to solving two linear systems consecutively, constructed fromEq. (48) by writing ∆Φ ( k ) as a function of ∆ χ ( k ) , as follows: (cid:98) H χχ ( k − · ∆ χ ( k ) = − (cid:98) R χ ( k − with (cid:98) H χχ ( k − : = H χχ ( k − − H χ Φ ( k − · H − ΦΦ ( k − · H Φ χ ( k − (cid:98) R χ ( k − : = R χ ( k − − H χ Φ ( k − · H − ΦΦ ( k − · R Φ ( k − ;(51a) H ΦΦ ( k − · ∆Φ ( k ) = − (cid:98) R Φ ( k − with (cid:98) R Φ ( k − : = R Φ ( k − + H Φ χ ( k − · ∆ χ ( k ) . (51b)From Eq. (51) it is clear that the descent and ascent directions are respectively identified by (cid:98) R χ ( k − and (cid:98) R Φ ( k − , i.e. the modified residuals which take into account the coupled nature of the enthalpy13otential. Therefore, the modified increments are computed as follows: α χ ( k ) = − (cid:98) R χ ( k − · ∆ χ ( k ) > , + α Φ ( k ) = − (cid:98) R Φ ( k − · ∆Φ ( k ) < , + β ( k ) = min + , γ max / (cid:115)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∆ χ ( k ) χ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∆Φ ( k ) Φ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ; (52c) ∆ χ ( k ) = α χ ( k ) β ( k ) ∆ χ ( k ) ; (52d) ∆Φ ( k ) = α Φ ( k ) β ( k ) ∆Φ ( k ) ; (52e)with χ and Φ characteristic factors of the problem for displacement and potential. In practice, γ max is treated as an adaptive heuristic parameter, tunable for proper convergence.Finally, the solution at the k -th iteration is updated with { χ , Φ } ( k ) = { χ , Φ } ( k − + { ∆ χ , ∆Φ } ( k ) . (53)The external loads are applied incrementally in a sequence of load steps, and the modified-step Newton-Raphson algorithm presented here is used to obtain converged solutions at every loadstep. Once convergence is reached, the stability of the solution is checked by assuring { χ , Φ } ( k ) is a saddle point in the enthalpy functional Π [ χ , Φ ] in accordance to the variational principle inEq. (33). By means of Eq. (50), stability of { χ , Φ } ( k ) is given by δ χ (cid:98) Π [ χ ( k ) ; ∆ χ ; ∆ χ ] > ∀ ∆ χ ∈ X , (54a) δ φ Π [ χ ( k ) , Φ ( k ) ; ∆Φ ; ∆Φ ] < ∀ ∆Φ ∈ P . (54b)Numerically, Eq. (54) is met by checking the sign of the extremal eigenvalues λ of (cid:98) H χχ ( k ) and H ΦΦ ( k ) as follows: λ min (cid:20) (cid:98) H χχ ( k ) (cid:21) > , λ max (cid:104) H ΦΦ ( k ) (cid:105) < . (55)We recognize convergence to unstable solutions by the violation of Eq. (55). In such case, the so-lution { χ , Φ } ( k ) is slightly perturbed and the iterative algorithm is run again until a stable solution isfound. In practice, we found that λ max (cid:104) H ΦΦ ( k ) (cid:105) remains always negative, and therefore the encoun-tered instabilities are given by λ min (cid:20) (cid:98) H χχ ( k ) (cid:21) becoming negative only (i.e. geometrical instabilities).The eigenvector associated to λ min (cid:20) (cid:98) H χχ ( k ) (cid:21) is an appropriate direction for numerical perturbationson χ ( k ) to reach stable solutions. 14 One-dimensional analytical models for flexoelectric rods un-dergoing large displacements and rotations
In this Section, we derive simplified closed-form solutions for planar bending and buckling offlexoelectric slender uniform rods undergoing large displacements and rotations under open circuitand close circuit conditions. A material point in the reference configuration is denoted by X = X E + X E + S E , where { E , E , E } is a global right-handed orthonormal basis of R , ( X , X )denotes the coordinates of the undeformed cross-section and X = S is the Lagrangian coordinatealong the undeformed arc-length, see Fig. 2. The rod is assumed to be extensible, and thus S is not arc-length of the deformed centerline, but unshearable, following the special Cosserat rodkinematics (Antman, 1995). We further assume that the cross-sections of the rod remain plane andrigid during the deformation. The corresponding deformation map can be defined as χ ( X , X , S ) = r ( S ) + X d + X d , where r ( S ) is the deformed position of the centerline at S E and ( d , d ) arethe director vectors associated with the cross section. For planar bending (Fig. 2), d = − sin θ E + cos θ E , (56) d = E , (57)where θ is the angle of deflection. The deformation gradient F iI = ∂ x i ∂ X I is obtained as ∂ χ ∂ X = d , (58) ∂ χ ∂ X = E , (59) ∂ χ ∂ S = r (cid:48) + X d (cid:48) . (60)where dd S = () (cid:48) . Now, following Antman (1995) and Gupta and Kumar (2017) we have r (cid:48) = ν d , (61) d (cid:48) = − θ (cid:48) d , (62)where ν is the stretch, and d = d × d = cos θ E + sin θ E . (63)Thus, Eq. (60) becomes ∂ χ ∂ S = ( ν − X θ (cid:48) ) d , (64)15he deformation gradient tensor can written as F = cos θ ν − X θ (cid:48) ) sin θ − sin θ ν − X θ (cid:48) ) cos θ , (65)and the Green-Lagrange strain tensor as E = (cid:16) F T F − I (cid:17) = ν − X θ (cid:48) ) − . (66)We rewrite the only non-vanishing component of E as2 E = ( ν − X θ (cid:48) ) − ≈ ν − X ν θ (cid:48) − , (67)where the term X θ (cid:48) has been neglected for thin rods. Expanding ν around ν = E = ζ − X θ (cid:48) , (68)where ζ = ν − S , the dominant strain gradientcomponent is (cid:101) E = − θ (cid:48) . (69) E E d d r ( S ) (cid:1) L H + ( X ,X ,S ) Figure 2: A typical schematic of deformed planar rod of length L and height H from its referencestraight configuration. For upward bending, θ > ff ect of E , theequilibrium condition Eq. (34), reduces to (cid:90) L (cid:34)(cid:90) A (cid:98) S δ E d A + (cid:90) A (cid:101) S δ (cid:101) E d A − (cid:90) A D δ E d A (cid:35) d X − δ ˆ T + δ ˆ W = . (70)16 L H (a) A cantilever rod subjected to an endpoint load N at the right end tip, and electrically grounded at themid-point in the right end cross-section. NHL (b) A clamped-clamped rod subjected to a com-pressive load N at the right end, and electricallygrounded at the mid-point in the left end cross-section. However, the right end is allowed to dis-place horizontally. L H V (c) A cantilever actuator sandwiched between twoelectrodes (blue in color) under voltage V . V L H (d) A clamped-clamped actuator sandwiched be-tween two electrodes (blue in color) under voltage V . Figure 3: A schematic of flexoelectric rod under external mechanical load or external voltage.where L is the undeformed length of the rod, A is area of the cross-section, and δ ˆ T and δ ˆ W are thevariations of the external work done by mechanical tractions and surface charges. Since bendingof slender rods can involve large displacements but typically small Lagrangian strains, all isotropicconstitutive models are very close. For convenience, we consider the isotropic Kirchho ff -Saint-Venant model, requiring two elastic constants, here Young’s modulus Y and Poisson’s ratio ν , seeEq. (A.1). The flexoelectric tensor µ is assumed to have cubic symmetry with three independentconstants µ L , µ T and µ S , namely the longitudinal, transversal and shear coe ffi cients (Eq. (A.3)). Weassume that all material properties are homogeneous in the cross-section.Using Eqs. (68) and (69) in Eq. (70), the corresponding local stress, higher order stress andelectric displacement relations in Eq. (37)-(39) reduce to (cid:98) S = ¯ Y ( ζ − X θ (cid:48) ) − (1 − ζ + X θ (cid:48) ) (cid:32) − E µ T θ (cid:48) + (cid:15) E (cid:33) , (71) (cid:101) S = − (1 + ζ − X θ (cid:48) ) µ T E , (72) D = (cid:0) + ζ − X θ (cid:48) (cid:1) (cid:0) (cid:15) E − µ T θ (cid:48) (cid:1) , (73)with ¯ Y = Y (1 − ν ) / (1 + ν )(1 − ν ). Note that in this reduced order theory, only transverse flexo-electricity is relevant. 17 .1 Flexoelectric rod in open circuit under mechanical load We consider now a flexoelectric rod in open circuit conditions, i.e. one of the rod end’s is groundedand all other boundaries are free of surface charges, i.e. they satisfy that D · n =
0, see Figs. 3a, 3b.Thus, at the top and bottom surfaces, the vertical electric displacement vanishes, D =
0, and forthin rods it can be assumed to vanish within the cross-section as well (Majdoub et al., 2008, 2009,Liang et al., 2014). In this case, the vertical electric field can be computed from Eq. (73) as E = µ T (cid:15) θ (cid:48) , (74)and then, Eqs. (71) and (72) reduce to (cid:98) S = ¯ Y ( ζ − X θ (cid:48) ) + (1 − ζ + X θ (cid:48) ) µ θ (cid:48) (cid:15) , (75) (cid:101) S = − (1 − ζ + X θ (cid:48) ) µ (cid:15) θ (cid:48) . (76)By substituting Eqs. (75) and (76) into Eq. (70), and using ˆ W = (cid:90) L (cid:40) ¯ Y A (cid:34) ζ +
12 (1 − ζ ) (cid:96) µ θ (cid:48) (cid:35) δζ + ¯ Y (cid:34) I (cid:32) − (cid:96) µ θ (cid:48) (cid:33) + (1 − ζ ) (cid:96) µ A (cid:35) θ (cid:48) δθ (cid:48) (cid:41) d X − δ ˆ T = , (77)where I = (cid:82) A X d A is the moment of inertia of the cross-section and (cid:96) µ = µ T / √ ¯ Y (cid:15) is a lengthscalearising from transversal flexoelectricity. Since the stretch in thin rods is expected to be small, evenif deformations are not, we can approximate 1 − ζ ≈ (cid:90) L (cid:40) ¯ Y A (cid:34) ζ + (cid:96) µ θ (cid:48) (cid:35) δζ + ¯ Y (cid:34) I (cid:32) − (cid:96) µ θ (cid:48) (cid:33) + (cid:96) µ A (cid:35) θ (cid:48) δθ (cid:48) (cid:41) d X − δ ˆ T = , (78)where we identify the axial force and the bending moment as N = ¯ Y A (cid:34) ζ + (cid:96) µ θ (cid:48) (cid:35) , (79) M = ¯ Y (cid:34) I − (cid:96) µ θ (cid:48) I + (cid:96) µ A (cid:35) θ (cid:48) . (80)Interestingly, Eq. (78) points out the two main size-dependent e ff ects of flexoelectricity. On onehand, flexoelectricity induces a positive size-dependent axial strain in the rod which dependsquadratically on the flexural strain θ (cid:48) . On the other hand, flexoelectricity modifies the e ff ectivebending sti ff ness by two size-dependent contributions of opposite sign. The first is a reduction inthe rod’s sti ff ness which depends quadratically on the flexural strain while the second makes therod sti ff er independent of deformation. 18n order to evaluate the relative importance of the di ff erent contributions, let us consider arectangular cross-section of unit width and thickness H , i.e. I = H /
12 and A = H . The secondand third contribution to the e ff ective bending sti ff ness are comparable in magnitude for a radiusof curvature R = ν /θ (cid:48) ≈ H /
5, which is unphysically small. For reasonable radii of curvature, weexpect 12 (cid:96) µ θ (cid:48) I << (cid:96) µ A , and thus Eq. (78) reduces to (cid:90) L (cid:34) ¯ Y A (cid:34) ζ + (cid:96) µ θ (cid:48) (cid:35) δζ + ¯ Y I e ff θ (cid:48) δθ (cid:48) (cid:35) d X − δ ˆ T = , (81)with I e ff = I + (cid:96) µ A . (82)Furthermore, the values of (cid:96) µ for typical flexoelectric polymers are in the order of 1 −
10 nm (Chuand Salem, 2012, Zhang et al., 2015, Zhou et al., 2017). The minimum radius of curvature for arectangular cross-section is R = H /
2, which implies that the maximum flexoelectrically-inducedaxial strain is approximately 2 (cid:96) µ / H , in the order of 10 − for a H = ζ . This is later verified inthe numerical examples in Section 5.1, with ζ in the order of 10 − for a H = ζ + (cid:96) µ θ (cid:48) /
2, we consider now a flexoelectric cantileverrod subjected to a point load N = N E + N E on one of its ends, Figs. 3a, 3b. The work done bythe external force isˆ T = N · r ( L ) = N · (cid:90) L (1 + ζ ) d ( S ) d S = (cid:90) L (1 + ζ ) ( N sin θ + N cos θ ) d S , (83)where we have used Eq. (61). Substituting the first variation of Eq. (83) in Eq. (81), and assumingthat the stretch is small, i.e. 1 + ζ ≈ (cid:90) L (cid:34) ¯ Y A (cid:34) ζ + (cid:96) µ θ (cid:48) (cid:35) δζ + ¯ Y I e ff θ (cid:48) δθ (cid:48) (cid:35) d S = (cid:90) L (cid:2) ( N sin θ + N cos θ ) δζ + ( N cos θ − N sin θ ) δθ (cid:3) d S . (84)Upon integration by parts, Eq. (84) becomes (cid:90) L (cid:34) ¯ Y A δ (cid:34) ζ + (cid:96) µ θ (cid:48) (cid:35) − ¯ Y I e ff θ (cid:48)(cid:48) δθ (cid:35) d S + ¯ Y I e ff θ (cid:48) δθ (cid:12)(cid:12)(cid:12) L = (cid:90) L (cid:2) ( N sin θ + N cos θ ) δζ + ( N cos θ − N sin θ ) δθ (cid:3) d S , (85)from where the Euler-Lagrange equations can be derived for all admissible δζ and δθ as¯ Y A ζ + ¯ Y A (cid:96) µ θ (cid:48) − N sin θ − N cos θ = , (86a)19 Y I e ff θ (cid:48)(cid:48) + N cos θ − N sin θ = , (86b)where we have assumed that the external force N is known. Equations (86) form a system of twocoupled equations for the two unknowns ζ and θ , where θ can be obtained from Eq. (86b) and usedin Eq. (86a) to compute ζ . Note that Eq. (86b) corresponds to bending moment balance of a purelymechanical non-linear Kirchho ff rod with modified (larger) bending rigidity I e ff (Antman, 1995).This e ff ective sti ff ness coincides with that identified by Majdoub et al. (2008, 2009), Liang et al.(2014) for linear flexoelectric rods. Equation (86b) can be rewritten in standard form as θ (cid:48)(cid:48) + ¯ β N · d = , (87)with ¯ β − = ¯ Y I e ff .We derive next the solution for bending of a cantilever flexoelectric rod under a vertical pointload, and buckling of a doubly clamped rod under axial compression, see Fig. 3a. We consider a cantilever flexoelectric rod subjected to a vertical force N = − N E , see Fig. 3a. Inthis case, Eq. (87) reduces to θ (cid:48)(cid:48) − β cos θ = , (88)with β = ¯ β N and boundary conditions θ (0) = , (89a) θ (cid:48) ( L ) = . (89b)The solution to this problem was obtained by Bisshopp and Drucker (1945). As derived in detailin Appendix C, the vertical displacement of the tip is r ( L ) = L + β (cid:104) ˜ E ( p , ψ ) − ˜ E ( p ) (cid:105) , (90)where ˜ E ( p ) and ˜ E ( p , ψ ) are the complete and incomplete elliptical integrals of the second kind,respectively, see Eq. (C.12), with p = √ (1 − sin θ max ) /
2, and 1 / sin ψ = √ − sin θ max , with θ max = θ ( L ). For a given load N , θ max is obtained by the shooting method, using Eq. (C.9). UsingEqs. (86a) and (C.2), the axial strain can be computed as ζ ( S ) = − N ¯ Y A sin θ − β (cid:96) µ (sin θ − sin θ max ) , (91)which attains its maximum value (in magnitude) at the free end ζ ( L ) = − N ¯ Y A sin θ max . (92)20ote that for N >
0, the rod bends downwards ( θ <
0) while for N <
0, the rod bends upwards( θ >
0) and thus in all cases ζ > E (0) = µ T (cid:15) θ (cid:48) (0) = − µ T (cid:15) β (cid:112) | θ max | = − (cid:96) µ (cid:114) N (cid:15) I e ff sin | θ max | . (93)In the limit case of small deflections, or small N , we recover the well-known flexoelectrictheory relying on linear Euler-Bernoulli beams (Majdoub et al., 2008, 2009), yielding the verticaldisplacement at the free end and the curvature at the fixed end as r ( L ) = − N L Y I e ff , (94) θ (cid:48) (0) ≈ − N L ¯ Y I e ff , (95)and thus the electric field at the fixed end for the linear Euler-Bernoulli beam is E (0) = µ T (cid:15) θ (cid:48) (0) ≈ − µ T N L (cid:15) ¯ Y I e ff . (96)Interestingly, for a given rod of length L , with cross-section area A and moment of inertia I , thevertical electric field at the fixed end in Eq. (96) attains a maximum for µ ∗ T = (cid:114) (cid:15) ¯ Y IA . (97)From a physical point of view, this maximum is a result of two competing e ff ects of flexoelectricity.On one hand, flexoelectricity increases the bending rigidity of the rod by increasing the e ff ectivemoment of inertia in (cid:96) µ A ∝ µ , see Eq. (82), and thus reduces the rod deflection, the curvatureand the resulting vertical electric field. On the other hand, for a given curvature, the verticalelectric field is proportional to µ T . The maximum electric field at the fixed end for this optimumflexoelectric coe ffi cient µ ∗ T becomes E max1 (0) = − N L √ (cid:15) ¯ Y IA . (98)Similarly, for a given material with properties Y , (cid:15), and µ T , one can find an optimal design thatmaximizes the flexoelectric response. For instance, considering a rod with square / rectangularcross-section, the optimal thickness is H ∗ = µ T (cid:114) Y (cid:15) = √ (cid:96) µ , (99)and the corresponding maximum vertical electric field at the fixed end for a rod with a unit widthis E max1 (0) = − √ (cid:15) ¯ Y N L √ µ . (100)21 .1.2 Buckling of a flexoelectric rod under axial compression We consider next an open-circuit flexoelectric rod clamped at the left end and with vertical dis-placement and rotation prevented at the right end, subjected to a compressive force N = − N E ,see Fig. 3b. We examine the buckling critical load N cr and the post-buckling behavior. In this case,the vertical reaction at the right end N is unknown. Thus, the Euler-Lagrange Eqs. (86) have to besupplemented with the constraint of vanishing vertical displacements at the right end given by0 = E · r ( L ) = E · (cid:90) L (1 + ζ ) d ( s ) d s = (cid:90) L (1 + ζ ) sin θ d s . (101)The bending moment balance Eq. (87) is written as θ (cid:48)(cid:48) + β sin θ + N ¯ β cos θ = , (102)with β = N / ¯ Y I e ff , subject to the boundary conditions θ (0) = , (103a) θ ( L ) = . (103b)Since the expected lowest buckling mode is symmetric, the constraint in Eq. (101) is fulfilled bysymmetry, and thus the reaction N =
0, and Eq. (102) reduces to θ (cid:48)(cid:48) + β sin θ = . (104)The expected lowest mode exhibits inflection points at S = L / S = L /
4, which requirespecial attention (Lin and Chiao, 1998). Instead, we invoke symmetry considerations and avoidthe inflection points by solving Eq. (104) over a quarter of the rod and replace Eq. (103b) with θ (cid:48) (cid:18) L (cid:19) = . (105)After solving the BVP, see Appendix D for a detailed derivation, the vertical displacement and thevertical electric field at the center of the rod for upward buckling are, r (cid:18) L (cid:19) = β (cid:112) − cos θ max ) , (106a) E (cid:18) L (cid:19) = − µ T (cid:15) (cid:114) N (1 − cos θ max )¯ Y I e ff . (106b)where θ max = θ ( L /
4) is computed for a given load N by the shooting method using Eq. (D.6). Sincethe right end of the rod is allowed to move horizontally, its length is assumed to remain unchangedand the stretch is ν ≈ N = β ¯ β = F (cid:32) sin θ max (cid:33) ¯ Y I e ff L , (107)22here F is the complete elliptical integral of first kind, cf. Eq. (C.10), and we have used Eq. (D.6).By letting θ max → N cr = π ¯ Y I e ff L , (108)which coincides with the buckling load for a linear flexoelectrically-sti ff ened Euler-Bernoulli beamwith a modified bending sti ff ness (Timoshenko and Gere, 2009), see Eq. (82). We consider a flexoelectric rod in closed circuit, i.e. electrodes are attached to the top and bottomsurfaces. Under actuation operation mode, i.e. the bottom electrode is grounded ( φ = φ = V is applied to the top electrode, two setups are studied, bending of a cantilever,Fig. 3c, and buckling of a doubly clamped rod, Fig. 3d.In these setups, neglecting the localized boundary e ff ects at the ends of the rod, the non-vanishing electric field component is E = − VH , (109)and then, Eqs. (71) – (73) reduce to (cid:98) S = ¯ Y ( ζ − X θ (cid:48) ) − (1 − ζ + X θ (cid:48) ) (cid:32) VH µ T θ (cid:48) + (cid:15) V H (cid:33) , (110) (cid:101) S = (1 + ζ − X θ (cid:48) ) µ T VH , (111) D = (cid:0) + ζ − X θ (cid:48) (cid:1) (cid:18) − (cid:15) VH − µ T θ (cid:48) (cid:19) . (112)Hence, the balance law in Eq. (70) becomes (cid:90) L (cid:40) (cid:34)(cid:32) ¯ Y + µ T VH θ (cid:48) + (cid:15) V H (cid:33) ζ − µ T VH θ (cid:48) − (cid:15) V H (cid:35) A δζ + (cid:34)(cid:32) ¯ Y + µ T VH θ (cid:48) + (cid:15) V H (cid:33) I θ (cid:48) − (1 + ζ ) µ T VH A (cid:35) δθ (cid:48) (cid:41) d s − δ ˆ T = , (113)where δ ˆ T is given by Eq. (83) for an external force N = N E + N E applied at the right end, andwe have used δ E = δ ˆ W =
0. Assuming again that the strain is small, integration by parts yields, (cid:90) L (cid:40)(cid:34)(cid:18) ˜ Y + µ T VH θ (cid:48) (cid:19) ζ − µ T VH θ (cid:48) − (cid:15) V H (cid:35) A δζ − (cid:20)(cid:18) ˜ Y + µ T VH θ (cid:48) (cid:19) I θ (cid:48) − (1 + ζ ) µ T VH A (cid:21) (cid:48) δθ (cid:41) d s − δ ˆ T = − (cid:20)(cid:18) ˜ Y + µ T VH θ (cid:48) (cid:19) I θ (cid:48) − (1 + ζ ) µ T VH A (cid:21) δθ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) θ ( L ) θ (0) , (114)where we have defined an e ff ective Young’s modulus modified by electrostriction as˜ Y = ¯ Y + (cid:15) V H . (115)23 .2.1 Bending of a flexoelectric cantilever under applied voltage We consider next a flexoelectric cantilever rod sandwiched between two electrodes as depicted inFig. 3c. In this case, there are no applied mechanical loads and there is no kinematical constraintat the right end, and thus N = . The Euler-Lagrange equations are identified as (cid:18) ˜ Y + µ T VH θ (cid:48) (cid:19) ζ − µ T VH θ (cid:48) − (cid:15) V H = , (116a) (cid:20)(cid:18) ˜ Y + µ T VH θ (cid:48) (cid:19) I θ (cid:48) − (1 + ζ ) µ T VH A (cid:21) (cid:48) = , (116b)Equation (116a) yields ζ = µ T VH θ (cid:48) + (cid:15) V H ˜ Y + µ T VH θ (cid:48) ≈ (cid:15) V H ˜ Y + µ T V ˜ Y H (cid:32) − (cid:15) V Y H (cid:33) θ (cid:48) , (117)where we have expanded ζ in a Taylor series around θ (cid:48) = (cid:34) ˜ Y (cid:32) I − µ AV ˜ Y H (cid:32) − (cid:15) V Y H (cid:33)(cid:33) θ (cid:48) − µ T VH A (cid:32) + (cid:15) V Y H (cid:33)(cid:35) (cid:48) = . (118)By defining, an e ff ective moment of inertia modified by flexoelectricity and electrostriction as˜ I = I − µ AV ˜ Y H (cid:32) − (cid:15) V Y H (cid:33) , (119)and an e ff ective cross-section area modified by electrostriction as˜ A = A (cid:32) + (cid:15) V Y H (cid:33) , (120)Eq. (118) reduces to (cid:20) ˜ Y ˜ I θ (cid:48) − µ T VH ˜ A (cid:21) (cid:48) = . (121)Equation (121) implies that the flexural strain θ (cid:48) is uniform along s . By Eq. (114), the correspond-ing boundary conditions are θ (0) = , (122)˜ Y ˜ I θ (cid:48) ( L ) − µ T VH ˜ A = . (123)The resulting uniform flexural strain in this case is θ (cid:48) = θ (cid:48) ( L ) = µ T VH ˜ A ˜ Y ˜ I , (124)24hich agrees with the expression given in Bursian and Trunov (1974) for a linear flexoelectric rodby replacing the e ff ective quantities ( ˜ Y , ¯ I , ¯ A ) by the nominal ones ( ¯ Y , I , A ). From Eq. (117) theaxial strain is obtained as ζ = µ ˜ A ˜ Y ˜ I V H (cid:32) − (cid:15) V Y H (cid:33) + (cid:15) V H ˜ Y . (125)We now examine Eqs. (124) and (125) by Taylor expansion of these expressions around V / H = θ (cid:48) = (cid:18) VH (cid:19) A µ T I ¯ Y + (cid:18) VH (cid:19) A µ T I ¯ Y + O (cid:32)(cid:18) VH (cid:19) (cid:33) , (126a) ζ = (cid:18) VH (cid:19) (cid:15) ¯ Y + A (cid:96) µ I − (cid:18) VH (cid:19) (cid:18) (cid:15) ¯ Y (cid:19) + A (cid:96) µ I − A (cid:96) µ I + O (cid:32)(cid:18) VH (cid:19) (cid:33) ≈ (cid:18) VH (cid:19) (cid:15) ¯ Y − (cid:18) VH (cid:19) (cid:18) (cid:15) ¯ Y (cid:19) + O (cid:32)(cid:18) VH (cid:19) (cid:33) (126b)According to Eq. (126a), under the application of a voltage V , the flexoelectric cantilever bendsupwards for V > V < V , due to both flexoelectricity and electrostriction from Eq. (126b).However, the contribution of the flexoelectric e ff ect on the axial strain is negligible as for typicalflexoelectric elastomers A (cid:96) µ / I ≈ − (cid:28) H = R = θ (cid:48) (1 + ζ ) ≈ µ T AI ˜ Y VH . (127)Integrating θ from Eq. (127) and accounting for the clamping condition Eq. (122), we have θ ( S ) = (1 + ζ ) µ T AI ˜ Y VH S . (128)Finally, the vertical deflection at the free end can be evaluated as r ( L ) = (cid:90) L (1 + ζ ) sin θ d S = (1 + ζ ) (cid:90) L sin (cid:18) (1 + ζ ) µ T AI ˜ Y VH S (cid:19) d S = ˜ Y I µ T A HV (cid:20) − cos (cid:18) (1 + ζ ) µ T ALI ˜ Y VH (cid:19)(cid:21) ≈ ˜ Y I µ T A HV (cid:20) − cos (cid:18) µ T ALI ˜ Y VH (cid:19)(cid:21) . (129) We consider now a doubly clamped flexoelectric rod in closed circuit conditions subjected to anexternal electrical bias V , Fig. 3d. Since as we have seen above, an applied bias leads to an elon-gation of the rod, if axially constrained this should lead to buckling, and hence here we study the25ritical buckling load V cr and the post-buckling behavior. Similarly to Section 4.1.2, the kinematicconstraints of vanishing vertical and horizontal displacements at the right end, give rise to a re-action force at the right end N = N E + N E , where now N and N are unknown quantities.Since the expected lowest buckling mode is symmetric, the vertical displacement at the right endvanishes by symmetry and thus N =
0. Hence, Eq. (83) reduces toˆ T = (cid:90) L (1 + ζ ) N cos θ d S , (130)and its variation is δ ˆ T = (cid:90) L (cid:2) N cos θ δζ − (1 + ζ ) N sin θ δθ (cid:3) d S . (131)Replacing Eq. (131) in Eq. (114), the Euler-Lagrange equations are derived as (cid:18) ˜ Y + µ T VH θ (cid:48) (cid:19) A ζ − µ T A VH θ (cid:48) − (cid:15) V H A − N cos θ = , (132a) (cid:20)(cid:18) ˜ Y + µ T VH θ (cid:48) (cid:19) I θ (cid:48) − (1 + ζ ) µ T VH A (cid:21) (cid:48) − (1 + ζ ) N sin θ = , (132b)and the constraint of vanishing horizontal displacement at the right end is E · r ( L ) − L = E · (cid:90) L (1 + ζ ) d ( S ) d S − L = (cid:90) L (1 + ζ ) cos θ d S − L = . (133)The unknown reaction force magnitude N is calculated by evaluating Eq. (132a) at the left end,with θ (0) =
0, as N = (cid:18) ˜ Y + µ T VH θ (cid:48) | (cid:19) A ζ | − µ T A VH θ (cid:48) | − (cid:15) V H A . (134)Furthermore, by assuming that the axial strain and all material parameters are uniform along S andneglecting the nonlinear term 2 µ T VH θ (cid:48) I θ (cid:48)(cid:48) , Eq. (132b) reduces to˜ Y I θ (cid:48)(cid:48) − (1 + ζ ) N sin θ = . (135)Finally, substituting N from Eq. (134), Eq. (135) simplifies to θ (cid:48)(cid:48) + (1 + ζ ) ˜ β sin θ = , (136)with ˜ β = A ˜ Y I (cid:32) − ˜ Y ζ + (cid:15) V H + (1 − ζ ) µ T VH θ (cid:48) | (cid:33) , (137)26nd subject to the boundary conditions θ (0) = , (138a) θ ( L ) = . (138b)Similarly to the problem in Section 4.1.2, the expected lowest mode exhibits inflection points at S = L / S = L /
4. To avoid having to deal with them, we consider only a quarter of the rodand replace Eq. (138b) with θ (cid:48) (cid:18) L (cid:19) = , (139)After solution of the above BVP, see Appendix E for a detailed derivation, the vertical displacementat the center of the rod is obtained as r (cid:18) L (cid:19) = − (cid:112) + ζ ˜ β sin θ max , (140)where θ max = θ ( L / β and ζ are computed from Eqs. (E.3) and (E.5) in terms of θ max . Thecurvature at the left end is θ (cid:48) (0) = ˜ β (cid:112) + ζ )(1 − cos θ max ) . (141)Finally, using Eq. (137) and (141), the postbuckling voltage can be obtained as V = H (cid:112) Y /(cid:15) (1 − ζ ) − ˜ β I / A (cid:32) (cid:113)(cid:16) ˜ β ( ζ + ζ − (cid:96) µ (1 − cos θ max ) − (cid:16) ( ζ − + ˜ β I / A (cid:17) (cid:16) ζ + ˜ β I / A (cid:17)(cid:17) + ( ζ −
1) ˜ β(cid:96) µ (cid:112) ( ζ + − cos θ max ) (cid:33) (142)The critical buckling voltage is determined from Eq. (142) in the limit θ max → V cr = π HL (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) Y (cid:15) (cid:32) AI − π L (cid:33) . (143)The critical electric field for a rectangular / square cross section becomes: E cr = π L (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) Y (cid:15) (cid:32) H − π L (cid:33) = (cid:18) HL (cid:19) (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) Y (cid:15) (cid:32) π − (cid:18) HL (cid:19) (cid:33) , (144)and for slender rods, the Taylor approximation around H / L → E cr ≈ (cid:18) HL (cid:19) π (cid:114) Y (cid:15) . (145)27 Numerical examples for general nonlinear flexoelectric rodproblems
In this Section, we present numerical results of our general nonlinear model of flexoelectricity forbending and buckling of flexoelectric rods, both in open-circuit and in closed-circuit conditions.We compare these results with the solutions of the 1D nonlinear analytical model for rods devel-oped in Section 4 and its linearized Euler-Bernoulli (E-B) counterpart, by considering materialparameters to match the assumptions of these models. This comparison allows us to validate ourcomputational approach. We then explore more general flexoelectric problems and establish thelimits of the simplified 1D flexoelectric rod models.To model standard elasticity, we consider isotropic hyperelastic potentials, either Saint-Venant–Kirchho ff (Eq. (A.1)) or Neo-Hookean (Eq. (A.2)) models, requiring two elastic constants, hereYoung’s modulus Y and Poisson’s ratio ν . Strain-gradient elasticity is modeled by the analo-gous isotropic hyperelastic Saint-Venant–Kirchho ff law (Eq. (A.4)), which additionally dependson the characteristic length scale (cid:96) . The flexoelectric tensor µ is assumed to have cubic symmetrywith three independent constants µ L , µ T and µ S , namely the longitudinal, transversal and shearcoe ffi cients (Eq. (A.3)). Isotropic flexoelectricity tensor is just a particular case with only twoindependent parameters, with 2 µ S = µ L − µ T .The dielectric strength (i.e. maximum electric field magnitude that a dielectric can sustain be-fore electric breakdown occurs) is typically around 1 − / µ m (Liu, 2014). Here, for simplicity,electrical breakdown is neglected, i.e. we assume an infinite dielectric strength in all the examples.In all simulations, we consider a cubic ( p =
3) spline mesh with square cells of size h = H / H the thickness of the rod. We consider here a flexoelectric cantilever rod under bending by a vertical point load in an opencircuit configuration with the mechanically free end electrically grounded, cf. Fig. 3a. Young’smodulus is chosen as Y = . (cid:15) = . / V m, whichcorrespond to polyvinylidene fluoride (PVDF) (Chu and Salem, 2012, Zhang et al., 2016b, Zhouet al., 2017). We first validate the full computational model in Section 2 and 3 against the 1D nonlinear modelfor flexoelectric rods presented in Section 4.1.1, and its linearized Euler-Bernoulli counterpart.For this, we choose a Saint-Venant–Kirchho ff mechanical constitutive law with ν = µ L = µ S = (cid:96) =
0. We consider a thickness H = L / H ≥
20. Fig. 4 collects allthe validation results. Typical computational solutions are shown in Fig. 4a, where the electric28otential φ is plotted on the deformed configuration. These simulations highlight the very largedeformations attained. In this figure, we show numerical calculations for a given force at the tip,and for several values of µ T . As in the linear case (Majdoub et al., 2008, 2009) and as expected bythe expression I e ff in the reduced theory, cf. Eq. (82), flexoelectricity leads to an e ff ective sti ff eningof the system even though the elastic constants are kept fixed. As anticipated in Section 4.1.1 for thelinearized Euler-Bernoulli beam, cf. Eq. (97), we find that also for the non-linear rod the maximumelectric field generated at the clamping cross-section exhibits a maximum for an intermediate valueof the flexoelectric constant. The existence of an optimal value of µ T , for which the flexoelectricresponse is maximized results from the competition of the two conflicting e ff ects of µ T : (1) thesti ff ening and (2) the flexoelectric coupling. For small values of µ T the structure is very compliantand larger strain gradients are attained but the generated field is small due to the small coupling,whereas for very large values of µ T the flexoelectric coupling is large but the sti ff er beam attainssmaller deformations and thus smaller strain gradients.To further analyze these e ff ects, we present in Fig. 4b the dependence of the cantilever rodvertical displacement at the tip on the endpoint load, and the vertical electric field on the clampededge, for di ff erent values of transversal flexoelectric coe ffi cient µ T . The results for the tip dis-placement show i) the sti ff ening as µ T increases, ii) the nonlinearity in the response of the system(particularly for the most deformable systems), iii) an excellent quantitative agreement with thenonlinear flexoelectric rod model given by the analytical expression in Eq. (90), and iv) an agree-ment with the linearized E-B model for small deformations, i.e. smaller loads or sti ff er cantilevers(large values of µ T ). Similarly, we find an excellent agreement between the numerical simulationsand the nonlinear rod model in the vertical electric field on the clamped end. Its behavior is nonlin-ear for large loads since the electric field is directly proportional to the curvature, cf. Eq. (74). Thenon-monotonicity in the maximum electric field as a function of µ T discussed above is apparentfrom this plot. To further examine this point, we represent in Fig. 4c a contour plot showing thedependence of the vertical electric field at the clamped cross-section on µ T and on the load. Wefind that the load for maximum electrical output depends on the value of the flexoelectric couplingin the nonlinear model, whereas it is independent of it according to the linearized E-B model, seeEq. (98).Finally, we examine the e ff ect of the slenderness on the load vs. deflection and the load vs. elec-tric field curves for a given µ T , see Fig. 4d. As the slenderness ratio increases, the rod becomesmore flexible and therefore nonlinearity is stronger and manifests for smaller loads, with a largeroverestimation of the vertical displacement by the linear E-B model. In contrast, the nonlinear 1Drod model closely follows our simulations even deep into the nonlinear regime. We investigate now more general flexoelectric conditions beyond the restrictive assumptions of thereduced model in Section 4.1.1. We consider an L = µ m by H = ν = . (cid:96) = . µ m and varying flexoelectric constants.29 T = 1 μ T = 5 μ T = 10 μ T = 50 (a) Deformed shape and electric potential [V] distribution of cantilever rods ofslenderness L / H =
20 under a point load of 20 nN, for di ff erent transversal flex-oelectric coe ffi cients µ T .
0 5 10
Load [nN]
15 20 r [ μ m ] a t S = L -1.5-1.0-0.50 E-B beam NL rod model μ T = 1 5 10 50 Numerical
Load [nN]
15 20 E [ V / μ m ] a t S = -150-100-500 E-B beam NL rod model μ T = 1 5 10 50 Numerical (b) Bending of a cantilever rod of slenderness L / H =
20 with varying transversalflexoelectric coe ffi cient µ T . The left plot shows the vertical displacement at theloaded end, and the right one shows the vertical electric field at the fixed end. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - E Y [V/μm]
10 20 30 40 μ T coefficient [nJ/Vm] L o a d [ n N ] -120-100-80-60-40-200E-B beamFit to numerical - - - - - (c) Countour plot the vertical electric field E at the fixed end of the rod of L / H =
20 as a function of the applied load and the transversal flexoelectriccoe ffi cient µ T . Load [nN] r [ μ m ] a t S = L -5-4-3-2-10 E-B beamNL rod model
L/H = 20 40 60 Numerical
Load [nN] E [ V / μ m ] a t S = -15-10-50 E-B beamNL rod model
L/H = 20 40 60 Numerical (d) Bending of a cantilever rod of µ T = / Vm with varying slenderness. Theleft panel shows the vertical displacement at the loaded end, and the right oneshows the vertical electric field at the fixed end.
Figure 4: Validation results for bending of open-circuit flexoelectric cantilever in sensor mode.The transversal flexoelectric coe ffi cient µ T in the legends is expressed in nJ / Vm = nC / m.30 Load [nN]
15 20 r [ μ m ] a t S = L -1.4-1.2-1.0-0.8-0.6-0.4-0.20 μ S = μ L =
0 , μ T = μ L = -10 , μ T = μ L =
0 , μ T = μ L =
10 , μ T = μ L = -10 , μ T =
10 0 10 -1-10 (a) Vertical displacement at the loaded end.
0 5 10
Load [nN]
15 20 E [ V / μ m ] a t S = -140-120-100-80-60-40-20020 (b) Vertical electric field at the fixed end. Figure 5: Electromechanical response of Neo-Hookean cantilever flexoelectric sensor under bend-ing, with di ff erent flexoelectric tensors (expressed in nJ / Vm).Fig. 5 represents the electromechanical response of the open circuit cantilever rod under pointload for varying flexoelectric constants µ L , µ T , µ S = {− , , } nJ / Vm . Fig. 5a shows the deflec-tion r of the loaded end, whereas Fig. 5b shows the vertical electric field E at the clamped end.For the sake of brevity, some combinations of flexoelectric tensors are omitted, since we found thatthe responses are analogous to the ones of other combinations as follows: r | µ = r | − µ ; (146a) E | µ = − E | − µ . (146b)From Fig. 5a, it is clear that flexoelectricity is always increasing the bending sti ff ness of therod. The largest sti ff ening is found with opposite µ T and µ L , followed by the case of vanishing µ L .On the contrary, the simulations with µ L ∼ µ T and the ones with vanishing µ T present a smallersti ff ening. In all cases, the e ff ect of the shear flexoelectric coe ffi cient µ S on bending sti ff ness ismuch smaller, and therefore less relevant.Fig. 5b shows the electric response of the rod at the clamped tip, revealing that all three flex-oelectric coe ffi cients are relevant here. Within the studied range, a larger flexoelectricity-inducedbending sti ff ness leads also to a larger electric field. However, in addition, the shear flexoelectrice ff ect µ S has a large influence on the electric field. In most cases, a non-vanishing µ S leads to asubstantial decrease in the reported electric field, which slightly depends also on the sign of µ S .The only case in which a non-vanishing µ S increases the electric field is the one where µ S is the only non-vanishing flexoelectric coe ffi cient. 31 a) Buckled rod geometry and resulting elec-tric potential [V]. Load [nN]
30 40 λ m i n ×10 -8 μ T = 0 μ T = 1 μ T = 5 μ T = 10 (b) Minimum eigenvalue λ min (cid:20) (cid:98) H χχ ( k ) (cid:21) . r [ μ m ] a t S = L / -2.5-2.0-1.5-1.0-0.50 0 10 20 Load [nN]
30 40 μ T = 0 μ T = 1 μ T = 5 μ T = 10 (c) Vertical displacement at the middle cross-section of the rod. Δ r [ μ m ] a t S = L -5-4-3-2-10 μ T = 0 μ T = 1 μ T = 5 μ T = 10 Load [nN]
30 40 (d) Horizontal displacement at the right endof the rod. E [ V / μ m ] a t S = L / Load [nN]
30 40 μ T = 0 μ T = 1 μ T = 5 μ T = 10 (e) Vertical electric field at the middle cross-section of the rod. Slenderness
L/H L o a d c r i t n u m / L o a d c r i t a n a μ T = 5 (f) Convergence of the buckling critical loadwith respect to L / H ratio. Figure 6: Force-controlled buckling of a flexoelectric rod of L / H =
60. Markers refer to thenumerical implementation and solid lines refer to the analytical nonlinear model for rods. Thetransversal flexoelectric coe ffi cient µ T is expressed in nJ / Vm.32 .2 Buckling of open-circuit flexoelectric rod under mechanical load
We now compress a slender flexoelectric rod ( L = µ m, H =
100 nm ) in open-circuit untilbuckling occurs, and also during the post-buckling stage. The left tip is clamped and a uniformhorizontal load is applied on the right cross-section, which can only move uniformly in axial di-rection, i.e. vertical displacement and rotation of the right end are prevented (see Fig. 3b). Weconsider an isotropic Saint-Venant–Kirchho ff model with Young’s modulus Y = . (cid:15) = . / V m and di ff erent transversal flexoelectric coe ffi cients: µ T = { , , , } nJ / Vm. The other material parameters are set to zero ( ν = µ L = µ S = (cid:96) = ff ect is not present yet since the rod is not bent, and hence the electric response iszero. Once the rod has buckled (see Fig. 6a), the vertical displacement at s = L / s = L (Fig. 6d) suddenly deviate from zero and evolve nonlinearly withrespect to the applied load. The flexoelectric e ff ect arises due to the curvature induced by buckling,leading to a measurable electric field at s = L /
2, which also evolves nonlinearly with applied load(Fig. 6e).The role of the magnitude of the flexoelectric coe ffi cient µ T is twofold. On the one hand, thecritical buckling load becomes larger with a larger µ T coe ffi cient, as suggested by the nonlinearrod model, cf. Eq. (108), for an e ff ectively sti ff er structure. Numerically, the precise value of thecritical buckling load is identified by the load at which the eigenvalue λ min (cid:20) (cid:98) H χχ ( k ) (cid:21) vanishes, asreported in Fig. 6b. On the other hand, the electric field at the post-buckling stage grows fasterwith a larger µ T coe ffi cient, which is also predicted by the nonlinear rod model, cf. Eq. (106b).Thus, the buckling-induced flexoelectric response is delayed but stronger when µ T is larger.We expect the agreement of the simplified rod model and the computational model to dete-riorate for thicker rods, and thus the assumptions of the rod model loose validity. In Fig. 6f weshow the e ff ect of the finite thickness of the rod on the buckling critical load by plotting the valuepredicted by the computational model normalized by that estimated by the nonlinear rod model fordi ff erent values of slenderness L / H . For all L / H values, the 1D nonlinear rod model overestimatesthe buckling load, as it provides a more constrained model. As expected, Fig. 6f shows that thebuckling critical load computed with the 2D computational model converges towards the approx-imated value given by the 1D nonlinear rod model as the slenderness L / H increases and thus the1D assumption is approached. We now consider a closed-circuit flexoelectric cantilever rod with Young’s modulus Y = . (cid:15) = . / V m, and dimensions L = µ m, H = µ m, which rolls up intoa circle upon electrical stimulus. The geometry and boundary conditions are depicted in Fig. 3c.33he left tip cross-section of the rod is clamped, while all other boundaries are traction-free. Theelectric potential at the top boundary is set to a certain non-zero value φ = V , and the bottomboundary is grounded ( φ = ff erence ∆ φ = V induces a transverse electric fieldacross the rod thickness, cf. Eq. (109), which triggers the flexoelectric and electrostrictive e ff ects,thereby generating a non-uniform strain that bends the rod, as shown in Fig. 8d. Depending onthe sign of the applied electric field the cantilever will bend upwards or downwards. This bendingactuator was first used by Bursian and Zaikovskii (1968) to experimentally demonstrate for the firsttime the flexoelectric e ff ect, which had been predicted theoretically by Mashkevich and Tolpygo(1957). Figure 7 shows the electromechanical response of an elastically isotropic Saint-Venant–Kirchho ff flexoelectric rod ( ν = l =
0) with the flexoelectric constants µ T = / Vm , µ L = µ S = / R (Fig. 7a) and the axial strain ζ (Fig. 7b) are captured very well by the closed-circuitflexoelectric rod model, where we have considered only the leading term in the expansions inEq. (126), up to a relatively large value of applied voltage V . Beyond this limit, the small strainsassumption of the 1D non-linear model loose validity. According to Eq. (126), the rod bends thanksto the flexoelectric coupling, and elongates mainly due to electrostriction, cf Section 4.2.1. Since the curvature is found to be uniform, cf. Eq. (127), the rod forms an arc of a circle, cf. Fig. 8d.Thus, a natural question that arises is which set of flexoelectric parameters achieve a fully-closedcircular shape more e ffi ciently (i.e. with a lower applied voltage). To address this question, we
0 0.2 0.4 0.6
Applied voltage V [kV] C u r v a t u r e [ mm - ] NL rod modelNumerical (a) Curvature R − ( V ). Applied voltage V [kV] A x i a l s t r a i n NL rod model Numerical (b) Axial strain ζ ( V ). Figure 7: Actuation of Saint-Venant–Kirchho ff cantilever rod with transversal flexoelectric coef-ficient µ T = / Vm. Numerically, the axial strain corresponds to the axial component of theGreen-Lagrangian strain tensor ( E ), whereas the value from the 1D model corresponds to itsTaylor approximation in Eq. (68), evaluated at X = ν = . (cid:96) = . µ m andvarying flexoelectric constants. To quantify the curvature of the rod relative to the curvature ofthe closed circle, we define the normalized curvature R − ( V ) = R − ( V ) / R − ◦ ( V ), where R − ◦ ( V ) = π/ ((1 + ζ ( V )) L ) is the curvature required to form a closed circular shape.Figure 8 shows the evolution of ζ ( V ), R − ( V ) and R − ( V ) for flexoelectric tensors with di ff erentcombinations of longitudinal ( µ L ), transversal ( µ T ) and shear ( µ S ) flexoelectric coe ffi cients. Thecases including a non-vanishing shear coe ffi cient are omitted, since the results do not change sig-nificantly, even when µ S is one order of magnitude larger than µ L or µ T . For the sake of brevity, thesimulations (i) with negative applied electric voltage V , and (ii) yielding negative curvatures, arealso omitted since the results are analogous to those simulations with (i) positive applied voltageand (ii) negative flexoelectric coe ffi cients, respectively, as ζ ( V ) | µ = ζ ( − V ) | µ = ζ ( V ) | − µ = ζ ( − V ) | − µ ; (147a) R − ( V ) | µ = − R − ( − V ) | µ = − R − ( V ) | − µ = R − ( − V ) | − µ ; (147b) R − ( V ) | µ = − R − ( − V ) | µ = − R − ( V ) | − µ = R − ( − V ) | − µ ; (147c)in accordance with Eqs. (126a) and (127).As expected, the axial strain of the rod (depicted in Fig. 8a) does not vary much with the dif-ferent flexoelectric parameters, since it is mainly a consequence of electrostriction. The curvature(Fig. 8b), instead, varies significantly for the di ff erent combinations of flexoelectric parameters.The dominant parameter is the transversal flexoelectric coe ffi cient µ T which leads to positive cur-vature, as shown in case B. The longitudinal flexoelectric coe ffi cient µ L is also relevant and leads tonegative curvature, as shown in case D. The largest response is found with positive µ T and negative µ L , as shown in case A. Finally, case C corresponds to positive µ L and µ T , and yields curvaturesinbetween cases B (purely transversal µ ) and D (purely longitudinal µ ).The normalized curvature is shown in Fig. 8c. For su ffi ciently large actuation, case A reaches R − >
1, which indicates that the actuator rolls up forming a closed circle. This process is shownin Fig. 8d, where the deformed configuration and electric potential distribution within the rod isdepicted at di ff erent applied voltages. We also show in Fig. 8e the resulting polarization field oncethe circle is formed, which remains normal to the bent rod. In the previous example, the rod undergoes elongation upon electrical actuation mainly due toelectrostriction. In this Section, we present a similar setup where the right tip is also clamped, asshown in Fig. 3d. In this case, an axial compressive force is expected at the clamped ends sincethe elongation of the rod is prevented. Restricting Eq. (134) in pre-buckling stage, the axial forcegrows quadratically with the applied voltage and, for a large enough applied (critical) voltage V cr ,cf. Eq. (144), a mechanical instability is reached, inducing buckling of the rod.35 pplied voltage V [kV] A x i a l s t r a i n ABC D -10 10 0-10 0 0 0 10 010 10 0
ABCD μ L μ T μ S (a) Axial strain ζ ( V ) C u r v a t u r e [ mm - ] ABCD
Applied voltage V [kV] -10 10 0-10 0 0 0 10 010 10 0 ABCD μ L μ T μ S (b) Curvature R − ( V ) Applied voltage V [kV] N o r m a l i z e d C u r v a t u r e ABCD -10 10 0-10 0 0 0 10 010 10 0
ABCD μ L μ T μ S (c) Normalized curvature R − ( V ) .
0 1 . ϕ [kV]
99 1 . (d) Deformed configuration and electric potential distributionin case A upon increasing voltage V [kV], indicated by thenumber at the free end. . (e) Distribution of polar-ization field in case A at V = . Figure 8: Actuation of Neo-Hookean cantilever rod with di ff erent flexoelectric tensors (expressedin nJ / Vm)Figure 9 shows numerical simulations of a flexoelectric Saint-Venant–Kirchho ff rod ( ν = (cid:96) = L = µ m, H = µ m, with Young’s modulus Y = . (cid:15) = . / V m and transversal flexoelectric coe ffi cient µ T = / Vm ( µ L = µ S = ϕ [kV] (a) Buckled shape and electric poten-tial distribution of the L / H = . Applied voltage V [kV] -1.4-1.2-1.0-0.8-0.6-0.4-0.20 V e r t i c a l d e f l e c t i o n [ μ m ] a t S = L / L / H =
NumericalNL rod model
20 40 60 (b) Vertical deflection r at the centerof the rod. Applied voltage [kV] V ×10 -3 Numerical NL rod model = 20 40 60 A x i a l s t r a i n a t S = L / L/H (c) Axial strain as a function of slen-derness. λ m i n Applied voltage V [kV] L / H =
NumericalAnalytical V cr
20 40 60 (d) Normalized minimum eigenvalue λ min (cid:20) (cid:98) H χχ ( k ) (cid:21) and critical voltage ofanalytical 1D model. C r i t i c a l e l e c t r i c f i e l d [ V / μ m ] L / H = 20 L / H = 40 L / H = 60 NL rod model (e) Critical electric field as a functionof slenderness. Figure 9: Actuation of Saint-Venant–Kirchho ff clamped-clamped rod with transversal flexoelectriccoe ffi cient µ T = / Vm and varying slenderness. In (d), ¯ λ min = λ min ( n DOF / n ) ,where n DOF isthe number of degrees of freedom of each simulation, and n =
312 is an arbitrary normalizationconstant, chosen such that ¯ λ min (0) ≈
1. 37
Conclusions and directions of future work
We have developed the material form of the balance equations for dielectric elastomers, includingthe flexoelectric e ff ect. Unlike previously considered models of the flexoelectric coupling, herewe formulate our model in terms of polarization, strain gradients and flexocoupling tensor in afully material frame. As a result, our formulation is objective by construction, and the flexocou-pling tensor has the same symmetries as that used in linearized theories. After partial Legendretransform, the equations are written in terms of the electric potential and the displacement fieldas a fourth order unconstrained system of partial di ff erential equations, which is convenient forfinding numerical and analytical solutions. A numerical implementation of the theory is developedusing open B-spline basis of su ffi cient smoothness on a uniform Cartesian grid, enabling robustsimulations deep into the nonlinear regime, for very large deformations, and including mechan-ical instabilities (Yvonnet and Liu, 2017). On the other hand, analytical closed-form solutionsare derived for open- and closed-circuit nonlinear extensible flexoelectric rods under bending andbuckling. Direct comparison of this model with direct numerical simulations of the full modelshows excellent agreement well into the nonlinear regime in conditions where the rod theory isexpected to apply. The analytical rod theory serves both as a means of validation of our nonlinearsimulations, and as fast and simple model to analyze and design nonlinear flexoelectric devices.The current model could be easily extended in several ways. For instance, rather than homo-geneous electric Neumann boundary conditions on the free surfaces, it may be more realistic todirectly model the surrounding medium as a dielectric when considering soft materials materialswith relatively low dielectric constant (Yvonnet and Liu, 2017, Thai et al., 2018). Our model canbe extended to account for converse flexoelectricity (Lifshitz and Landau, 1951, Sharma et al.,2010, Landau and Lifshitz, 2013), for polarization gradient dielectricity (Mindlin, 1968), for ma-terial incompressibility, and coupled with flexible discretization methods, e.g. based on immersedboundaries (Codony et al., 2019), to model domains of general, and possibly complex, geometrythat might enhance field gradients. Acknowledgments
This work was supported by the Generalitat de Catalunya (“ICREA Academia” award for ex-cellence in research to I.A., and Grant No. 2017-SGR-1278), and the European Research Coun-cil (StG-679451 to I.A.). CIMNE is recipient of a Severo Ochoa Award of Excellence from theMINECO.
Appendix A Material characterization
The material is fully characterized by specifying the elastic energy density Ψ Elast ( C ) and the mate-rial tensors of flexoelectricity µ and strain gradient elasticity h .38 sotropic Saint-Venant–Kirchho ff model. It corresponds to the extension of the linear isotropic elastic material model to the non-linearregime, and depends on the Lam´e parameters λ = Y ν/ (1 + ν )(1 − ν ) and µ = Y / + ν ) as follows: Ψ Elast ( C ) = λ E )] + µ Tr( E ) , (A.1a) ∂ Ψ Elast ( C ) ∂ C IJ = λ E )] δ IJ + µ E IJ , (A.1b) ∂ Ψ Elast ( C ) ∂ C IJ C KL = λ δ IJ δ KL + µ δ IK δ JL . (A.1c) Isotropic Neo-Hookean model
The Neo-Hookean model is adequate for describing nonlinear stress-strain behavior of cross-linkedpolymers at moderate strains. It is mathematically defined as Ψ Elast ( C ) = λ (cid:2) log( J ) (cid:3) + µ C ) − , (A.2a) ∂ Ψ Elast ( C ) ∂ C IJ = λ J ) C − IJ + µ (cid:16) δ IJ − C − IJ (cid:17) , (A.2b) ∂ Ψ Elast ( C ) ∂ C IJ C KL = λ C − IJ C − KL + (cid:2) µ − λ log( J ) (cid:3) (cid:16) C − IK C − JL + C − IL C − JK (cid:17) . (A.2c) Flexoelectricity tensor µ . The cubic flexoelectric tensor depends on the longitudinal µ L , transversal µ T and shear µ S param-eters (Le Quang and He, 2011, Codony et al., 2019). In the Cartesian axes, it takes the followingform: µ LIJK = µ L , for L = I = J = K ,µ T , for I = J (cid:44) K = L ,µ S , for L = I (cid:44) J = K or L = J (cid:44) I = K , . (A.3) Strain gradient elasticity tensor h.
We consider an isotropic simplified strain gradient elasticity tensor (Altan and Aifantis, 1997),which depends on λ , µ and the length scale (cid:96) in the following form: h IJKLMN = ( λδ IJ δ LM + µδ IL δ JM ) (cid:96) δ KN . (A.4)39 ppendix B Second variation of the enthalpy functional The second variation of the enthalpy functional, required in our solution method, is given by δ Π [ χ , φ ; δ χ , δφ ; ∆ χ , ∆ φ ] = δ (cid:0) R [ χ , φ ; δ χ , δφ ] (cid:1) [ ∆ χ , ∆ φ ] = ∂ R [ χ , φ ; δ χ , δφ ] ∂ χ ∆ χ + ∂ R [ χ , φ ; δ χ , δφ ] ∂φ ∆ φ = (cid:90) Ω (cid:40) δ E IJ ∆ E KL (cid:32) ∂ ¯ Ψ Elast ( C ) ∂ C IJ ∂ C KL (cid:33) + (cid:32) ∂ Ψ Elast ( C ) ∂ C IJ (cid:33) ( ∆ δ ) E IJ + h IJKLMN δ (cid:101) E IJK ∆ (cid:101) E IJK + (cid:16) h IJKLMN (cid:101) E LMN (cid:17) ( ∆ δ ) (cid:101) E IJK − (cid:15) JC − MF δ E F ∆ E M + (cid:15) J C MFIJ E F (cid:32) E M ( ∆ δ ) E IJ + δ E IJ ∆ E M + δ E M ∆ E IJ (cid:33) + (cid:15) J (cid:101) C MFIJKL E M E F δ E IJ ∆ E KL − µ FABK JC − MF (cid:16) E M ( ∆ δ ) (cid:101) E ABK + δ E M ∆ (cid:101) E ABK + δ (cid:101) E ABK ∆ E M (cid:17) + µ FABK J C MFIJ (cid:16)(cid:101) E ABK ( δ E IJ ∆ E M + δ E M ∆ E IJ ) + E M (cid:16) δ E IJ ∆ (cid:101) E ABK + δ (cid:101) E ABK ∆ E IJ (cid:17) + E M (cid:101) E ABK ( ∆ δ ) E IJ (cid:17) + µ FABK J (cid:101) C MFIJPQ E M (cid:101) E ABK δ E IJ ∆ E PQ (cid:41) d Ω = (cid:90) Ω (cid:40)(cid:98) S IJ ( ∆ δ ) E IJ + (cid:101) S IJK ( ∆ δ ) (cid:101) E IJK + (cid:16) A Elast
IJKL + A Diele
IJKL + A Flexo
IJKL (cid:17) δ E IJ ∆ E KL + (cid:101) A SGEla
IJKLMN δ (cid:101) E IJK ∆ (cid:101) E LMN + (cid:101) A Flexo
IJKLM (cid:16) δ E IJ ∆ (cid:101) E KLM + δ (cid:101) E KLM ∆ E IJ (cid:17) + B Diele IJ ( δ E I ∆ E J ) + (cid:16) C Diele
IJK + C Flexo
IJK (cid:17) ( δ E IJ ∆ E K + δ E K ∆ E IJ ) + (cid:101) C Flexo
IJKL (cid:16) δ E L ∆ (cid:101) E IJK + δ (cid:101) E IJK ∆ E L (cid:17) (cid:41) d Ω , (B.1)40here ∆ χ and ∆ φ are variations of χ and φ , respectively, and ∆ E L (cid:66) − ∂ ( ∆ φ ) ∂ X L , (B.2a) ∆ F iI (cid:66) ∂ ( ∆ x i ) ∂ X I , (B.2b) ∆ (cid:101) F iIJ (cid:66) ∂ ( ∆ x i ) ∂ X I ∂ X J , (B.2c) ∆ E IJ = ∆ C IJ (cid:66) symm IJ ( ∆ F kI F kJ ) , (B.2d) ∆ (cid:101) E IJK = ∆ (cid:101) C IJK (cid:66) symm IJ (cid:16) ∆ F kI (cid:101) F kJK + F kI ∆ (cid:101) F kJK (cid:17) , (B.2e)( ∆ δ ) E IJ = ( ∆ δ ) C IJ (cid:66) symm IJ ( ∆ F kI δ F kJ ) , (B.2f)( ∆ δ ) (cid:101) E IJK = ( ∆ δ ) (cid:101) C IJK (cid:66) symm IJ (cid:16) ∆ F kI δ (cid:101) F kJK + δ F kI ∆ (cid:101) F kJK (cid:17) . (B.2g)The material tensors in the right hand side of Eq. (B.1) are defined as follows: A Elast
IJKL ( C ) (cid:66) ∂ ¯ Ψ Elast ∂ E IJ ∂ E KL (B.3a) A Diele
IJKL ( C , E ) (cid:66) ∂ ¯ Ψ Diele ∂ E IJ ∂ E KL = J (cid:101) C MFIJKL E M (cid:15) E F (B.3b) A Flexo
IJKL ( C , (cid:101) C , E ) (cid:66) ∂ ¯ Ψ Flexo ∂ E IJ ∂ E KL = J (cid:101) C MFIJKL E M µ FABC (cid:101) E ABC (B.3c) (cid:101) A SGEla
IJKLMN (cid:66) ∂ ¯ Ψ SGEla ∂ (cid:101) E IJK ∂ (cid:101) E LMN = h IJKLMN , (B.3d) (cid:101) A Flexo
IJKLM ( C , E ) (cid:66) ∂ ¯ Ψ Flexo ∂ E IJ ∂ (cid:101) E KLM = J C ABIJ µ BKLM E A (B.3e) B Diele IJ ( C ) (cid:66) ∂ ¯ Ψ Diele ∂ E I ∂ E J = − (cid:15) JC − IJ (B.3f) C Diele
IJK ( C , E ) (cid:66) ∂ ¯ Ψ Diele ∂ E IJ ∂ E K = (cid:15) J C KMIJ E M (B.3g) C Flexo
IJK ( C , (cid:101) C ) (cid:66) ∂ ¯ Ψ Flexo ∂ E IJ ∂ E K = µ MABC J C KMIJ (cid:101) E ABC (B.3h) (cid:101) C Flexo
IJKL ( C ) (cid:66) ∂ ¯ Ψ Flexo ∂ (cid:101) E IJK ∂ E L = − µ MIJK JC − ML (B.3i)41he tensor (cid:101) C in Eq. (B.1) is defined as (cid:101) C ABCDEF (cid:66) J ∂ ( J C ABCD ) ∂ C EF = ( D ACBDEF + D BDACEF + D ADBCEF + D BCADEF − D ABCDEF − D CDABEF ) , (B.4)where D ABCDEF (cid:66) C − AB (cid:16) C − CD C − EF − C − CE C − DF − C − CF C − DE (cid:17) . Appendix C Analytical solutions for the displacement and theelectric field in flexoelectric rods under bending
Following Bisshopp and Drucker (1945), Eq. (88) is integrated as12 (cid:32) d θ d S (cid:33) + β (sin θ max − sin θ ) = , (C.1)where θ ( L ) = θ max ≤ S = − d θβ √ θ − sin θ max ) , (C.2)since θ ≤ θ/ d S ≤ L = (cid:90) L d S = (cid:90) θ ( L ) θ (0) d S d θ d θ = (cid:90) θ max d θβ √ θ − sin θ max ) , (C.3)and thus β L = (cid:90) θ max d θ √ θ − sin θ max ) . (C.4)In order to evaluate this integral, let us assumesin θ max = − p , sin θ = − p sin ψ, ψ ∈ [ ψ , π , (C.5)with ψ = sin − p √ = sin − (cid:32) √ − sin θ max (cid:33) . (C.6)Using cos θ = (cid:112) − sin θ = p sin ψ (cid:112) − p sin ψ , we obtaind θ = − p cos ψ (cid:113) − p sin ψ d ψ, (C.7)42nd substituting in Eq. (C.4) yields β L = (cid:90) π/ ψ d ψ (cid:113) − p sin ψ , (C.8)which can be written as β L = F ( p ) − F ( p , ψ ) , (C.9)where F ( p ) = (cid:90) π/ (cid:113) − p sin ψ d ψ, and F ( p , ψ ) = (cid:90) ψ (cid:113) − p sin ψ d ψ. (C.10)are the complete and incomplete elliptical integrals of the first kind, respectively (Jahnke, 1945).Hence, for a given value of θ max , β can be determined from Eq. (C.9) using Eqs. (C.5) and (C.6),and the corresponding applied vertical load producing the rotation θ max at the free end is then N = ¯ Y I e ff β . For a given N , the problem is thus solved by the shooting method.Using Eq. (61), the vertical displacement of the rod is r ( S ) = (cid:90) S (1 + ζ ) sin θ d ˜ S ≈ (cid:90) θ sin θ d θβ √ θ max − sin θ ) = (cid:90) ψψ p sin ˜ ψ − β (cid:113) − p sin ˜ ψ d ˜ ψ = β (cid:2) F ( p , ψ ) − F ( p , ψ ) (cid:3) + β (cid:104) ˜ E ( p , ψ ) − ˜ E ( p , ψ ) (cid:105) , (C.11)where ˜ E ( p ) = (cid:90) π/ (cid:113) − p sin ψ d ψ, and ˜ E ( p , ψ ) = (cid:90) ψ (cid:113) − p sin ψ d ψ, (C.12)are the complete and incomplete elliptical integrals of the second kind, respectively (Jahnke, 1945).Thus, the deflection of the rod at its loaded end is r ( L ) = L + β (cid:104) ˜ E ( p , ψ ) − ˜ E ( p ) (cid:105) . (C.13)Finally, the vertical electric field is computed from Eq. (74) as E ( S ) = µ T (cid:15) θ (cid:48) = − βµ T (cid:15) (cid:112) θ − sin θ max ) = − µ T (cid:15) (cid:114) N ¯ Y I e ff (sin θ − sin θ max ) . (C.14)Therefore, the electric field at the fixed end is E (0) = − µ T (cid:15) (cid:114) N ¯ Y I e ff sin | θ max | . (C.15)43 ppendix D Analytical solutions for displacement and electricfield in flexoelectric rods under compressive axialload Integration of Eq. (104) yields 12 (cid:32) d θ d S (cid:33) − β (cos θ − cos θ max ) = , (D.1)where we assume upward buckling without loss of generality, and θ ( L / = θ max >
0. Equivalently,d S = d θβ √ θ − cos θ max ) . (D.2)Since the right end of the rod is allowed to move horizontally under the action of the compressiveload, the length of the rod is assumed to remain approximately unaltered after buckling. Hence,using Eq. (D.2), L = θ ( L ) (cid:90) θ (0) d S d θ d θ = θ max (cid:90) d θβ √ θ − cos θ max ) , (D.3)and thus β L = θ max (cid:90) (cid:114) sin θ max − sin θ θ. (D.4)To compute this integral, we definesin θ max = p , sin θ = p sin ψ, ψ ∈ [0 , π . (D.5)Hence, β L = F ( p ) = F (cid:32) sin θ max (cid:33) , (D.6)where again F ( p ) is the complete elliptical integral of the first kind, see Eq. (C.10). So, for a givenload N , θ max is determined by the shooting method, i.e. by giving values to θ max and computing thecorresponding loading parameter β from Eq. (D.6) until the target β = (cid:112) N / ¯ Y I e ff is reached.44imilarly, the change in the horizontal displacement, ∆ r , can be evaluated by the di ff erence ofactual length ( L ) and the length projected over axial direction upon buckling as ∆ r ≈ L − θ ( L ) (cid:90) θ (0) cos θ d S d θ d θ = L − θ max (cid:90) θβ (cid:114) sin θ max − sin θ θ = L − β π/ (cid:90) − p sin ψ (cid:113) − p sin ψ d ψ = L − E ( p ) − F ( p ) β = F ( p ) − ˜ E ( p )] β , (D.7)where we have used (cid:112) sin ( θ max / − sin ( θ/ = p cos ψ , cos θ/ cos ( θ/ = − p sin ψ (cid:113) − p sin ψ , andEq. (D.6), and again ˜ E ( p ) is the complete elliptical integral of the second kind, see Eq. (C.12).Since, the deformations in the half-rod are antisymmetric with respect to S = L /
4, we split thevertical deflection into two parts. Hence, assuming that the rod buckles upwards without loss ofgenerality, S ∈ (cid:20) , L (cid:21) : r ( S ) ≈ (cid:90) S sin θ ( u ) d u = (cid:90) θ sin γ d γβ (cid:112) γ − cos θ max ) = (cid:90) ψ p sin ξ d ξβ = p β (1 − cos ψ ) , ψ ∈ [0 , π S ∈ (cid:20) L , L (cid:21) : r ( S ) ≈ p β − (cid:90) θθ max sin γ d γβ (cid:112) γ − cos θ max ) = p β + (cid:90) ψ π p sin (cid:16) π − ξ (cid:17) d ξβ = p β (cid:18) + cos (cid:18) π − ψ (cid:19)(cid:19) , ψ ∈ [0 , π , (D.8b)(D.8c)where we have used (cid:112) sin ( θ max / − sin ( θ/ = p cos ψ , and sin θ/ cos ( θ/ = p sin ψ . Finally,the electric field can be evaluated as S ∈ (cid:20) , L (cid:21) : E ( S ) = µ T (cid:15) θ (cid:48) = µ T β(cid:15) (cid:112) θ − cos θ max ) = µ T (cid:15) (cid:114) N ¯ Y I e ff (cos θ − cos θ max ) , θ ∈ (cid:2) , θ max (cid:3) (D.9a) S ∈ (cid:20) L , L (cid:21) : E ( S ) = − µ T (cid:15) (cid:114) N ¯ Y I e ff (cid:18) cos (cid:18) θ (cid:18) L − S (cid:19)(cid:19) − cos θ max (cid:19) , θ ∈ (cid:2) , θ max (cid:3) (D.9b)45herefore, the vertical deflection and electric field at the center of the rod are r (cid:18) L (cid:19) = β sin (cid:32) θ max (cid:33) , (D.10) E (0) = − E (cid:18) L (cid:19) = E ( L ) = µ T (cid:15) (cid:114) N ¯ Y I e ff (1 − cos θ max ) . (D.11) Appendix E Analytical solutions for displacement and voltagein flexoelectric rods under transversal voltage ac-tuation
Similarly to Appendix D, integration of the moment balance Eq. (136) yields12 (cid:32) d θ d S (cid:33) − (1 + ζ ) ˜ β (cos θ − cos θ max ) = , (E.1)where θ ( L / = θ max and upon integration L = (cid:90) L / d S = θ ( L ) (cid:90) θ (0) d S d θ d θ = (cid:112) + ζ θ max (cid:90) d θ ˜ β √ θ − cos θ max ) = F ( p )˜ β (cid:112) + ζ . (E.2)Thus ˜ β (cid:112) + ζ L = F ( p ) = F (cid:32) sin θ max (cid:33) . (E.3)In this case, the right end of the rod is clamped and thus the length of the rod after buckling isunknown, but its projection on the horizontal axis is the undeformed length L , therefore with thehelp of constraint Eq. (133) L = (cid:90) L / d r = θ max (cid:90) (1 + ζ ) cos θ d S d θ d θ = (cid:112) (1 + ζ ) (cid:104) E ( p ) − F ( p ) (cid:105) ˜ β . (E.4)Therefore, by using Eqs. (E.3) and (E.4) ζ = F ( p )2 ˜ E ( p ) − F ( p ) − . (E.5)Once ζ is known, ˜ β can be evaluated using Eq. (E.3) for any θ max .46ow, similar to Appendix D, S ∈ (cid:20) , L (cid:21) : r ( S ) = − (cid:90) θ (1 + ζ ) sin θ d θ ˜ β (cid:112) + ζ √ θ − cos θ max ) = − (cid:90) ψ p (cid:112) + ζ sin ψ d ψ ˜ β = − p (cid:112) + ζ (1 − cos ψ )˜ β , S ∈ (cid:20) L , L (cid:21) : r ( S ) = − (cid:112) + ζ ˜ β (cid:32)(cid:90) θ max sin θ d θ √ θ − cos θ max ) − (cid:90) θθ max sin θ d θ √ θ − cos θ max ) (cid:33) = − p (cid:112) + ζ (1 + cos ψ )˜ β , (E.6)with sin θ = p sin ψ . Hence, the deflection at the center of the rod and the curvature at the left endfor downward buckling are r (cid:18) L (cid:19) = − p (cid:112) + ζ ˜ β , (E.7) θ (cid:48) (0) = ˜ β (cid:112) + ζ )(1 − cos θ max ) . (E.8) References
Amir Abdollahi and Irene Arias. Constructive and destructive interplay between piezoelectricityand flexoelectricity in flexural sensors and actuators.
Journal of Applied Mechanics , 82(12),2015. URL https://doi.org/10.1115/1.4031333 .Amir Abdollahi, Christian Peco, Daniel Mill´an, Marino Arroyo, and Irene Arias. Computationalevaluation of the flexoelectric e ff ect in dielectric solids. Journal of Applied Physics , 116(9):093502, 2014. URL https://doi.org/10.1063/1.4893974 .Amir Abdollahi, Daniel Mill´an, Christian Peco, Marino Arroyo, and Irene Arias. Revisiting pyra-mid compression to quantify flexoelectricity: A three-dimensional simulation study.
Phys. Rev.B , 91:104103, 2015a. URL https://doi.org/10.1103/PhysRevB.91.104103 .Amir Abdollahi, Christian Peco, Daniel Mill´an, Marino Arroyo, Gustau Catalan, and Irene Arias.Fracture toughening and toughness asymmetry induced by flexoelectricity.
Phys. Rev. B , 92:094101, 2015b. URL https://doi.org/10.1103/PhysRevB.92.094101 .F. Ahmadpoor and P. Sharma. Flexoelectricity in two-dimensional crystalline and biological mem-branes.
Nanoscale , 7:16555, 2015. 47. Ahmadpoor, Q. Deng, L. P. Liu, and P. Sharma. Apparent flexoelectricity in lipid bilayer mem-branes due to external charge and dipolar distributions.
Physical Review E , 88:050701, 2013.BS Altan and EC Aifantis. On some aspects in the special theory of gradient elasticity.
Journalof the Mechanical Behavior of Materials , 8(3):231–282, 1997. URL https://doi.org/10.1515/JMBM.1997.8.3.231 .Li Anqing, Zhou Shenjie, Qi Lu, and Chen Xi. A flexoelectric theory with rotation gradient e ff ectsfor elastic dielectrics. Modelling and Simulation in Materials Science and Engineering , 24(1):015009, 2015. URL https://doi.org/10.1088/0965-0393/24/1/015009 .S.S. Antman.
Nonlinear Problems of Elasticity . Springer-Verlag, New York, 1995. ISBN9780387276496. URL https://doi.org/10.1007/0-387-27649-1 .G. Barbero, I. Dozov, J. F. Palierne, and G. Durand. Order electricity and surface orientation innematic liquid crystals.
Physical Review Letters , 56(19):2056–2059, 1986.S Baroudi and F Najar. Dynamic analysis of a nonlinear nanobeam with flexoelectric actua-tion.
Journal of Applied Physics , 125(4):044503, 2019. URL https://doi.org/10.1063/1.5057727 .Sivapalan Baskaran, Xiangtong He, Qin Chen, and John Y Fu. Experimental studies on the directflexoelectric e ff ect in α -phase polyvinylidene fluoride films. Applied Physics Letters , 98(24):242901, 2011. URL https://doi.org/10.1063/1.3599520 .Sivapalan Baskaran, Xiangtong He, Yu Wang, and John Y Fu. Strain gradient induced electric po-larization in α -phase polyvinylidene fluoride films under bending conditions. Journal of AppliedPhysics , 111(1):014109, 2012. URL https://doi.org/10.1063/1.3673817 .S Bauer and F Bauer. Piezoelectric polymers and their applications. In
Piezoelectricity , pages157–177. Springer, 2008. URL https://doi.org/10.1007/978-3-540-68683-5_6 .KE Bisshopp and DC Drucker. Large deflection of cantilever beams.
Quarterly of Applied Math-ematics , 3(3):272–275, 1945. URL .Lance Breger, Takeo Furukawa, and Eiichi Fukada. Bending piezoelectricity in polyvinylidenefluoride.
Japanese Journal of Applied Physics , 15(11):2239, 1976. URL https://doi.org/10.1143/JJAP.15.2239 .´E. V. Bursian and N. N. Trunov. Nonlocal piezoelectric e ff ect. Sov. Phys. Solid State ,16(4):760 – 762, 1974. URL .J M Bursian and O I Zaikovskii. Changes in curvature of a ferroelectric film due to polarization.
Soviet Physics Solid State , 10(5):1121–1124, 1968.48ojca ˇCepiˇc, Barbara Rovˇsek, and Boˇstjan ˇZekˇs. Flexoelectrically induced polarization in polarsmectic films.
Ferroelectrics , 244(1):59–66, 2000.Baojin Chu and DR Salem. Flexoelectricity in several thermoplastic and thermosetting poly-mers.
Applied Physics Letters , 101(10):103905, 2012. URL https://doi.org/10.1063/1.4750064 .David Codony, Onofre Marco, Sonia Fern´andez-M´endez, and Irene Arias. An immersed boundaryhierarchical b-spline method for flexoelectricity.
Computer Methods in Applied Mechanics andEngineering , 354:750–782, 2019. URL https://doi.org/10.1016/j.cma.2019.05.036 .C. de Boor.
A Practical Guide to Splines . Applied Mathematical Sciences. Springer New York,2001. ISBN 9780387953663. URL .P. G. de Gennes and J. Prost.
The Physics of Liquid Crystals . Number 83 in International Series ofMonographs on Physics. Oxford Science Publications, second edition edition, 1993.Feng Deng, Qian Deng, Wenshan Yu, and Shengping Shen. Mixed finite elements for flexoelectricsolids.
Journal of Applied Mechanics , 84(8), 2017. URL https://doi.org/10.1115/1.4036939 .Feng Deng, Qian Deng, and Shengping Shen. A three-dimensional mixed finite element for flexo-electricity.
Journal of Applied Mechanics , 85(3), 2018. URL https://doi.org/10.1115/1.4038919 .Qian Deng, Mejdi Kammoun, Alper Erturk, and Pradeep Sharma. Nanoscale flexoelectric energyharvesting.
International Journal of Solids and Structures , 51(18):3218–3225, 2014a. URL https://doi.org/10.1016/j.ijsolstr.2014.05.018 .Qian Deng, Liping Liu, and Pradeep Sharma. Electrets in soft materials: Nonlinearity, size e ff ects,and giant electromechanical coupling. Physical Review E , 90(1):012603, 2014b. URL https://doi.org/10.1103/PhysRevE.90.012603 .Qian Deng, Liping Liu, and Pradeep Sharma. Flexoelectricity in soft materials and biologicalmembranes.
Journal of the Mechanics and Physics of Solids , 62:209–227, 2014c. URL https://doi.org/10.1016/j.jmps.2013.09.021 .A. Derzhanski, A. G. Petrov, A. T. Todorov, and K. Hristova. Flexoelectricity of lipid bilayers.
Liquid Crystals , 7(3):439–449, 1990.A. F. Devonshire. Theory of barium titanate. part i.
Philosophical Magazine , 40:1040, 1949.A. F. Devonshire. Theory of barium titanate. part ii.
Philosophical Magazine , 42:1065, 1951.A. F. Devonshire. Theory of ferroelectrics.
A quarterly Supplement of the Philosophical Magazine ,3(10):85–130, 1954. 49 Dorfmann and RW Ogden. Nonlinear electroelasticity.
Acta Mechanica , 174(3-4):167–183,2005. URL https://doi.org/10.1007/s00707-004-0202-2 .Luis Dorfmann and Ray W Ogden.
Nonlinear theory of electroelastic and magnetoelastic interac-tions , volume 1. Springer, 2014. URL https://doi.org/10.1007/978-1-4614-9596-3 .Luis Dorfmann and Ray W Ogden. Nonlinear electroelasticity: material properties, continuumtheory and applications.
Proceedings of the Royal Society A: Mathematical, Physical and En-gineering Sciences , 473(2204):20170311, 2017. URL https://doi.org/10.1098/rspa.2017.0311 .John Y Fu, Wenyi Zhu, Nan Li, and L Eric Cross. Experimental studies of the converse flexo-electric e ff ect induced by inhomogeneous electric field in a barium strontium titanate composi-tion. Journal of Applied Physics , 100(2):024112, 2006. URL https://doi.org/10.1063/1.2219990 .L.T. Gao, X-Q Feng, Y-J Yin, and H. Gao. An electromechanical liquid crystal model of vesicles.
Journal of the Mechanics and Physics of Solids , 56:2844–2862, 2008.Hamid Ghasemi, Harold S Park, and Timon Rabczuk. A level-set based iga formulation for topol-ogy optimization of flexoelectric materials.
Computer Methods in Applied Mechanics and En-gineering , 313:239–258, 2017. URL https://doi.org/10.1016/j.cma.2016.09.029 .Hamid Ghasemi, Harold S Park, and Timon Rabczuk. A multi-material level set-based topologyoptimization of flexoelectric composites.
Computer Methods in Applied Mechanics and Engi-neering , 332:47–62, 2018. URL https://doi.org/10.1016/j.cma.2017.12.005 .Prakhar Gupta and Ajeet Kumar. E ff ect of material nonlinearity on spatial buckling of nanorodsand nanotubes. Journal of Elasticity , 126(2):155–171, 2017. URL https://doi.org/10.1007/s10659-016-9586-1 .Ali R Hadjesfandiari. Size-dependent piezoelectricity.
International Journal of Solids and Struc-tures , 50(18):2781–2791, 2013. URL https://doi.org/10.1016/j.ijsolstr.2013.04.020 .Khader M Hamdia, Hamid Ghasemi, Xiaoying Zhuang, Naif Alajlan, and Timon Rabczuk. Sen-sitivity and uncertainty analysis for flexoelectric nanostructures.
Computer Methods in AppliedMechanics and Engineering , 337:95–109, 2018. URL https://doi.org/10.1016/j.cma.2018.03.016 .John Harden, Badel Mbanga, Nandor ´Eber, Katalin Fodor-Csorba, Samuel Sprunt, James T Glee-son, and Antal Jakli. Giant flexoelectricity of bent-core nematic liquid crystals.
Physical reviewletters , 97(15):157802, 2006. URL https://doi.org/10.1103/PhysRevLett.97.157802 .50iawang Hong and David Vanderbilt. First-principles theory of frozen-ion flexoelectricity.
PhysicalReview B , 84(18):180101, 2011. URL https://doi.org/10.1103/PhysRevB.84.180101 .ShuLing Hu and ShengPing Shen. Variational principles and governing equations in nano-dielectrics with the flexoelectric e ff ect. Science China Physics, Mechanics and Astronomy , 53(8):1497–1504, 2010. URL https://doi.org/10.1007/s11433-010-4039-5 .Shujin Huang, Lu Qi, Wenbin Huang, Longlong Shu, Shenjie Zhou, and Xiaoning Jiang. Flex-oelectricity in dielectrics: Materials, structures and characterizations.
Journal of AdvancedDielectrics , 8(02):1830002, 2018. URL https://doi.org/10.1142/S2010135X18300025 .VL Indenbom, EB Loginov, and MA Osipov. Flexoelectric e ff ect and structure of crystals. Kristal-lografiya , 28:1157–1162, 1981a.VL Indenbom, EB Loginov, and MA Osipov. Flexoelectric e ff ect and crystal-structure. Kristallo-grafiya , 26(6):1157–1162, 1981b.Eugene Jahnke. Tables of functions with formulae and curves.
New York: Dover Publications,—c1945, 4th ed. , 1945. URL https://ui.adsabs.harvard.edu/abs/1945tfwf.book.....J .S. A. Jewell. Living systems and liquid crystals.
Liquid Crystals , 38(11-12):1699–1714, 2011.Xiaoning Jiang, Wenbin Huang, and Shujun Zhang. Flexoelectric nano-generator: Materials, struc-tures and devices.
Nano Energy , 2(6):1079–1092, 2013. URL https://doi.org/10.1016/j.nanoen.2013.09.001 .Sh M Kogan. Piezoelectric e ff ect during inhomogeneous deformation and acoustic scattering ofcarriers in crystals. Soviet Physics-Solid State , 5(10):2069–2070, 1964.Sana Krichen and Pradeep Sharma. Flexoelectricity: a perspective on an unusual electromechani-cal coupling.
Journal of Applied Mechanics , 83(3), 2016. URL https://doi.org/10.1115/1.4032378 .W. Kuczynski and J. Ho ff mann. Determination of piezoelectric and flexoelectric polarization inferroelectric liquid crystals. Physical Review E , 72(4):041701, 2005.S.T. Lagerwall and I. Dahl. Ferroelectric liquid crystals.
Molecular Crystals and Liquid Crystals ,114(1-3):151–187, 1984.Lev Davidovich Landau and Evgenii Mikhailovich Lifshitz.
Course of theoretical physics . Else-vier, 2013. URL https://books.google.es/books?id=LuBbAwAAQBAJ .M Lax and DF Nelson. Maxwell equations in material form.
Physical Review B , 13(4):1777, 1976.URL https://doi.org/10.1103/PhysRevB.13.1777 .51. Le Quang and Q.-C. He. The number and types of all possible rotational symmetries forflexoelectric tensors.
Proceedings of the Royal Society of London A: Mathematical, Physi-cal and Engineering Sciences , 467(2132):2369–2386, 2011. ISSN 1364-5021. URL https://doi.org/10.1098/rspa.2010.0521 .Xu Liang, Shuling Hu, and Shengping Shen. E ff ects of surface and flexoelectricity on a piezo-electric nanobeam. Smart materials and structures , 23(3):035020, 2014. URL https://doi.org/10.1088/0964-1726/23/3/035020 .EM Lifshitz and LD Landau. Statistical physics (course of theoretical physics, volume 5), 1951.Liwei Lin and Mu Chiao. Electro, thermal and elastic characterizations of suspended microbeams.
Microelectronics journal , 29(4-5):269–276, 1998. URL https://doi.org/10.1016/S0026-2692(97)00066-9 .M.E. Lines and A.M. Glass.
Principles and applications of ferroelectrics and related materials .Oxford University Press, 1979.Liping Liu. An energy formulation of continuum magneto-electro-elasticity with applications.
Journal of the Mechanics and Physics of Solids , 63:451–480, 2014. URL https://doi.org/10.1016/j.jmps.2013.08.001 .L.P. Liu and P. Sharma. Flexoelectricity and thermal fluctuations of lipid bilayer membranes:renormalization of flexoelectric, dielectric, and elastic properties.
Physical Review E , 87:032715, 2013.Wenhui Ma and L Eric Cross. Large flexoelectric polarization in ceramic lead magnesium nio-bate.
Applied Physics Letters , 79(26):4420–4422, 2001a. URL https://doi.org/10.1063/1.1426690 .Wenhui Ma and L Eric Cross. Observation of the flexoelectric e ff ect in relaxor Pb (Mg / Nb / )O ceramics. Applied Physics Letters , 78(19):2920–2921, 2001b. URL https://doi.org/10.1063/1.1356444 .Wenhui Ma and L Eric Cross. Flexoelectric polarization of barium strontium titanate in the para-electric state.
Applied Physics Letters , 81(18):3440–3442, 2002. URL https://doi.org/10.1063/1.1518559 .Wenhui Ma and L Eric Cross. Strain-gradient-induced electric polarization in lead zirconate ti-tanate ceramics.
Applied Physics Letters , 82(19):3293–3295, 2003. URL https://doi.org/10.1063/1.1570517 .Wenhui Ma and L Eric Cross. Flexoelectric e ff ect in ceramic lead zirconate titanate. AppliedPhysics Letters , 86(7):072905, 2005. URL https://doi.org/10.1063/1.1868078 .52enhui Ma and L Eric Cross. Flexoelectricity of barium titanate.
Applied Physics Letters , 88(23):232902, 2006. URL https://doi.org/10.1063/1.2211309 .M. S. Majdoub, P. Sharma, and T. C¸ a˘gin. Erratum: Enhanced size-dependent piezoelectricity andelasticity in nanostructures due to the flexoelectric e ff ect [phys. rev. b 77, 125424 (2008)]. Phys.Rev. B , 79:119904, 2009. doi: 10.1103 / PhysRevB.79.119904. URL https://doi.org/10.1103/PhysRevB.79.119904 .MS Majdoub, P Sharma, and T C¸ a˘gin. Enhanced size-dependent piezoelectricity and elasticity innanostructures due to the flexoelectric e ff ect. Physical Review B , 77(12):125424, 2008. URL https://doi.org/10.1103/PhysRevB.77.125424 .Sheng Mao and Prashant K. Purohit. Insights into flexoelectric solids from strain-gradient elas-ticity.
ASME Journal of Applied Mechanics , 81(8):1–10, 2014. URL https://doi.org/10.1115/1.4027451 .Sheng Mao, Prashant K Purohit, and Nikolaos Aravas. Mixed finite-element formulations in piezo-electricity and flexoelectricity.
Proceedings of the Royal Society A: Mathematical, Physical andEngineering Sciences , 472(2190):20150879, 2016. URL https://doi.org/10.1098/rspa.2015.0879 .R Maranganti, ND Sharma, and P Sharma. Electromechanical coupling in nonpiezoelectric materi-als due to nanoscale nonlocal size e ff ects: Green’s function solutions and embedded inclusions. Physical Review B , 74(1):014110, 2006. URL https://doi.org/10.1103/PhysRevB.74.014110 .J. P. Marcerou and J. Prost. The di ff erent aspects of flexoelectricity in nematics. Molecular Crystalsand Liquid Crystals , 58(3-4):259–284, 1980.M Marvan and A Havr´anek. Flexoelectric e ff ect in elastomers. In Relationships of PolymericStructure and Properties , pages 33–36. Springer, 1998. URL https://doi.org/10.1007/BFb0114342 .V.S. Mashkevich and K.B. Tolpygo. Electrical, optical and elastic properties of diamond typecrystals. 1.
Soviet Physics JETP-USSR , 5(3):435–439, 1957. URL .Andrew McBride, Denis Davydov, and Paul Steinmann. Modelling the flexoelectric e ff ect insolids: a micromorphic approach. arXiv preprint , 2019. URL https://arxiv.org/abs/1909.08695 .Robert B Meyer. Piezoelectric e ff ects in liquid crystals. Physical Review Letters , 22(18):918,1969. URL https://doi.org/10.1103/PhysRevLett.22.918 .53aymond David Mindlin. Micro-structure in linear elasticity.
Archive for Rational Mechanics andAnalysis , 16(1):51–78, 1964.Raymond David Mindlin. Polarization gradient in elastic dielectrics.
International Jour-nal of Solids and Structures , 4(6):637–642, 1968. URL https://doi.org/10.1016/0020-7683(68)90079-6 .R.D. Mindlin and N.N. Eshel. On first strain-gradient theories in linear elasticity.
InternationalJournal of Solids and Structures , 4(1):109–124, 1968. URL https://doi.org/10.1016/0020-7683(68)90036-X .P. Mohammadi, L.P. Liu, and P. Sharma. A theory of flexoelectric membranes and e ff ective prop-erties of heterogeneous membranes. Journal of Applied Mechanics , 81:011007, 2014.Anna N Morozovska, Eugene A Eliseev, Christian M Scherbakov, and Yulian M Vysochanskii.Influence of elastic strain gradient on the upper limit of flexocoupling strength, spatially modu-lated phases, and soft phonon dispersion in ferroics.
Physical Review B , 94(17):174112, 2016.URL https://doi.org/10.1103/PhysRevB.94.174112 .Anna N. Morozovska, Victoria V. Khist, Maya D. Glinchuk, Christian M. Scherbakov, Maxim V.Silibin, Dmitry V. Karpinsky, and Eugene A. Eliseev. Flexoelectricity induced spatially modu-lated phases in ferroics and liquid crystals.
Journal of Molecular Liquids , 267:550–559, 2018.SS Nanthakumar, Xiaoying Zhuang, Harold S Park, and Timon Rabczuk. Topology optimizationof flexoelectric structures.
Journal of the Mechanics and Physics of Solids , 105:217–234, 2017.URL https://doi.org/10.1016/j.jmps.2017.05.010 .BH Nguyen, X Zhuang, and Timon Rabczuk. Nurbs-based formulation for nonlinear electro-gradient elasticity in semiconductors.
Computer Methods in Applied Mechanics and Engineer-ing , 346:1074–1095, 2019. URL https://doi.org/10.1016/j.cma.2018.08.026 .Thanh D. Nguyen, Sheng Mao, Yao-Wen Yeh, Prashant K. Purohit, and Michael C. McAlpine.Nanoscale flexoelectricity.
Advanced Materials , 25(7):946–974, 2013. URL https://doi.org/10.1002/adma.201203852 .Ailish O’Halloran, Fergal O’malley, and Peter McHugh. A review on dielectric elastomer actua-tors, technology, applications, and challenges.
Journal of Applied Physics , 104(7):9, 2008. URL https://doi.org/10.1063/1.2981642 .M. A. Osipov and S. A. Pikin. Dipolar and quadrupolar ordering in ferroelectric iquid crystals.
J.Phys. II France , 5:1223–1240, 1995.Ronald E Pelrine, Roy D Kornbluh, and Jose P Joseph. Electrostriction of polymer dielectrics withcompliant electrodes as a means of actuation.
Sensors and Actuators A: Physical , 64(1):77–85,1998. URL https://doi.org/10.1016/S0924-4247(97)01657-9 .54G Petrov, RL Ramsey, and PNR Usherwood. Curvature-electric e ff ects in artificial and naturalmembranes studied using patch-clamp techniques. European Biophysics Journal , 17(1):13–17,1989. URL https://doi.org/10.1007/BF00257141 .Alexander G Petrov. Flexoelectric model for active transport. In
Physical and Chemical Basesof Biological Information Transfer , pages 111–125. Springer, 1975. URL https://doi.org/10.1007/978-1-4684-2181-1_9 .Alexander G. Petrov. Liquid crystal physics and the physics of living matter.
Molecular Crystalsand Liquid Crystals Science and Technology. Section A. Molecular Crystals and Liquid Crystals ,332(1):577–584, 1999.Alexander G Petrov. Flexoelectricity of model and living membranes.
Biochimica et Biophys-ica Acta (BBA)-Biomembranes , 1561(1):1–25, 2002. URL https://doi.org/10.1016/S0304-4157(01)00007-7 .L. Piegl and W. Tiller.
The NURBS Book . Monographs in Visual Communication. Springer BerlinHeidelberg, 2012. ISBN 9783642973857. doi: 10.1007 / https://books.google.es/books?id=58KqCAAAQBAJ .S.A. Pikin and V.L. indenbom. piezoe ff ects and ferroelectric phenomena in smectic liquid crystals. Ferroelectrics , 20:151–153, 1978.Roman Poya, Antonio J Gil, Rogelio Ortigosa, and Roberto Palma. On a family of numericalmodels for couple stress based flexoelectricity for continua and beams.
Journal of the Mechanicsand Physics of Solids , 125:613–652, 2019. URL https://doi.org/10.1016/j.jmps.2019.01.013 .Jacques Prost and JP Marcerou. On the microscopic interpretation of flexoelectricity.
Journal de Physique , 38(3):315–324, 1977. URL https://doi.org/10.1051/jphys:01977003803031500 .Ra ff aele Resta. Towards a bulk theory of flexoelectricity. Physical review letters , 105(12):127601,2010. URL https://doi.org/10.1103/PhysRevLett.105.127601 .A. D. Rey. Liquid crystal model of membrane flexoelectricity.
Physical Reviw E , 74(1):011710,2006.D.F. Rogers.
An Introduction to NURBS: With Historical Perspective . Morgan Kaufmann Seriesin Computer Graphics and Geometric Modeling. Morgan Kaufmann Publishers, 2001. ISBN9781558606692. URL https://doi.org/10.1016/B978-1-55860-669-2.X5000-3 .Samuel Rosset and Herbert R Shea. Small, fast, and tough: Shrinking down integrated elastomertransducers.
Applied Physics Reviews , 3(3):031105, 2016. URL https://doi.org/10.1063/1.4963164 . 55 Sahin and S Dost. A strain-gradients theory of elastic dielectrics with spatial dispersion.
Inter-national Journal of Engineering Science , 26(12):1231–1245, 1988. URL https://doi.org/10.1016/0020-7225(88)90043-2 .Andrea Schia ffi no, Cyrus E Dreyer, David Vanderbilt, and Massimiliano Stengel. Metric waveapproach to flexoelectricity within density functional perturbation theory. Physical Review B ,99(8):085107, 2019. URL https://doi.org/10.1103/PhysRevB.99.085107 .N.D. Sharma, C.M. Landis, and P. Sharma. Piezoelectric thin-film superlattices without us-ing piezoelectric materials.
Journal of Applied Physics , 108(2):1–25, 2010. doi: 10.1063 / http://dx.doi.org/10.1063/1.3443404 .Shengping Shen and Shuling Hu. A theory of flexoelectricity with surface e ff ect for elastic di-electrics. Journal of the Mechanics and Physics of Solids , 58(5):665 – 677, 2010. ISSN 0022-5096. URL https://doi.org/10.1016/j.jmps.2010.03.001 .Paul Steinmann and Duc Khoi Vu. Computational challenges in the simulation of nonlinear elec-troelasticity.
Computer Assisted Methods in Engineering and Science , 19(3):199–212, 2017.URL https://cames.ippt.pan.pl/index.php/cames/article/view/90 .Kai Sun. Toward molecular mechanoelectric sensors: Flexoelectric sensitivity of lipid bilayers tostructure, location, and orientation of bound amphiphilic ions.
The Journal of Physical Chem-istry B , 101(33):6327–6330, 1997.A. K. Tagantsev. Piezoelectricity and flexoelectricity in crystalline dielectrics.
Phys. Rev. B , 34:5883–5889, 1986. URL https://doi.org/10.1103/PhysRevB.34.5883 .AK Tagantsev. Theory of flexoelectric e ff ect in crystals. Zhurnal Eksperimental’noi i Teoretich-eskoi Fiziki , 88(6):2108–22, 1985. URL .Alexander K Tagantsev. Electric polarization in crystals and its response to thermal and elas-tic perturbations.
Phase Transitions: A Multinational Journal , 35(3-4):119–203, 1991. URL https://doi.org/10.1080/01411599108213201 .Tran Quoc Thai, Timon Rabczuk, and Xiaoying Zhuang. A large deformation isogeometric ap-proach for flexoelectricity and soft materials.
Computer Methods in Applied Mechanics andEngineering , 341:718–739, 2018. URL https://doi.org/10.1016/j.cma.2018.05.019 .Stephen P Timoshenko and James M Gere.
Theory of elastic stability . Courier Corporation, 2009.URL https://books.google.es/books?id=98B6JOW2HiUC .A Todorov, A Petrov, Michael O Brandt, and Janos H Fendler. Electrical and real-time strobo-scopic interferometric measurements of bilayer lipid membrane flexoelectricity.
Langmuir , 7(12):3127–3137, 1991. URL https://doi.org/10.1021/la00060a036 .56T Todorov, AG Petrov, and JH Fendler. First observation of the converse flexoelectric e ff ect inbilayer lipid membranes. The Journal of Physical Chemistry , 98(12):3076–3079, 1994. URL https://doi.org/10.1021/j100063a004 .KB Tolpygo. Long wavelength oscillations of diamond-type crystals including long range forces.
Soviet Physics-Solid State , 4(7):1297–1305, 1963.Richard Toupin. The elastic dielectric.
Journal of Rational Mechanics and Analysis , 5(6):849–915,1956.CL Trabi, CV Brown, AAT Smith, and NJ Mottram. Interferometric method for determining thesum of the flexoelectric coe ffi cients (e + e ) in an ionic nematic material. Applied PhysicsLetters , 92(22):223509, 2008. URL https://doi.org/10.1063/1.2938722 .DK Vu, P Steinmann, and G Possart. Numerical modelling of non-linear electroelasticity.
International Journal for Numerical Methods in Engineering , 70(6):685–704, 2007. URL https://doi.org/10.1002/nme.1902 .Bo Wang, Yijia Gu, Shujun Zhang, and Long-Qing Chen. Flexoelectricity in solids: Progress,challenges, and perspectives.
Progress in Materials Science , 2019. URL https://doi.org/10.1016/j.pmatsci.2019.05.003 .P.V. Yudin and A.K. Tagantsev. Fundamentals of flexoelectricity in solids.
Nanotechnology , 24(43):1–36, 2013. URL https://doi.org/10.1088/0957-4484/24/43/432001 .PV Yudin, R Ahluwalia, and AK Tagantsev. Upper bounds for flexoelectric coe ffi cients in ferro-electrics. Applied Physics Letters , 104(8):082913, 2014. URL https://doi.org/10.1063/1.4865208 .PV Yudin, R Ahluwalia, and AK Tagantsev. Erratum:“upper bounds for flexoelectric coe ffi cients inferroelectrics”[appl. phys. lett. 104, 082913 (2014)]. Applied Physics Letters , 106(18):189902,2015. URL https://doi.org/10.1063/1.4919883 .Julien Yvonnet and LP Liu. A numerical framework for modeling flexoelectricity and Maxwellstress in soft dielectrics at finite strains.
Computer Methods in Applied Mechanics and Engi-neering , 313:450–482, 2017. URL https://doi.org/10.1016/j.cma.2016.09.007 .Runzhi Zhang, Xu Liang, and Shengping Shen. A timoshenko dielectric beam model with flex-oelectric e ff ect. Meccanica , 51(5):1181–1188, 2016a. URL https://doi.org/10.1007/s11012-015-0290-1 .Shuwen Zhang, Minglong Xu, Xu Liang, and Shengping Shen. Shear flexoelectric coe ffi cient µ Journal of Applied Physics , 117(20):204102, 2015. URL https://doi.org/10.1063/1.4921444 . 57huwen Zhang, Minglong Xu, Guoliang Ma, Xu Liang, and Shengping Shen. Experimentalmethod research on transverse flexoelectric response of poly (vinylidene fluoride).
JapaneseJournal of Applied Physics , 55(7):071601, 2016b. URL https://doi.org/10.7567/JJAP.55.071601 .Yang Zhou, Jie Liu, Xinping Hu, Baojin Chu, Shutao Chen, and David Salem. Flexoelectric e ff ectin PVDF-based polymers. IEEE Transactions on Dielectrics and Electrical Insulation , 24(2):727–731, 2017. URL https://doi.org/10.1109/TDEI.2017.006273 .Xiaoying Zhuang, SS Nanthakumar, and Timon Rabczuk. A meshfree formulation for large de-formation analysis of flexoelectric structures accounting for the surface e ff ects. arXiv preprint ,2019. URL https://arxiv.org/abs/1911.06553 .Xiaoying Zhuang, Binh Huy Nguyen, Subbiah Srivilliputtur Nanthakumar, Thai Quoc Tran, NaifAlajlan, and Timon Rabczuk. Computational modeling of flexoelectricity—a review. Energies ,13(6):1326, 2020. URL https://doi.org/10.3390/en13061326 .P Zubko, G Catalan, A Buckley, PRL Welche, and JF Scott. Strain-gradient-induced polarizationin SrTiO single crystals. Physical Review Letters , 99(16):167601, 2007. URL https://doi.org/10.1103/PhysRevLett.99.167601 .Pavlo Zubko, Gustau Catalan, and Alexander K. Tagantsev. Flexoelectric e ff ect in solids. AnnualReview of Materials Research , 24(43):387–421, 2013. URL https://doi.org/10.1146/annurev-matsci-071312-121634https://doi.org/10.1146/annurev-matsci-071312-121634