We first revisit the mathematical modeling of the flexoelectric effect in the context of continuum mechanics at infinitesimal deformations. We establish and clarify the relation between the different formulations, point out theoretical and numerical issues related to the resulting boundary value problems, and present the natural extension to finite deformations. We then present a simple B-spline based computational technique to numerically solve the associated boundary value problems, which can be extended to handle unfitted meshes, hence allowing for arbitrarily-shaped geometries. Several numerical examples illustrate the flexoelectric effect in simple benchmark setups, as well as in new flexoelectric devices and metamaterials engineered for sensing or actuation.

Flexoelectricity generally refers to the two-way coupling between strain gradient and electric polarization (direct and inverse effects). Converse flexoelectricity couples polarization gradients with strain. All these effects have been demonstrated experimentally in cantilever thin films or nanobeams and truncated pyramids, under inhomogeneous mechanical fields (direct effect)1,2 and under homogeneous applied electric field (inverse effect),3–6 and in truncated pyramids and in PFM,7 under inhomogeneous applied electric field8 (converse effect). The setups mobilizing flexoelectricity, which necessarily involve non-uniform fields, demand accurate theoretical models for the quantitative evaluation of the effect and, in particular, for the correct interpretation of characterization experiments.9 

Phenomenological models for flexoelectricity in crystalline dielectrics were first proposed by Kogan10 after the early studies by Mashkevich and Tolpygo11 and Tolpygo.12 First comprehensive theoretical works by Tagantsev13,14 clarified the distinction between piezoelectricity and flexoelectricity. In the mechanics community, Mindlin15 formalized the converse flexoelectric effect in elastic dielectrics. A complete unified continuum framework, including strain gradient elasticity, both direct and converse flexoelectric couplings, and the polarization inertia effect was proposed later by Sahin and Dost.16 More recently, a simplified framework for isotropic dielectrics was proposed by Maranganti et al.17 

Nowadays, various continuum theories of flexoelectricity co-exist in the literature. Some of them are based on extensions to the electromechanical setting of gradient elasticity and other enriched continuum theories.18–24 Other authors consider the couplings with further physics, such as the flexoelectric effect in ferroelectrics,25,26 the coupling with magnetic fields27,28 or photovoltaics,29,30 and the contributions of surface effects.31 General variational principles for flexoelectric materials can be found in Liu,28 Shen and Hu,31 Hu and Shen.32 The reader is referred to Yudin and Tagantsev,33 Nguyen et al.,34 Zubko et al.,35 Krichen and Sharma,36 Wang et al.37 for comprehensive reviews of flexoelectricity in solids. Another focus of recent research is the modeling of flexoelectricity for soft materials (e.g., polymers and elastomers), which requires a finite deformation framework.22,28,38–41

All the aforementioned theories of flexoelectricity can be classified depending on the following considerations:

  • The choice of state variables describing the flexoelectric effect. For the mechanics, either the displacement gradient or its symmetrized form (i.e., strain) can be used, which give rise to type-I or type-II flexoelectricity, respectively. For the dielectrics, the most natural variable is the electric polarization. However, there exist many theories taking the electric field or the electric displacement instead.

  • The flexoelectric coupling considered explicitly, either the direct one, the converse, both, or the Lifshitz-invariant form. Section II A further elaborates on this topic.

  • The dielectric media surrounding the flexoelectric solid being included in the modeling or not. Since Maxwell equations hold in the whole R3, the electric potential must be solved as well outside the flexoelectric solid, in principle. However, if the dielectric permittivity of the flexoelectric material is large enough as compared to that of the surrounding media, zero-electric displacement around the flexoelectric solid is a reasonable assumption.

We present here an overview of various theoretical continuum formulations for flexoelectricity and the resulting computational models based on the discretization of the corresponding governing equations with smooth approximants. The paper is structured as follows. Section II presents the continuum modeling of type-II flexoelectricity models neglecting the surrounding media. It includes the discussion of flexoelectric couplings at infinitesimal deformations proposed in the literature, establishing the relationships between them, the different variational principles governing the behavior of flexoelectric solids, their associated boundary value problems and corresponding boundary conditions. The extension of the continuum modeling of flexoelectricity to finite deformations is stated and discussed. Section III presents a computational approach based on smooth B-spline basis functions, which can be combined with the immersed boundary method to allow for arbitrarily-shaped geometries. h-adaptivity is easily achieved by means of hierarchical refinement. The non-dimensionalization of the resulting equations is briefly commented. In Sec. IV, we use the computational approach discussed in Sec. III to solve the BVPs describing the electromechanical response of a flexoelectric cantilever beam in bending, under both mechanical and electrical input, based on the different formulations presented in Sec. II. A detailed comparison of the electromechanical output is given and several numerical issues are pointed out. In Sec. V, we explore the cantilever beam actuator in the regime of finite deformations, and study the interplay between flexoelectricity and electrostriction for different actuation levels. Then, in Sec. VI, we propose a new flexoelectric device and flexoelectricity-based piezoelectric metamaterials engineered to harness the flexoelectric effect toward sensing/actuation applications. Section VII concludes this paper.

In this section, we review the continuum modeling of the flexoelectric effect considering type-II flexoelectricity models at infinitesimal deformations and neglecting the surrounding media. In Sec. II A, we discuss the direct and converse flexoelectric couplings and their relation via null Lagrangians. Section II B presents the variational principles for flexoelectric solids, taking either the electric polarization or the electric potential as primal state variables, and explores their relation via a partial Legendre transform. Section II C formulates the boundary value problems of potential-based variational principles using the direct and the Lifshitz-invariant flexoelectricity models and introduces the weak enforcement of essential boundary conditions via Nitsche’s method. Non-standard boundary conditions, such as sensing electrode conditions, interface conditions and generalized periodicity conditions are briefly discussed as well. Finally, the modeling of flexoelectricity at finite deformations is introduced in Sec. II D.

Following the works by Maranganti et al.17 Yudin and Tagantsev,33 the most general (type-II) expression for the internal energy density describing the bulk static flexoelectric effect in centrosymmetric dielectrics can be written, under the assumption of infinitesimal deformations, in terms of the strain tensor ε, the electric polarization field P, and their corresponding spatial gradients in the following form:

ψ(0)(ε,ε,P,P)=12cijklεijεkl+12hijklmnεij,kεlm,n+12aklPkPl+12bijklPi,kPj,lflijk(1)εij,kPlflijk(2)εijPl,k,
(1)

where

  • c is the usual fourth-order elasticity tensor,

  • a is the usual second-order reciprocal dielectric susceptibility tensor,

  • h is the sixth-order strain gradient elasticity tensor, representing the purely non-local elastic effects,

  • b is the fourth-order polarization gradient tensor, representing the purely non-local effects of polarization,

  • f(1) is the direct flexocoupling tensor,

  • f(2) is the polarization gradient-strain coupling tensor introduced by Mindlin in his theory of polarization gradient.15 In a more modern context, it is also known as the converse flexocoupling tensor.

In non-centrosymmetric dielectrics, other relevant physics arise, such as piezoelectricity, the linear two-way coupling between strains and polarization. This coupling can be incorporated to Eq. (1) as an additional contribution to the internal energy density, ψpiezo=dijkεijPk, with the third-order piezoelectric tensor d. The interplay between flexoelectricity and piezoelectricity is very rich and worth studying (see, for instance, Abdollahi and Arias5 and Bhaskar et al.42). Here, we focus on the formulation of flexoelectricity in centrosymmetric dielectrics for the sake of simplicity. More complete phenomenological formulations for flexoelectricity have also been proposed and studied in the literature,43,44 which consider as well the linear couplings between strain and strain gradients, polarization and polarization gradients, and polarization gradients and strain gradients.

Assuming uniform material constants, the two latter terms in Eq. (1) can be rewritten as follows:

flijk(1)εij,kPlflijk(2)εijPl,k=flijkεij,kPlflijk(2)(εijPl)xk
(2a)
=flijkεijPl,kflijk(1)(εijPl)xk
(2b)
=12flijk(εij,kPlεijPl,k)12(flijk(1)+flijk(2))(εijPl)xk,
(2c)

with the (effective) flexocoupling tensor

flijk=flijk(1)flijk(2).
(3)

The first terms in (2a)(2c) are known, respectively, as the direct, converse, and Lifshitz invariant flexocouplings,45–47 and all of them represent both the direct and converse flexoelectric effects, as discussed later on. The second terms in (2a)(2c) are null-Lagrangians,48 in the sense that their bulk integral can be written as a surface integral by means of the divergence theorem, e.g., for (2c),

Ω12(flijk(1)+flijk(2))(εijPl)xkdΩ=Γ12(flijk(1)+flijk(2))εijPlnkdΓ.
(4)

Even though they affect boundary conditions, null Lagrangians are often omitted in the literature,33,46 yielding different internal energy densities to Eq. (1) as

ψ(Dir)(ε,ε,P,P)=12cijklεijεkl+12hijklmnεij,kεlm,n+12aklPkPl+12bijklPi,jPk,lflijkεij,kPl,
(5)
ψ(Con)(ε,ε,P,P)=12cijklεijεkl+12hijklmnεij,kεlm,n+12aklPkPl+12bijklPi,jPk,l+flijkεijPl,k,
(6)

and

ψ(Lif)(ε,ε,P,P)=12cijklεijεkl+12hijklmnεij,kεlm,n+12aklPkPl+12bijklPi,jPk,l12flijk(εij,kPlεijPl,k).
(7)

As derived later in Eq. (14), the bulk constitutive equations for the physical stress σ and physical electric field E are

El=ψPlxk(ψPl,k)=aklPkbijlkPi,jkflijkεij,k,
(8a)
σij=ψεijxk(ψεij,k)=cijklεklhijklmnεlm,nk+flijkPl,k,
(8b)

regardless of the choice of internal energy density, i.e., with ψ representing either the originalψ(0), directψ(Dir), converseψ(Con), or Lifshitz-invariantψ(Lif) energy densities of flexoelectric materials. Since formulations differing in null Lagrangians result in the same constitutive equations up to divergence-free fields, they also yield the same Euler–Lagrange equations. Null Lagrangians are thus viewed as modeling choices.49 In the absence of a clear physical interpretation of the higher-order boundary conditions, assuming natural boundary conditions to close the problem seems a reasonable choice. However, the expressions of the natural boundary conditions are actually different in each case, resulting in boundary value problems that are not equivalent, a fact that is commonly overlooked. See computational examples illustrating this fact in Secs. IV A and IV B.

The direct and converse flexoelectric effects are present in all the aforementioned formulations, as revealed in the constitutive equations for mechanics and dielectrics in Eq. (8). On the one hand, the direct flexoelectric effect in Eq. (8a) induces an electric field (or polarization) proportional to the strain gradient, whereas the converse effect in Eq. (8b) consists on a contribution to the mechanical stress proportional to the gradient of polarization. Note that both flexoelectric effects are governed by the same flexocoupling tensor f.

Frequently, the direct ψ(Dir) and converse ψ(Con) flexoelectricity models are simplified, so that only strain gradients or polarization gradients are present in the expression of the internal energy density. In the case of ψ(Dir), the polarization gradient tensor b is typically neglected, yielding

ψ(Dir)(ε,ε,P)=ψ(Dir)(ε,ε,P,P)|b=0.
(9)

In turn, the strain gradient elasticity tensor h can be neglected in ψ(Con), resulting in

ψ(Con)(ε,P,P)=ψ(Con)(ε,ε,P,P)|h=0.
(10)

Such simplifications are convenient (and hence popular) in order to facilitate the derivations of closed-form analytical solutions for simple flexoelectric devices, e.g., Euler–Bernoulli beams50–54 and Timoshenko beams55 and also to ease the implementation of numerical solution methods, e.g., the finite element (or related) methods. However, as reported in several references,15,18,28,56,57 either strain gradient or polarization gradient terms are required for a stable formulation.

The physics governing the behavior of dielectric materials can be modeled at a continuum level by different variational principles, depending on the state variables that are chosen to describe them. The usual choice for the mechanical state variable is the displacement field u, whose symmetric gradient yields the strain field, ε=12(u+(u)T). However, different options are frequently considered for the electrical state variable: either the electric polarization P, the electric displacement D, or the electric field E. Sections focus on two of them, P and E, which give rise to two different variational principles.

1. Minimization of the free energy Π[u, P]

The most natural choice is taking the polarization field P as the electrical state variable. This results in a variational formulation in terms of the free energy of the system such that, upon minimization over the admissible states, Euler–Lagrange equations and boundary conditions follow as necessary conditions.28 Such free energy takes the form17,28

Π[u,P]=Ω(ψ(u,P)+12ε0E2)dΩWext,
(11)

where Ω represents the flexoelectric material with internal energy density ψ (either ψ(0), ψ(Dir), ψ(Con), or ψ(Lif)), E is the electric field —i.e., work-conjugate to the polarization field P w.r.t. the internal energy density of the system, 12ε0E2 is the corresponding electrostatic energy density, with ε0 being the vacuum permittivity, and Wext is the external work.

The corresponding variational principle is stated as a constrained minimization problem of the form

(u,P)=argminuminPΠ[u,P],
(12a)

subject to

×E=0inΩ,
(12b)
D=qinΩ,
(12c)

where Eqs. (12b) and (12c) are the stationary Maxwell equations in dielectrics (i.e., Maxwell–Faraday’s and Gauss’s laws), D is the electric displacement—i.e., work-conjugate to the electric field E w.r.t. the total (internal + electrostatic) energy density of the system— and q is the external electric free charges per unit volume. Generally, Maxwell’s equations are stated over the whole space, including the materials and the ambient medium. Here, for the sake of simplicity, we restrict Eqs. (12b) and (12c) to the material domain Ω under the assumption that the flexoelectric material is surrounded by a medium with a much lower permittivity, e.g., air and thus sustains a vanishingly small electric field. Consideration of the full electrostatics including the ambient medium implies stating Maxwell’s equations in the whole space and setting electric continuity at the boundary of the material.

Such minimization leads to the following Euler–Lagrange equations:17,28

σ=fextinΩ,
(13a)
(ε0E+P)=qinΩ,
(13b)
×E=0inΩ,
(13c)

where fext represents the external body forces per unit volume, and

σ(u,P)=dψ(u,P)dε=[ψεψε](u,P),
(14a)
E(u,P)=dψ(u,P)dP=[ψPψP](u,P).
(14b)

The system of equations in (13) is composed by the differential Eqs. (13a) and (13b), corresponding to a (fourth-order) coupled elliptic problem and by Eq. (13c) which represents an additional constraint, requiring the irrotationality of the resulting electric field E. This constraint hinders the solution of the system in Eq. (13) to find u and P, either by means of analytical closed forms or numerical approximations.28 Hence, many authors prefer electric field-based models instead of polarization-based ones, since the electric field can be irrotational by construction, as explained next.

2. Optimization of the free enthalpy Π¯[u,ϕ]

Instead of taking the polarization field to be the electrical state variable, one (perhaps, less intuitive) alternative consists on considering the electric field E. Note that the Maxwell–Faraday’s law in Eq. (12b) implies the existence of an electric potential ϕ such that

E=ϕ.
(15)

Hence, by considering ϕ as the actual electrical state variable, Maxwell–Faraday’s law is automatically fulfilled, without the need of including it as a constraint.

A key ingredient required to proceed with this approach is the analogous expression to ψ(u,P) written in terms of the electric field E (and u) instead of P. Such an expression is known as the free enthalpy density—also called electrical Gibbs free energy density—of the system, which we denote by the symbol ψ¯(u,ϕ). Analogously to Sec. II A, one can define several enthalpy densities for flexoelectric materials, namely,

ψ¯(0)(ε,ε,E,E)=12c¯ijklεijεkl+12h¯ijklmnεij,kεlm,n12κklEkEl12MijklEi,jEk,lμlijk(1)εij,kElμlijk(2)εijEl,k,
(16)
ψ¯(Dir)(ε,ε,E,E)=12c¯ijklεijεkl+12h¯ijklmnεij,kεlm,n12κklEkEl12MijklEi,jEk,lμlijkεij,kEl,
(17)
ψ¯(Con)(ε,ε,E,E)=12c¯ijklεijεkl+12h¯ijklmnεij,kεlm,n12κklEkEl12MijklEi,jEk,l+μlijkεijEl,k,
(18)

and

ψ¯(Lif)(ε,ε,E,E)=12c¯ijklεijεkl+12h¯ijklmnεij,kεlm,n12κklEkEl12MijklEi,jEk,l12μlijk(εij,kElεijEl,k),
(19)

where c¯ is the elasticity tensor, κ is the electric permittivity tensor, h¯ is the strain gradient elasticity tensor, M is the gradient dielectricity tensor, and μ=μ(1)μ(2) is the flexoelectricity tensor. The direct and converse forms can be simplified as

ψ¯(Dir)(ε,ε,E)=ψ¯(Dir)(ε,ε,E,E)|M=0
(20)

and

ψ¯(Con)(ε,E,E)=ψ¯(Con)(ε,ε,E,E)|h¯=0,
(21)

where the gradient dielectricity tensor M or the gradient elasticity tensor h¯ is, respectively, vanished.

The associated free enthalpy of the system reads52,58,59

Π¯[u,ϕ]=Ωψ¯(u,ϕ)dΩWext,
(22)

where ψ¯ represents one of the aforementioned free enthalpy densities. The corresponding variational principle is stated as

(u,ϕ)=argminumaxϕΠ¯[u,ϕ],
(23)

which is now an unconstrained optimization (minmax) problem given that Eq. (12b) is automatically fulfilled and that Eq. (12c) follows directly from the solution of the variational principle (23) as one of the two associated Euler–Lagrange equations59 

σ=fextinΩ,
(24a)
D=qinΩ,
(24b)

where

σ(u,ϕ)=dψ¯(u,ϕ)dε=[ψ¯εψ¯ε](u,ϕ),
(25a)
D(u,ϕ)=dψ¯(u,ϕ)dE=[ψ¯E+ψ¯E](u,ϕ).
(25b)

Equations (24) and (25) describe a (fourth-order) coupled elliptic problem equivalent to that in Eqs. (13) and (14), but with no additional constraints, thus facilitating (in general) its solution via numerical or analytical method. If desired, the polarization field P is easily retrieved in terms of the obtained solution (u, ϕ) as

P(u,ϕ)=D(u,ϕ)+ε0ϕ.
(26)

3. Equivalence between different variational principles via partial Legendre transform

The variational models based on free energy Π(u,P) minimization and the ones based on free enthalpy Π¯(u,ϕ) optimization are related by means of a partial Legendre transform. That is, given one form, the other one is uniquely determined, and it can be directly obtained—under some assumptions that are next discussed—by a partial Legendre transform, thus revealing the relations between the material parameters of each form.

In particular, given an expression for the internal energy density ψ, the corresponding free enthalpy density ψ¯ is obtained by substracting the electrostatic energy density to the internal enthalpy density,28,60,61ψ¯Int, which is found via the following partial Legendre transform:

ψ¯(u,ϕ)=ψ¯Int(u,ϕ)12ε0E2=minP(ψ(u,P)PE)12ε0E2.
(27)

In the case of simplified direct internal energy density ψ(Dir)(ε,ε,P) in Eq. (9), i.e., without gradient polarization term (b=0), Eq. (27) reads

ψ¯(Dir)(u,ϕ)=minP(12cijklεijεkl+12hijklmnεij,kεlm,n+12amlPmPlflijkεij,kPlPlEl)12ε0E2.
(28)

The stationarity condition of the minimization results in

El=dψdPl=amlPmflijkεij,k,
(29)

which is explicitly invertible to

Pm=alm1(El+flijkεij,k)=alm1(El+ElFlexo),
(30)

where ElFlexo=flijkεij,k is well-known as the flexoelectric field. By inserting Eq. (30) into (28), the free enthalpy density is obtained as

ψ¯(Dir)(u,ϕ)=12cijklεijεkl12(aij1+ε0δij)EiEj+12(hijklmnalm1fmijkfblmn)εij,kεlm,n(alm1fmijk)Elεij,k.
(31)

The relation between the tensors c,a,f,h in the direct internal energy density (9) and their counterparts c¯,κ,μ,h¯ in the free enthalpy density (20) is revealed by Eq. (31) as follows:

c¯ijkl=cijkl,
(32a)
κij=aij1+ε0δij,
(32b)
μlijk=alm1fmijk,
(32c)
h¯ijklmn=hijklmnaab1faijkfblmn.
(32d)

Further, by assuming the standard expression for isotropic reciprocal dielectric susceptibility tensors a=(χeε0)1I=(εε0)1I, Eq. (32) simplifies to

c¯ijkl=cijkl,
(33a)
κij=εδij,
(33b)
μlijk=χeε0flijk,
(33c)
h¯ijklmn=hijklmnχeε0faijkfalmn,
(33d)

which yields the standard definition of κ for isotropic dielectrics, as a function of its electric permittivity ε, and reveals a well-known feature of flexoelectricity: its linear growth with the dielectric susceptibility χe. This is the reason why materials with high dielectric constant (e.g., ferroelectric perovskites) typically feature also large flexoelectric constants.34,35,59 Equation (33d) is also noticeable, since it shows that the strain gradient elasticity tensor in the free enthalpy is reduced due to the flexocoupling tensor.

It is easy to show that the same partial Legendre transform procedure, and hence Eq. (32), hold also for other density forms without gradient polarization/dielectricity terms (b,M=0), differing from (9) and (20) in null Lagrangians only.62 However, considering the gradient polarization/dielectricity terms (b,M0) breaks this equivalence, since the resulting stationarity condition of the minimization, analogous to Eq. (29), is rewritten as

El=dψdPl=amlPmflijkεij,kbijlkPi,jk,
(34)

which does not allow an explicit inversion to find P(E). An alternative consists on expressing the electric field E as a sum of local Elocal and nonlocal Enl fields as follows:

E=Elocal+Enl,
(35a)
Ellocal=amlPmflijkεij,k,
(35b)
Elnl=bijlkPi,jk
(35c)

and invert Eq. (35b), analogously to Eq. (30),

Pm=alm1(Ellocal+flijkεij,k).
(36)

With these definitions, both E(ε,Elocal) and P(ε,Elocal) can be expressed in terms of ε and Elocal, and introduced in Eq. (27) to yield

ψ¯(Dir)(u,Elocal)=12c¯ijklεijεkl+12h¯ijklmnεij,kεlm,n12κklEklocalEllocalμlijkεij,kEllocal12(aip1akq1bpjql)Ei,jlocalEk,llocal,
(37)

where higher order terms—i.e., terms depending on ε, Elocal or higher derivatives—have been neglected. Upon inspection of Eqs. (17) and (37), the following relation is revealed:

Mijkl=aip1akq1bpjql.
(38)

However, the equivalence between energy and enthalpy forms is not complete, since Eq. (37) is written in terms of the local part of the electric field only. In particular, this implies that Maxwell–Faraday’s law in Eq. (12b) is rewritten as

×E=×Elocal+×Enl=0,
(39)

hence ×Elocal0, preventing writing Elocal in terms of a scalar potential field ϕ. Therefore, the associated variational formulation remains constrained by Maxwell–Faraday’s law in Eq. (12b). This fact is sometimes overlooked in the literature, where there is no distinction made between Elocal and E, expressed in terms of ϕ via Eq. (15), regardless whether the enthalpy density features gradient dielectricity or not.

In this section, we consider the minmax variational principle on Π¯[u,ϕ] in Sec. II B2 and derive the corresponding boundary value problems for the Direct ψ¯(Dir)(ε,ε,E) and Lifshitz ψ¯(Lif)(ε,ε,E,E) free enthalpy densities in Eqs. (19) and (20).

Let Ω be a physical domain in RD,D={2,3}. The domain boundary, Ω, can be conformed by several smooth portions as Ω=fΩf [see Figs. 1(a) and 2(a)]. At each point xΩf, we define nf as the outward unit normal vector.

FIG. 1.

Sketch of the geometry of ΩR3. (a) Detail of Ω subdivided in smooth portions Ωi and Ωj, with their corresponding normal vectors ni and nj, (b) detail of Ωi, with the triplet {mi,si,ni} defined on Ωi, and (c) detail of Ωj, with the triplet {mj,sj,nj} defined on Ωj. The orientation of the tangent vector s is arbitrary and not relevant in this context. This figure was published in Codony et al.63 Copyright 2019 Elsevier.

FIG. 1.

Sketch of the geometry of ΩR3. (a) Detail of Ω subdivided in smooth portions Ωi and Ωj, with their corresponding normal vectors ni and nj, (b) detail of Ωi, with the triplet {mi,si,ni} defined on Ωi, and (c) detail of Ωj, with the triplet {mj,sj,nj} defined on Ωj. The orientation of the tangent vector s is arbitrary and not relevant in this context. This figure was published in Codony et al.63 Copyright 2019 Elsevier.

Close modal
FIG. 2.

Sketch of the geometry of ΩR2. (a) Ω and the boundary Ω subdivided in several smooth portions Ωf, e.g., Ωi and Ωj with their corresponding normal vectors ni and nj. The intersection between Ωi and Ωj is the corner Cij. (b) Detail of Ωi, with the pair {mi,ni} defined on Cij, and (c) detail of Ωj, with the pair {mj,nj} defined on Cij. This figure was published in Codony et al.63 Copyright 2019 Elsevier.

FIG. 2.

Sketch of the geometry of ΩR2. (a) Ω and the boundary Ω subdivided in several smooth portions Ωf, e.g., Ωi and Ωj with their corresponding normal vectors ni and nj. The intersection between Ωi and Ωj is the corner Cij. (b) Detail of Ωi, with the pair {mi,ni} defined on Cij, and (c) detail of Ωj, with the pair {mj,nj} defined on Cij. This figure was published in Codony et al.63 Copyright 2019 Elsevier.

Close modal

The boundary of the fth portion of Ω is a closed curve (D=3) or a pair of points (D=2) denoted as Ωf. The union C=fΩf represents the regions at which Ω is not smooth, and Cij=ΩiΩj, (ij), represents the part of C that is adjacent to Ωi and Ωj. In case D=3, at each point xΩf, we define mf as the unit co-normal vector pointing outward of Ωf, which is orthogonal to the normal vector nf and to the tangent vector sf of the curve Ωf [see Figs. 1(b) and 1(c)]. In case D=2, mf is defined on Ωf as the outward unit tangent vector to Ωf [see Figs. 2(b) and 2(c)].

The jump operator [[]] acting on a given quantity a is defined on Cij as the sum of that quantity evaluated at both sides of Cij, namely, [[a]]=ai+aj, where ak is the value of a from Ωk [see Fig. 1(a)].

1. Direct flexoelectricity

From Eqs. (20) and (22), it follows that

Π¯[u,ϕ]=Ω(12cijklεijεkl+12h¯ijklmnεij,kεlm,n12κklEkElμlijkεij,kEl)dΩΩWΩdΩΩWΩdΓCWCds,
(40)

where WΩ, WΩ, and WC represent the work densities exerted by external sources in the domain Ω, on its boundary Ω and at the non-smooth region C of its boundary, respectively. In this case, the admissible external sources of work are63 

WΩ(u,ϕ)=fextiuiqϕ,
(41a)
WΩ(u,ϕ)=tiui+rinuiwϕ,
(41b)
WC(u,ϕ)=jiui,
(41c)

where fext is the external body force per unit volume, q is the external free electric charge per unit volume, t and j are the forces per unit area (i.e., traction) and length, w is the surface charge density (i.e., electric charge per unit area), and r is the double traction (i.e., moment per unit area). The symbol n/n denotes the normal derivative.

Both the boundary of the domain Ω and its non-smooth region C are split into several disjoint Dirichlet and Neumann sets as

Ω=ΩuΩt=ΩvΩr=ΩϕΩw,
(42a)
C=CuCj,
(42b)

where Ωu, Cu, Ωv, and Ωϕ are the Dirichlet regions where either the displacement field, its normal derivative or the electric potential are prescribed,

u=u¯onΩuandCu,
(43a)
n(u)=v¯onΩv,
(43b)
ϕ=ϕ¯onΩϕ.
(43c)

The Neumann regions Ωt, Cj, Ωr, and Ωw correspond to prescribed values of their enthalpy conjugates, i.e., the traction t, the forces per unit length j, the double traction r or the surface charge density w. Whenever a given boundary patch Ωf belongs to Ωu, its boundaries Ωf must belong to Cu as well, namely, Cu=CΩu¯, and hence Cj=CCu.

Introducing the definitions (41) and (42) into Eq. (40) yields

Π¯[u,ϕ]=Ω(12cijklεijεkl+12h¯ijklmnεij,kεlm,n12κklEkElμlijkεij,kElfiextui+qϕ)dΩΩttiuidΓΩrrinuidΓ+ΩwwϕdΓCjjiuids.
(44)

According to Eq. (23), the equilibrium states (u,ϕ) correspond to the saddle points in the enthalpy potential (44) meeting the corresponding Dirichlet boundary conditions (43), implying that all the admissible variations with respect to the displacement and potential fields must vanish at (u,ϕ),

δΠ¯[u,ϕ]=δuΠ¯[u,ϕ]+δϕΠ¯[u,ϕ]=0,(δu,δϕ)U0P0.
(45)

The functional spaces of the admissible variations U0, P0 must meet homogeneous Dirichlet conditions, i.e.,

U0={δu[H2(Ω)]D|δu=0onΩuandCu,nδu=0onΩv},
(46a)
P0={δϕH1(Ω)|δϕ=0onΩϕ}.
(46b)

Introducing Eq. (44) into (45) yields the weak form of the problem.

Find (u,ϕ)UDPD such that (δu,δϕ)U0P0,

Ω(σ^ijδεij+σ~ijkδεij,kD^lδEl)dΩ=Ω(fextiδuiqδϕ)dΩ+ΩttiδuidΓ+ΩrrinδuidΓΩwwδϕdΓ+Cjjiδuids,
(47)

with functional spaces UD and PD fulfilling the Dirichlet boundary conditions in Eq. (43), that is,

UD={u[H2(Ω)]D|u=u¯onΩuandCu,nu=v¯onΩv},
(48a)
PD={ϕH1(Ω)|ϕ=ϕ¯onΩϕ}.
(48b)

The Cauchy stress σ^, the double stress σ~, and the electric displacement D^ fields in Eq. (47) are the conjugate quantities to the strain ε, the strain gradient ε, and the electric field E, respectively, as follows:

σ^ij=ψ¯(Dir)(ε,ε,E)εij=cijklεkl,
(49a)
σ~ijk=ψ¯(Dir)(ε,ε,E)εij,k=h¯ijklmnεlm,nμlijkEl,
(49b)
D^l=ψ¯(Dir)(ε,ε,E)El=κlmEm+μlijkεij,k.
(49c)
Remark 1:

In practice, in a finite element (FE) context, the functional spaces in Eq. (46) and (48) are approximated by means of a set of linear combinations of basis and test functions. However, UD and U0 are in general difficult to approximate since they require fulfilling second order Dirichlet conditions for the displacement field (i.e., prescribing its normal derivative). The typical approach to overcome this difficulty is considering Ωv=, which implies that only second order Neumann boundary conditions are allowed, i.e., Ω=Ωr.58,64–67 This choice is further justified by the unclear physical interpretation of second order Dirichlet boundary conditions.40 

The strong form of the problem is found by integrating the LHS of Eq. (47) by parts (twice for some terms) and invoking the divergence theorem on Ω, as follows:

Ω(σ^ijδεij+σ~ijkδεij,kD^lδEl)dΩ=Ω(σ^ij,jδuiσ~ijk,kδεijD^l,lδϕ)dΩ+Ω(σ^ijnjδui+D^lnlδϕ)dΓ+Ωσ~ijknkδui,jdΓ=Ω[(σ^ijσ~ijk,k),jδui+D^l,lδϕ]dΩ+Ω[(σ^ijσ~ijk,k)njδui+D^lnlδϕ]dΓ+Ωσ~ijknk(jS+njn)δuidΓ=Ω[(σ^ijσ~ijk,k),jδui+D^l,lδϕ]dΩ+Ω[(σ^ijσ~ijk,k)njδui+σ~ijknjnknδui+D^lnlδϕ]dΓ+Ωσ~ijknkjS(δui)dΓ,
(50)

where the gradient of the variations of the displacement field δui,j in the last integral on Ω is expressed in terms of its tangential jSδui and normal njnδui components, with S being the surface gradient operator. This separation is required since the tangential components of δu on Ω are dependent on δu, whereas the normal component is not. The last term in Eq. (50) is expressed in terms of δu via integration by parts and the surface divergence theorem on Ω,

Ωσ~ijknkjS(δui)dΓ=ΩjS(σ~ijknk)δuidΓ+ΩjS(σ~ijknkδui)dΓ=ΩjS(σ~ijknk)δuidΓ+ΩlS(nl)σ~ijknjnkδuidΓ+C[[σ~ijknkmk]]δuids.
(51)

Combining Eqs. (46), (47), (50), and (51), we obtain

Ω[[(σ^ijσ~ijk,k),j+fexti]δui+[D^l,lq]δϕ]dΩ=Ωt[(σ^ijσ~ijk,k+lS(nl)σ~ijknk)njjS(σ~ijknk)ti]δuidΓ+Ωr[σ~ijknjnkri]nδuidΓ+Ωw[D^lnl+w]δϕdΓ+Cj[[[σ~ijknkmk]]ji]δuids,
(52)

which reveals the strong form of the problem, composed by the Euler–Lagrange equations

{(σ^ijσ~ijk,k),j+fexti=0inΩ,D^l,lq=0inΩ,
(53)

Dirichlet boundary conditions (43) and Neumann boundary conditions

(σ^ijσ~ijk,k+lS(nl)σ~ijknk)njjS(σ~ijknk)=tionΩt,
(54a)
σ~ijknjnk=rionΩr,
(54b)
D^lnl=wonΩw,
(54c)
[[σ~ijkmjnk]]=jionCj.
(54d)

In view of Eq. (53), the definition of the physical stress σ arises naturally as

σij=σ^ijσ~ijk,k=cijklεklh¯ijklmnεlm,nk+μlijkEl,k,
(55)

whereas the physical electric displacement D is simply D^.

Remark 2:
The traction t can alternatively be written as63 
ti=(σ^ijσ~ijk,kkSσ~ikj)nj+σ~ijkN~jk,
(56)
where N~:=STr(S)nn is the second-order geometry tensor written in terms of the shape operator S of the surface Ω. The first term involves stress measures dotted with the normal vector n (a first-order measure of Ω, i.e., the orientation), whereas the second term involves the double stress dotted with the second-order geometry tensor (a second-order measure of Ω involving its curvature).

2. Lifshitz-invariant flexoelectricity

With this model, the admissible external sources of work are the ones corresponding to the direct flexoelectricity form in Eq. (41), plus the high order dielectric quantities (electric charge density per unit length) and r (double charge density, i.e., charge moment per unit area),62 analogous to j and r from mechanics

WΩ(u,ϕ)=fextiuiqϕ,
(57a)
WΩ(u,ϕ)=tiui+rinuiwϕrnϕ,
(57b)
WC(u,ϕ)=jiuiϕ.
(57c)

Accordingly, the boundary of the domain Ω is split into several disjoint Dirichlet and Neumann boundaries as

Ω=ΩuΩt=ΩvΩr=ΩϕΩw=ΩφΩr,
(58a)
C=CuCj=CϕC,
(58b)

where the high-order nature of the dielectrics leads to the definition of the Dirichlet Ωφ and Neumann Ωr boundaries corresponding to prescribed values for the normal derivative of the electric potential and its conjugate, i.e., the double charge density r. The edges C are also split into corresponding to the Dirichlet Cϕ and Neumann C edge sets, respectively, where either the electric potential ϕ or the electric charge density per unit length are prescribed. In this case, analogously to the mechanical partition of C, the sets Cϕ=CΩϕ¯ and C=CCϕ. The corresponding Dirichlet conditions of the problem are mathematically written as

u=u¯onΩuandCu,
(59a)
n(u)=v¯onΩv,
(59b)
ϕ=ϕ¯onΩϕandCϕ,
(59c)
n(ϕ)=φ¯onΩφ.
(59d)

With the definitions in (57) and (58) and the enthalpy density in (19), the enthalpy potential in (22) reads

Π¯[u,ϕ]=Ω(12cijklεijεkl+12h¯ijklmnεij,kεlm,n12κklEkEl12μlijk(εij,kElεijEl,k)12MijklEi,jEk,lfiextui+qϕ)dΩΩttiuidΓΩrrinuidΓCjjiuids+ΩwwϕdΓ+ΩrrnϕdΓ+Cϕds.
(60)

The variational principle in Eq. (23) implies that feasible equilibrium states uUD, ϕPD must meet Dirichlet boundary conditions (59), and their corresponding admissible variations δuU0, δϕP0 must vanish as stated in Eq. (45), where

UD={u[H2(Ω)]D|u=u¯onΩuandCu,nu=v¯onΩv},PD={ϕH2(Ω)|ϕ=0onΩϕandCϕ,
(61a)
nϕ=φ¯onΩφ},
(61b)

and

U0={δu[H2(Ω)]D|δu=0onΩuandCu,nδu=0onΩv},P0={δϕH2(Ω)|δϕ=0onΩϕandCϕ,
(62a)
nδϕ=0onΩφ}.
(62b)

Introducing Eq. (60) into (45) yields the weak form of the Lifshitz-invariant flexoelectricity problem

Find (u,ϕ)UDPD such that, (δu,δϕ)U0P0,

Ω(σ^ijδεij+σ~ijkδεij,kD^lδElD~lmδEl,m)dΩ=Ω(fextiδuiqδϕ)dΩ+ΩttiδuidΓ+ΩrrinδuidΓ+CjjiδuidsΩwwδϕdΓΩrrnδϕdΓCδϕds.
(63)

The Cauchy stress σ^, the double stress σ~, the electric displacement D^, and the double electric displacement D~ in Eq. (63) are the conjugate quantities to the strain ε, the strain gradient ε, the electric field E, and the electric field gradient E, respectively, as follows:

σ^ij=ψ¯(Lif)εij=cijklεkl+12μlijkEl,k,
(64a)
σ~ijk=ψ¯(Lif)εij,k=h¯ijklmnεlm,n12μlijkEl,
(64b)
D^l=ψ¯(Lif)El=κlmEm+12μlijkεij,k,
(64c)
D~lk=ψ¯(Lif)El,k=MijlkEi,j12μlijkεij.
(64d)

The strong form of the problem is found analogously to the Direct flexoelectricity problem and results in the following Euler–Lagrange equations:

{(σ^ijσ~ijk,k),j+fexti=0inΩ,(D^lD~lk,k),lq=0inΩ,
(65)

complemented with Dirichlet conditions (59) and Neumann conditions

(σ^ijσ~ijk,k+lS(nl)σ~ijknk)njjS(σ~ijknk)=tionΩt,
(66a)
[[σ~ijkmjnk]]=jionCj,
(66b)
σ~ijknjnk=rionΩr,(D^lD~lk,k+iS(ni)D~lknk)nl+lS(D~lknk)=wonΩw,
(66c)
[[D~jkmjnk]]=onC,
(66d)
D~jknjnk=ronΩr.
(66e)

The physical stress σ and physical electric displacement D arise from Eq. (65) as

σij=σ^ijσ~ijk,k=cijklεklh¯ijklmnεlm,nk+μlijkEl,k,
(67a)
Dl=D^lD~lk,k=κlmEmMijlkEi,jk+μlijkεij,k.
(67b)

3. Weak enforcement of Dirichlet boundary conditions via Nitsche’s method

So far, the functional spaces of the state variables (48) or (61) are chosen such that Dirichlet boundary conditions are automatically fulfilled. In computational flexoelectricity, this implies having an interpolant space on the boundaries Ω where Dirichlet conditions are to be enforced, e.g., by means of body-fitted meshes in a FE approach. However, it is typically difficult to achieve such functional spaces while fulfilling the C1 continuity requirement.59 

An alternative to enforce essential boundary conditions without constraining the functional spaces consists on incorporating them into the enthalpy functional, in a way that equilibrium states (u,ϕ) satisfying the corresponding variational principle necessarily fulfill (in a weak sense) the Dirichlet boundary conditions. We show here how to do so via Nitsche’s method68 due to its simple form and convenient numerical properties (i.e., self-consistency, symmetry, optimal error convergence rates, and preservation of the number of degrees of freedom69) as compared to other alternatives such as the Lagrange multipliers or penalty methods.

Following Nitsche’s approach, an additional term Π¯Nitsche[u,ϕ] acting on the Dirichlet boundaries is introduced in the enthalpy functional in Eq. (44) or Eq. (60) to weakly enforce Dirichlet boundary conditions Eq. (43) or Eq. (59), respectively. In the case of Lifshitz-invariant formulation, it has the following form:

Π¯Nitsche[u,ϕ]=Ωu(12βu(uiu¯i)2(uiu¯i)ti(u,ϕ))dΓ+Ωv(12βv(nuiv¯i)2(nuiv¯i)ri(u,ϕ))dΓ+Cu(12βCu(uiu¯i)2(uiu¯i)ji(u,ϕ))dsΩϕ(12βϕ(ϕϕ¯)2(ϕϕ¯)w(u,ϕ))dΓΩφ(12βφ(nϕφ¯)2(nϕφ¯)r(u,ϕ))dΓCu(12βCϕ(ϕϕ¯)2(ϕϕ¯)(u,ϕ))ds,
(68)

with the numerical parameters βu,βv,βCu,βϕ,βφ,βCϕR+ and Neumann terms explicitly depending on the state variables as defined in Eq. (66). For the direct formulation, the last two lines in Eq. (68) must be omitted according to the Dirichlet conditions (43). The penalty terms inserted in each boundary integral are quadratic in the Dirichlet boundary conditions, and its only purpose is to ensure equilibrium states (u,ϕ) being, respectively, actual minima and maxima of the energy functional with respect to u and ϕ.63 

Remark 3:
Nitsche’s formulations are self-consistent for any value of the penalty parameters βu,βv,βCu,βϕ,βφ,βCϕ. However, in the discrete space of numerical approximation of the state variables, they must be large enough to ensure stability, i.e., maintain the minmax nature of the corresponding variational principle. Nevertheless, arbitrarily large values are not suitable since the conditioning of the linear system is deteriorated. The analytical derivation of lower bounds of the penalty parameters can be found in Codony et al.63 for the discrete case using direct flexoelectricity, but moderate values of the penalty parameters are typically enough to ensure convergence and enforce boundary conditions properly.63,69,70 Thus, the explicit computation of stability lower bounds can be avoided by writing the penalty parameters in terms of a mesh-independent dimensionless parameter ζR+ as follows:
βu=Yhζ,βv=(mech2+μ2)βu,βCu=βvh,βϕ=εhζ,βφ=(elec2+μ2)βϕ,βCϕ=βφh,
(69)
where h denotes the characteristic length of the discretization (typically, the mesh size), Y is the Young modulus, ε is the electric permittivity, mech is the mechanical length scale, elec is the dielectric length scale, μμ/Yε is the flexoelectric length scale, and μ is the flexoelectric tensor (see  Appendix A for further details on material parameters). A suitable value of ζ can be determined empirically, e.g., in smooth B-spline approximations ζ=10100 typically provides stable results.63 

The variational principle associated with the enthalpy functional including Nitsche’s term leads to the same Euler–Lagrange equations and definition of Neumann terms.63 However, the functional spaces for the state variables and their variations are now unconstrained, at the cost of requiring more regularity.63 For the direct flexoelectricity form,

UD=U0=[H2(Ω)]DwithL2integrablethirdderivativesontheboundaryΩu,
(70a)
PD=P0=H1(Ω),
(70b)

where UD and U0 account for the integrals involving t in Eq. (54a), requiring the computation of third order derivatives. In practice, UD=U0=[H3(Ω)]D are typically considered.63 Analogously, the Lifshitz-invariant form requires

UD=U0=[H2(Ω)]DwithL2integrablethirdderivativesontheboundaryΩu,
(71a)
PD=P0=H2(Ω)withL2integrablethirdderivativesontheboundaryΩϕ.
(71b)
Remark 4:

Many authors in computational flexoelectricity neglect the edge Dirichlet conditions (43a) on Cu9,57,58,71 [or analogously (59a) on Cϕ]. In the cases where essential boundary conditions are enforced strongly, as in conforming FE or meshless discretizations, this fact has no practical relevance since the strong imposition on the surface Ωu (Ωϕ) automatically implies the strong imposition on the adjacent edges in Cu (Cϕ) as well. However, it is important to underline that in frameworks where boundary conditions are enforced weakly, dismissing edge conditions is equivalent to considering homogeneous Neumann edge conditions, which is wrong on Dirichlet edges, and can completely spoil the quality of the numerical results.63 

4. Beyond standard boundary conditions

So far, we have mentioned that the standard conditions along the boundary of the flexoelectric domain, arising naturally from the derivation of the boundary value problems, are either Dirichlet or Neumann. However, it is worth mentioning other types of boundary conditions that are also frequent in electromechanical boundary value problems:

a. Sensing electrode condition

The electric potential along electrodes placed on the surface of the devices is uniform, since they are made of conducting material. Electrodes in actuation mode correspond to a standard Dirichlet boundary condition for ϕ [see Eq. (43c) or (59c)], where the prescribed voltage is constant.

In the case of electrodes in sensing mode, the electric potential is uniform but a priori unknown. The corresponding boundary condition at each electrode ΩΦi and its edge ΩΦi is mathematically written as

ϕ=ΦionΩΦiandΩΦi,i=1,,NΦ,
(72)

where Φ1,,ΦNΦR are additional state variables and NΦ denotes the number of sensing electrodes in the system. Since each Φi is not prescribed on ΩΦi nor ΩΦi, their energy conjugates—i.e., the total surface charge density on the electrode63 and on their edges—must vanish, that is,

ΩΦiwdΩ=0onΩΦi,i=1,,NΦ,
(73a)
ΩΦids=0onΩΦi,i=1,,NΦ.
(73b)

In the direct flexoelectricity model, high-order electrostatics are not present, hence Eqs. (72) and (73) must hold on ΩΦi only.

In a discrete boundary value problem with body-fitted discretization, Eq. (72) can be enforced strongly by assigning the same degree of freedom for ϕ to all the nodes on the electrode, and Eq. (73) is then automatically fulfilled. Otherwise, they must be enforced weakly, e.g., via Lagrange multipliers72 or a modified version of Nitsche’s method.63 

b. Interface condition

Composite materials, i.e., domains composed by multiple materials with different properties, can be used to achieve a specific functionality of an electromechanical device. The paradigmatic example consists on cantilevers bimorphs for piezoelectric transduction,5,73,74 but one can also think of flexoelectricity-based devices with multimaterial stacks46,75,76 or structured materials with geometrically polarized cavities2,77–79 and even conducting/insulating inclusions.80,81

In the context of discrete flexoelectricity boundary value problems, high-order interface conditions must be considered.81 Such conditions correspond to high-order continuity of the state variables and high-order equilibrium, (i) across the interface I between two adjacents subdomains, as well as (ii) at the sharp regions S of the interfaces, shared by multiple—possibly more than two—subdomains (see Fig. 3).

FIG. 3.

(a) Composite domain Ω=Ω1Ω2Ω3. (b) Detail of the interface I between Ω1 and Ω2. (c) Detail of the point S shared by all three subdomains.

FIG. 3.

(a) Composite domain Ω=Ω1Ω2Ω3. (b) Detail of the interface I between Ω1 and Ω2. (c) Detail of the point S shared by all three subdomains.

Close modal

In the Lifshitz-invariant model, high-order continuity conditions read

uL=uRonI,
(74a)
ui=ujonS,
(74b)
uL/nL=uR/nLonI,
(74c)
ϕL=ϕRonI,
(74d)
ϕi=ϕjonS,
(74e)
ϕL/nL=ϕR/nLonI,
(74f)

and high-order equilibrium conditions read

tL+tR=0onI,
(75a)
kjk=0onS,
(75b)
rLrR=0onI,
(75c)
wL+wR=0onI,
(75d)
kk=0onS,
(75e)
rLrR=0onI,
(75f)

where L and R represent either sides across the interface I, ij represents all the pairwise combinations of subdomains adjacent to S, and k sums over all the subdomains adjacent to S. In the direct flexoelectricity model, high order electromechanics are not present, i.e., Eqs. (74e), (74f), (75e), and (75f) must not be considered.

The conditions in Eqs. (74) and (75) are implicitly fulfilled in C1-conforming body-fitted discretizations. Otherwise, they must be explicitly enforced, e.g., weakly via Nitsche’s method.81 

c. Generalized periodicity condition

Architected materials exhibit an intrinsic microstructure that is replicated periodically in space. In order to numerically compute their overall bulk response, free of finite edge effects, undergoing macroscopic displacement (sensor) or electric field (actuator) conditions, an efficient strategy consists on exploiting the lattice periodicity by studying a single representative volume element (RVE) or unit cell under (high-order) generalized-periodic conditions.82,83 These set of conditions are a generalization of classical (high-order) periodicity conditions to account for non-periodic solution fields with periodic gradients.

Consider an infinite lattice with translational periodicity (L1e1,L2e2,L3e3) in R3. Without loss of generality, the flexoelectricity boundary value problem can be reduced to a single RVE defined as Ω(0,L1)×(0,L2)×(0,L3) with the proper generalized periodic conditions (i) between a given periodic boundary A of the RVE (i.e., x1=0, x2=0, or x3=0) and its corresponding translational image A (i.e., x1=L1, x2=L2, or x3=L3), as well as (ii) between the intersections of A and A with the actual boundary of the infinite lattice, denoted as S and S, respectively (see Fig. 4).

FIG. 4.

(a) Architechted material with translational periodicity. A unit cell (or RVE) is highlighted. (b) Periodic boundaries of the RVE. (c) Intersections of the periodic boundaries of the RVE with the actual boundaries of the architected material.

FIG. 4.

(a) Architechted material with translational periodicity. A unit cell (or RVE) is highlighted. (b) Periodic boundaries of the RVE. (c) Intersections of the periodic boundaries of the RVE with the actual boundaries of the architected material.

Close modal

In the Lifshitz-invariant model, they read as follows:

uA=uA[[u]]AonAandS,
(76a)
uA/nA=uA/nAonA,
(76b)
ϕA=ϕA[[ϕ]]AonAandS,
(76c)
ϕA/nA=ϕA/nAonA,
(76d)
tA+tA=0onA,
(76e)
jS+jS=0onS,
(76f)
rArA=0onA,
(76g)
wA+wA=0onA,
(76h)
S+S=0onS,
(76i)
rArA=0onA,
(76j)

where [[]]A are either known (enforced) or unknown (resolved) constant quantities denoting the difference of the state variable within the brackets () across either sides of the RVE along the direction orthogonal to the plane A. For the direct flexoelectricity model, Eqs. (76d), (76i) and (76j) are disregarded and Eq. (76c) is enforced on A only.

In the context of flexoelectricity, generalized periodicity conditions allow the systematic study and design of flexoelectric metamaterials exhibiting an apparent piezoelectric behavior.83 An elegant way of enforcing high-order generalized-periodic conditions in a discrete boundary value problem consists on simply constructing a high-order generalized-periodic approximation space.82 Alternatively, one can resort to enforcing generalized periodicity conditions (76) weakly.81 

The interest on computational flexoelectricity in soft materials has increased in recent years22,38–41,64,84 due to different reasons. On the one hand, a large flexoelectric response is expected, since flexoelectric coefficients of polymers are at least the same order of magnitude as those of hard crystalline materials,85–87 but, in turn, they are much more compliant. On the other hand, electromechanical actuation of polymers by flexoelectricity overcomes the current limitations of traditional actuation (based on the Maxwell stress effect, i.e., an electromechanical coupling inducing a contribution in the stress field quadratic to the electric field): (i) one-way coupling, (ii) dielectric breakdown due to too large electric fields, and (iii) irreversibility of the deformation sign.36,88–90 We present here a brief introduction of computational flexoelectricity in the regime of finite deformations and review its most relevant features.

Consider a deformable dielectric body described by Ω0 in the reference (or undeformed) configuration and by Ω in the current (or deformed) configuration. The deformation map χ:Ω0Ω maps every material point XΩ0 to the spatial point x=χ(X)Ω. Relevant deformation measures in this context are the deformation gradient F=0χ, the Jacobian determinant J=det(F), and the right Cauchy–Green deformation tensor C=FTF. In the reference configuration, the strain is represented by the Green–Lagrangian strain tensor E=(CI)/2, and the electric field in by E=0Φ.

Analogously to the derivations for infinitesimal deformations, the Lagrangian description of the electromechanical enthalpy density per unit undeformed volume corresponding to a dielectric material with simplified direct flexoelectric coupling can be derived via a partial Legendre transform of the energy density. The detailed derivation can be found in Codony et al.,41 resulting in

Ψ¯Enth(E,0E,E)=ΨMech(E,0E)12Eκ¯EEμ¯0E,
(77)

where ΨMech is an hyperelastic potential accounting for the local and non-local mechanics, and the tensors

κ¯(E)=JC1ε,
(78a)
μ¯(E)=JC1μ
(78b)

are the effective dielectricity and flexoelectricity tensors,41 respectively, which explicitly depend on the Green–Lagrangian strain E.

The enthalpy functional governing the physics of a flexoelectric body is written, analogously to Eq. (44), as41 

Π¯[χ,Φ]=Ω0(Ψ¯Enth(E,0E,E)Fextiχi+QΦ)dΩ0Ω0TTiχidΓ0Ω0RRi0NχidΓ0C0JJiχids0+Ω0WWΦdΓ0,
(79)

where Fext is the external body force per unit reference volume, considered constant in the deformed frame (dead load), Q is the external free electric charge per unit reference volume, T is the external traction in deformed frame per unit reference area, R is the external double traction in deformed frame per unit reference area, J is the external force per unit reference length, and W is the external charge density per unit reference area.

The equilibrium states {χ,Φ} are located at the saddle points of Π¯ satisfying

{χ,Φ}=argminχXmaxΦPΠ¯[χ,Φ],
(80)

where X and P are the functional spaces for χ and Φ with sufficient regularity fulfilling Dirichlet boundary conditions. Vanishing the first variations of Π¯[χ,Φ] leads to the weak form of the boundary value problem, whose corresponding Euler–Lagrange equations read as follows:41 

(FiISIJ),J+Fexti=0,
(81a)
DK,KQ=0.
(81b)

The physical second Piola–Kirchhoff stress S and the electric displacement D are defined in terms of χ and Φ as

SIJ(χ,Φ)=SMechIJ(χ,Φ)+SMaxwellIJ(χ,Φ),
(82a)
DL(χ,Φ)=JCKL1(εEK+μKIJMEIJ,M),
(82b)

where

SMechIJ=ΨMechEIJXK(ΨMechEIJ,K)+JCLM1EM,KμLIJK
(83)

and D correspond to the extensions of Eqs. (49c) and (55) to finite deformations, respectively, and the additional term

SMaxwellIJ=J2(CMI1CLJ1+CLI1CMJ1CML1CIJ1)EMεEL
(84)

is the total second Piola–Maxwell stress tensor,41 arising as a direct consequence of the strain-dependence of the effective dielectricity tensor in Eq. (78a). It is also known as the electrostrictive effect, a stress contribution intrinsically present in the model at finite deformations (as long as the material is dielectric) that is proportional to the square of the electric field.

In summary, the regime of finite deformations presents a richer electromechanical behavior as compared to the one in the limit of infinitesimal deformations. On the one hand, one can study the interplay between the electrostrictive and the flexoelectric effects in different scenarios, as illustrated in Sec. V with the actuation of a soft flexoelectric cantilever. On the other hand, the nonlinear definition of the strain field implies geometrical nonlinearity, allowing the occurrence of mechanical instabilities such as buckling, wrinkling or creasing, which can be harnessed toward the design of tunable flexoelectricity-based devices.41 

The equations of flexoelectricity can only be solved analytically in very simple settings, such as simplified Euler–Bernoulli beam,51,52,91 Timoshenko beam,55 or Cosserat rod41 models. Otherwise, it is necessary to resort to computational flexoelectricity.

The major challenge is to handle the C1 continuity of the state variables required by the fourth-order PDE system. To address this, several numerical alternatives have been proposed, i.e., mesh-free approximations,5,9,58,64,92 isogeometric analysis,39,41,65,84,93 the C1 Argyris triangular element approximation,38 the immersed boundary B-spline method,63 or the C0-interior penalty method,72 among others. Another family of numerical methods are those circumventing the C1 continuity requirement by introducing additional variables, such as mixed formulations,66,67,94 or those based on micromorphic theories of continua.22,40 We refer to Zhuang et al.59 for a comprehensive review on computational approaches to solve flexoelectricity boundary value problems.

In this tutorial, in Sec. III A, we focus first on open B-spline approximations in body-fitted meshes (i.e., a particular case within isogeometric analysis). It is a simple method to successfully resolve the flexoelectricity equations in a smooth approximation space, albeit only rectangular geometries (or cuboidal ones in 3D) are allowed. It is, therefore, a suitable method for the computational modeling of flexoelectricity in beams. In Sec. III B, in order to circumvent this limitation and allow arbitrarily-shaped geometries, we combine this approximation with the immersed boundary method,95,96 involving non-conforming discretizations of the geometry and non-interpolant discretizations of the state variables. Dirichlet boundary conditions must be weakly enforced, and some further numerical treatments in cut cells must be accounted for, such as defining a numerical integration for cut elements and alleviating ill-conditioning produced by degrees of freedom with small support within the domain. Section III C presents the extension to hierarchical B-spline basis, allowing h-refinement of the computational mesh, and the adimensionalization of the resulting equations is shown in Sec. III D.

Using these approaches, as compared to the aforementioned ones, yields multiple benefits, i.e., (i) trivial mesh generation, (ii) easy handling of arbitrarily-shaped geometries, (iii) optimal convergence rates for the approximation error, that can be systematically increased using higher-order B-spline bases, (iv) cheap and robust evaluation of basis functions and derivatives, (v) relatively small number of degrees of freedom for a given spatial resolution, and (vi) good scalability and relatively low matrix fill-in due to a compact and structured set of basis functions.

B-spline functions97–99 are smooth, positive-valued piece-wise polynomials with compact support. Being p the polynomial degree, they are by construction Cp1-continuous throughout the domain and form a partition of unity. Since they are polynomials in nature, they are explicitly evaluated and exactly integrated by rich enough numerical quadratures. However, they do not satisfy in general the Kronecker delta property for p>1.

A univariate B-spline basis of degree p consisting of n basis functions on m=n+p cells is defined on the parametric space ξΞ=[k0,km] in terms of the knot vector k=[k0,k1,k2,,km], where ki<ki+1. The ith function of this basis is defined recursively as97 

Bi0(ξ)={1ki1ξ<ki,0otherwise,Bip(ξ)=ξki1ki+p1ki1Bip1(ξ)+ki+pξki+pkiBi+1p1(ξ),i=1,,n.
(85)
Remark 5:

In Eq. (85) and throughout the rest of the paper, we follow the convention that the first basis function corresponds to index i=1 and spans from ξ=k0 to ξ=kp+1. Note that other references might consider 0-based indexing of the basis functions or 1-based indexing for the knot vector instead.

Different kinds of B-spline bases are obtained depending on the choice of the knot vector. For simplicity, we consider equispaced knots ki=i, yielding a uniform B-spline basis where the ith B-spline function can be expressed as a translation of the first one as Bip(ξ)=B1p(ξi+1) (see Fig. 5).

FIG. 5.

First univariate uniform B-spline basis function B1p(ξ) of degree p: (a) linear (p=1), (b) quadratic (p=2), (c) cubic (p=3), (d) quartic (p=4).

FIG. 5.

First univariate uniform B-spline basis function B1p(ξ) of degree p: (a) linear (p=1), (b) quadratic (p=2), (c) cubic (p=3), (d) quartic (p=4).

Close modal

In order to enforce boundary conditions in a strong way, B-spline bases must satisfy the Kronecker delta property at the boundary, corresponding to the knots k0 and km in parametric space. To do so, the basis is modified by knot multiplicity, i.e., k=[k0,,k0,k1,k2,,km,,km], where k0 and km are repeated p+1 times. In doing so, the continuity of the B-spline basis at these knots is decreased to C1 (i.e., discontinuous), yielding a boundary-interpolant (or open) basis suitable to enforce essential boundary conditions strongly (see Fig. 6).

FIG. 6.

Univariate open uniform cubic (p=3) B-spline basis with knot vector k=[0,0,0,0,1,2,3,4,5,6,7,8,8,8,8].

FIG. 6.

Univariate open uniform cubic (p=3) B-spline basis with knot vector k=[0,0,0,0,1,2,3,4,5,6,7,8,8,8,8].

Close modal

In a D-dimensional space, the ith B-spline function Bip(ξ) of a D-variate B-spline basis (where i is the D-variate index [i1,,iD]) is defined as the tensor product of D univariate B-spline functions as

Bip(ξ)=B[i1,,iD]p([ξ1,,ξD])=d=1DBidp(ξd),id=1,,nd,
(86)

which is defined on the D-dimensional parametric space ξΞ=[k1,0,k1,m1][kd,0,kd,md].

In order to use B-spline basis (defined in the parametric space Ξ) to approximate the state variables u and ϕ (defined in the physical space Ω), let us define the geometrical map

φ:ΞΩξx=φ(ξ),
(87)

which maps a given point ξΞ in the parametric space to a given point xΩ in the physical space. The basis functions N(x)Ω are defined as N=[Bp°φ1], in such a way that

[u(x)]d[uh(x)]d=Ni(x)auid=Bip(ξ)auid,
(88a)
ϕ(x)ϕh(x)=Ni(x)aϕi=Bip(ξ)aϕi,
(88b)

for d=1,,D, where {au,aϕ} are the degrees of freedom (DOF) of the approximations uh(x) and ϕh(x). Since we restrict ourselves to the case of Ω being a rectangle in case D=2, or cuboid in case D=3, i.e., Ω=x0+[0,L1][0,LD], the map φ simply corresponds to an affine transformation composed by a non-uniform scaling and a translation,

x=φ(ξ)=x0+Jξ,J=[h1000000hD],
(89)

where hd is the size of each cell in physical space along the dimension d, determined as hd=Ld/md.

This simple expression, with a constant and diagonal Jacobian matrix J, facilitates the computation of (high-order) gradients of the basis functions N(x)Ω, since the tensor product structure of the basis is preserved, e.g.,

xjNi(x)=hj1(ξjBijp(ξj))djDBidp(ξd),
(90a)
xj[xjNi](x)=hj2(22ξjBijp(ξj))djDBidp(ξd),xk[xjNi](x)=hj1hk1(ξjBijp(ξj))(ξkBikp(ξk))dj,dkDBidp(ξd),
(90b)
forjk,
(90c)

and so on.

Using the multivariate B-spline basis functions in Eq. (86) and the geometrical map in Eq. (89), the state variables in the boundary value problem are approximated by means of Eq. (88). Volume, surface, and line integrals are numerically computed by standard Gaussian quadrature rules, and boundary conditions are enforced strongly. In univariate open bases, the DOF corresponding to the first and last basis functions can be directly prescribed with the value of the boundary condition. In multivariate open bases, however, the values of the control variables on the boundary are computed by means of the L2 projection of the boundary condition onto the space spanned by the corresponding B-spline basis functions.

The discretization and numerical approximation of the weak form yield an algebraic linear system of equations for au and aϕ, and hence, the approximated fields uh(x) and ϕh(x) are obtained.

The aforementioned approximation is suitable for rectangular/cuboidal geometries only. In order to overcome this limitation, we consider the immersed boundary method,95,96 where the arbitrarily-shaped domain Ω is embedded into a larger domain Ω with rectangular/cuboidal shape, i.e., ΩΩ. The geometrical map in Eq. (87) is redefined as

φ:ΞΩξx=φ(ξ),
(91)

which is independent on Ω, and hence, arbitrary geometries are allowed. The basis functions N(x)Ω are now defined on the embedding domain, which is discretized by means of a uniform Cartesian grid as Ω=cΩc, cf. Fig. 7(a). Hence, B-spline basis constitute a suitable approximation space on Ξ, and Eqs. (88)(90) hold.

FIG. 7.

Sketch of the immersed boundary method. (a) Physical domain ΩR2 (red) embedded in Ω (gray), which is discretized by a Cartesian grid into several cells Ωc. Inner cells are depicted in black stroke, cut cells in blue stroke, and outer cells are not depicted. (b) Detail of an inner cell. (c) Detail of a cut cell.

FIG. 7.

Sketch of the immersed boundary method. (a) Physical domain ΩR2 (red) embedded in Ω (gray), which is discretized by a Cartesian grid into several cells Ωc. Inner cells are depicted in black stroke, cut cells in blue stroke, and outer cells are not depicted. (b) Detail of an inner cell. (c) Detail of a cut cell.

Close modal
Remark 6:

Since the approximation space is unfitted to the geometry, it is not interpolant on Ω, even if B-splines of degree p=1 or any other family of basis functions are considered. Therefore, in the context of B-spline basis functions, knot multiplicity does not provide any benefit, and uniform knot vectors are typically considered instead.

The physical boundary Ω is allowed to intersect the cells Ωc of the embedding mesh arbitrarily, leading to a unfitted discretization of Ω. The cells with a nonempty intersection with Ω are labeled as active, whereas outer (or inactive) cells are neglected. Active cells can be either inner [ΩcΩ, see Fig. 7(b)] or cut [ΩcΩ, see Fig. 7(c)].

4. Cell classification

Labeling a cell depending on its intersection with Ω is usually accomplished by checking whether all vertices of each cell (and possibly more points within the cell) lie within the domain Ω (inner cell), only part of them (cut cell), or none of them (outer cell). In the case of implicit boundary representation (e.g., level set approaches), it is enough to evaluate the level set function on the evaluation points (see Fries and Omerović,100 Fries,101 Kudela et al.,102 Legrain et al.103). For explicit boundary representation (e.g., CAD descriptions), ray-tracing procedures are required, as explained in Marco et al.104,105

5. Enforcement of Dirichlet boundary conditions

Immersed boundary methods permit considering arbitrary geometries and involve trivial mesh generation, at the cost of having to deal with a non-conforming discretization. Since the approximation space is defined independently of Ω, it is not interpolant on its boundary Ω, hence, Dirichlet boundary conditions must be weakly enforced, e.g., by means of Nitsche’s method (cf. Sec. II C3).

6. Numerical integration

An additional challenge, associated with the presence of cut cells, consists on defining a good-enough numerical integration on cut cells. This implies defining special quadrature rules for each cell Ωc since only their intersection with Ω, i.e., ΩcΩ, is to be considered. The standard approach consists on dividing ΩcΩ into several non-overlapping sub-domains, each of them being easily integrated with predefined quadrature rules.

On the one hand, conforming subdivision approaches,101,102,104 most of them based on the marching cubes algorithm,106 attempt to subdivide cut cells into smaller conforming entities (e.g., squares or triangles in 2D, see Fig. 8), even allowing in some cases curved faces or edges.107–109 On the other hand, non-conforming subdivision schemes110,111 are based on quadtree/octree approaches, where each cell is recursively subdivided into 2D pieces until some predefined recursion level (see Fig. 9). Inner pieces are completely integrated, whereas at those that remain cut in the finest subdivision level, only the Gauss points that lie within Ω are considered.

FIG. 8.

Conforming sub-division of cut cells into non-overlapping triangles to perform numerical integration.

FIG. 8.

Conforming sub-division of cut cells into non-overlapping triangles to perform numerical integration.

Close modal
FIG. 9.

Non-conforming numerical integration of cut cells via recursive quadtree sub-division.

FIG. 9.

Non-conforming numerical integration of cut cells via recursive quadtree sub-division.

Close modal

Conforming approaches are typically more accurate but harder to implement, whereas non-conforming approaches are typically easier to implement and as accurate as required by increasing the recursivity.

Numerical integration on the faces and edges of cut cells requires an explicit boundary representation, available only if a CAD (or other explicit) description of Ω is considered. In the case of implicit boundary representations such as level sets, an explicit approximation of Ω must be precomputed at each cell.

7. Cut-cell stabilization

Another challenge related with cut cells consists on alleviating ill-conditioning produced by degrees of freedom with small intersection with Ω. The resulting algebraic system of equations suffers from ill-conditioning, i.e., large condition number of the system matrix, in the presence of cut cells with a small portion in the domain, i.e., when the volume fraction ηc=|ΩcΩ|/|Ωc|1. In particular, for the boundary value problems of flexoelectricity considered in this paper, the condition number scales with the minimum volume fraction ηmin=minc(ηc) in meshes of fixed size h at a rate of ηmin(2p+12/D),62 which implies that ill-conditioning is more severe for high-order basis (see Fig. 10).

FIG. 10.

Condition number against the minimum volume fraction using different cut-cell stabilization techniques for (a) cubic spline mesh (p=3) and (b) quartic spline mesh (p=4).

FIG. 10.

Condition number against the minimum volume fraction using different cut-cell stabilization techniques for (a) cubic spline mesh (p=3) and (b) quartic spline mesh (p=4).

Close modal

Several strategies have been proposed to alleviate ill-conditioning of trimmed cells, such as the ghost penalty method,112 the artificial stiffness approach,110,111 the extended B-spline method,113–116 or special preconditioning techniques specifically designed for immersed boundary methods,117 among others. We recommend here the extended B-spline approach by Höllig et al.113–116 due to its simplicity, ease of implementation, and good performance, cf. Fig. 10. The main idea consists on removing the critical basis functions (the ones with smaller intersection with Ω) from the approximation space and extrapolating well-behaved basis functions from neighboring cells toward the cut cell. Figure 11 illustrates this process for the univariate case. The modified basis has less degrees of freedom, but the condition number and error converge rates are independent of ηmin, and equivalent to those of body-fitted (untrimmed) methods.113 The extended B-spline basis stabilization can be easily implemented as a linear constraint on the approximation space of cut cells based on the uniform Cartesian structure of the discretization.62,113–116

FIG. 11.

Extended B-spline approach on the right end of a univariate mesh of degree p=2. The resulting functional space contains the extended basis functions {B~32(ξ), B~42(ξ), B~52(ξ)}, as a result of the extrapolation of the basis {B32(ξ), B42(ξ), B52(ξ)} in Ω4 toward Ω5. The critical basis B62(ξ) is removed.

FIG. 11.

Extended B-spline approach on the right end of a univariate mesh of degree p=2. The resulting functional space contains the extended basis functions {B~32(ξ), B~42(ξ), B~52(ξ)}, as a result of the extrapolation of the basis {B32(ξ), B42(ξ), B52(ξ)} in Ω4 toward Ω5. The critical basis B62(ξ) is removed.

Close modal

So far, uniform meshes have been considered. However, in practice, spatial h-adaptivity is required to increase the resolution of the approximation where needed. In the context of multivariate B-spline bases, refining the knot vectors around a precise location is not effective since the tensor product structure of the basis yields a delocalized refinement that spans along the whole mesh.

Alternatively, one can resort to hierarchical B-spline refinement,118 where the approximation space is locally enriched by replacing selected coarse B-spline basis functions (parents) with finer ones (children), see Fig. 12. This process can be performed recursively, leading to a parent-children hierarchy spanning several levels of refinement, where each refinement level is defined by its own uniform knot vector set.

FIG. 12.

Univariate uniform cubic (p=3) hierarchical B-spline basis. The basis functions Bi3,1(ξ1), i=3,,6 in the coarsest level of refinement (=1) have been replaced by their children Bj3,2(ξ2), j=5,,15 in the next (=2) refinement level.

FIG. 12.

Univariate uniform cubic (p=3) hierarchical B-spline basis. The basis functions Bi3,1(ξ1), i=3,,6 in the coarsest level of refinement (=1) have been replaced by their children Bj3,2(ξ2), j=5,,15 in the next (=2) refinement level.

Close modal

This process is based on a remarkable property of uniform B-splines: their natural refinement by subdivision, by which a uniform B-spline function can be expressed as a linear combination of contracted, translated and scaled copies of itself (see Fig. 13). Mathematically, the subdivision property is expressed as the two-scale relation119 

Bip,(ξ)=j=1p+2[sp]jB2(i1)+jp,(2ξ)=j=2i12i+p[sp]j2(i1)Bjp,+1(ξ+1),
(92)

where N stands for the refinement level, ξ+1=2ξ and the scaling coefficients

[sp]j=12p(p+1j1)=2p(p+1)!(j1)!(p+2j)!,
(93)

for j={1,,p+2}. This process maintains the properties of the spanned functional space and, in particular, its regularity.

FIG. 13.

Two-scale relation of univariate B-splines. Top: original (parent) B-spline B1p,1(ξ1) at level =1. Bottom: the children B-spline basis functions Bjp,2(ξ2) at level =2 scaled by the factors [sp]j, for j={1,,p+2}, (a) linear (p=1), s1=12[1,2,1], (b) quadratic (p=2), s2=14[1,3,3,1], (c) cubic (p=3), s3=18[1,4,6,4,1].

FIG. 13.

Two-scale relation of univariate B-splines. Top: original (parent) B-spline B1p,1(ξ1) at level =1. Bottom: the children B-spline basis functions Bjp,2(ξ2) at level =2 scaled by the factors [sp]j, for j={1,,p+2}, (a) linear (p=1), s1=12[1,2,1], (b) quadratic (p=2), s2=14[1,3,3,1], (c) cubic (p=3), s3=18[1,4,6,4,1].

Close modal

In practical computations, the goal is to refine some cells in the mesh. In a hierarchical B-spline context, it corresponds to refining some of the non-vanishing basis functions on those cells. There exist different hierarchical refinement strategies, depending on the relation between cells and basis to be refined,120–124 yielding to more or less localized refinement.

A subtle implementation detail of major practical relevance is considering the nondimensionalization of all the variables and parameters in the equations, specially since one is dealing with a coupled problem with multiple physics. To understand its relevance, think of the Kuu and Kϕϕ blocks in the diagonal of the system matrix K. If one considers the International System of Units (SI), their magnitudes are expressed in Pa and F/m, corresponding to Young’s modulus Y and electric permittivity ε of the material. In a typical dielectric, they are in the order of 109 (or higher) and 109 (or lower), respectively, which makes at least 18 orders of magnitude difference! This would preclude the solution of the numerical problem due to a very large condition number of K.

To avoid this issue, a simple strategy consists on dividing all the variables in the problem by some normalizing factors. In the case of flexoelectricity boundary value problems, there are three independent magnitudes, e.g., length (or displacement), stress, and polarization. Defining some suitable normalizing factors fL, fS, and fP for each of them, the normalized aforementioned magnitudes are defined as

u0=u/fL,σ0=σ/fS,P0=P/fP.
(94)

All the other quantities in the problem are divided by normalization factors that are derived as a combination of fL, fS, and fP, e.g.,

  • the electric potential ϕ0=ϕ/(fLfSfP1),

  • Young’s modulus Y0=Y/fS,

  • the flexoelectric coefficients μ0=μ/(fPfL),

  • the electric permittivity ε0=ε/(fP2fS1),

and so on. Typically, suitable values for the normalizing factors are fLLgeom, fSY, and fPYε, where Lgeom is a characteristic length of the geometry.

Remark 7:

The normalized flexoelectric coefficients μ0 scale inversely proportional to the size of the geometry, whereas the normalized Young modulus and electric permittivity remain constant. This fact reveals the well-known size-dependent nature of the flexoelectric coupling, relevant only at sub-micrometer scales.

In this section, we make use of the numerical methods described in Sec. III to solve different boundary values problems (BVP) stated in Sec. II C. We focus on cantilever beam bending, the most well-known benchmark for flexoelectricity, widely used by experimentalists to characterize the transversal flexoelectric effect.85,87,125–127 It has also been studied numerically58,59,63 and analytically.128,129 In the experiments in Secs. IV A and IV B, we analyze the electromechanical response of microscopic cantilevers under applied force and applied electric potential in closed circuit, respectively, considering the Direct and Lifshitz-invariant flexoelectricity formulations.

We consider a 2D (plane strain) cantilever beam of length L=8μm and thickness H=0.4μm. The material properties are simple enough to isolate the transversal flexoelectric effect, i.e., a Young modulus Y=100 GPa, electric permittivity ε=11 nC/V m, and transversal flexoelectric coefficient μT=1μC/Vm. The other material parameters are set to 0. For a complete description of material tensors, we refer to  Appendix A.

The left tip is clamped and a vertical force F=1μN/μm is applied on the top right corner. The right tip is electrically grounded and the other boundaries are free, which corresponds to open-circuit electrical boundary conditions [Fig. 14(a)]. The transversal flexoelectric effect is triggered due to the mechanically-induced gradient of axial strains along the beam’s cross section.

FIG. 14.

Cantilever beam bending using direct and Lifshitz-invariant flexoelectricity forms. (a) Geometry and boundary conditions. (b) Electric potential distribution within the beams. (c) Axial strains at the cross section x=L/2. (d) Transversal electric fields at the cross section x=L/2.

FIG. 14.

Cantilever beam bending using direct and Lifshitz-invariant flexoelectricity forms. (a) Geometry and boundary conditions. (b) Electric potential distribution within the beams. (c) Axial strains at the cross section x=L/2. (d) Transversal electric fields at the cross section x=L/2.

Close modal

The problem is numerically solved by means of the body-fitted B-spline method in Sec. III. The results are shown in Fig. 14. As pointed out in Sec. II A, the Lifshitz-invariant and direct flexoelectricity boundary value problems are not equivalent, even though their associated Euler–Lagrange Eqs. (53) and (65) coincide. The difference is due to the different boundary terms (e.g., tractions and surface charges) in Eqs. (54) and (66) given the different expressions in (49) and (64). In particular, traction-free or charge-free boundary conditions have different meaning in each formulation and thus the resulting fields are different.

Two main differences are pointed out next. First, the mechanical results are quite similar. The axial strains vary linearly along the cross sections of the bent beam in both cases [Fig. 14(c)], as expected. However, the flexoelectricity-induced stiffening of the beam41 is different in each case. Comparing the maximum deflection of the direct-flexoelectric beam (0.30μm) and the Lifshitz-invariant-flexoelectric beam (0.24μm) with respect to a standard elastic one (0.32μm), it becomes apparent that the effective stiffness is increased around 7% in the former and a 33% in the latter.

Nevertheless, the most interesting difference arises in the electrical response, as shown in Figs. 14(b) and 14(d). While the direct flexoelectricity form presents an electric potential varying linearly along a cross section of the beam, the Lifshitz-invariant form features boundary layers with opposite sign than that of the bulk. This phenomenon is highlighted in Fig. 14(d), which depicts the transversal electric fields along the middle cross section.

The electric response in the Lifshitz-invariant model resembles the predicted theoretical response of flexoelectricity in the presence of surface piezoelectricity.35 Since the fundamental difference between the two studied BVPs is just the null Lagrangian in Eq. (1), which is also acting on the surface, a relation may exist between the modeling of null Lagrangians and surface effects in this context.

Remark 8:

As discussed later on in  Appendix B, the appearance of a boundary layer may lead to mesh-dependent spurious oscillations in the numerical solutions. Hence, some stabilization technique is required in that case. The simulations presented in Secs. IV A and IV B are properly stabilized according to the approach described in  Appendix B.

In this experiment, we explore the transversal flexoelectric effect triggered by electrical actuation. This device was first used by Bursian and Zaikovskii3 to experimentally demonstrate for the first time the flexoelectric effect, which had been predicted theoretically by Mashkevich and Tolpygo.11 Computational studies are also present in Abdollahi et al.58 and Zhuang et al.59 

We consider the same geometry and material properties as in previous experiment. Here, an electric field across the beam thickness is enforced by attaching an electrode on the top boundary at prescribed voltage V=5 V, while the bottom boundary is grounded, which corresponds to closed-circuit electrical boundary conditions [Fig. 15(a)]. Mechanically, the left tip is clamped, and no force is applied. Due to the transversal flexoelectric effect, the electric field will yield an axial strain gradient along the thickness of the beam, inducing a constant curvature.

FIG. 15.

Cantilever beam actuator using direct and Lifshitz-invariant flexoelectricity forms. (a) Geometry and boundary conditions. (b) Electric potential at the cross section x=L/2. (c) Axial strains at the cross section x=L/2. (d) Transversal electric fields at the cross section x=L/2.

FIG. 15.

Cantilever beam actuator using direct and Lifshitz-invariant flexoelectricity forms. (a) Geometry and boundary conditions. (b) Electric potential at the cross section x=L/2. (c) Axial strains at the cross section x=L/2. (d) Transversal electric fields at the cross section x=L/2.

Close modal

The results are shown in Fig. 15 and are quite similar to the ones reported in the cantilever bending case. The two differences between the direct and Lifshitz-invariant flexoelectricity models are also present here. The latter presents more stiffening, in view of the maximum deflections obtained: 0.30μm for the direct case and only 0.12μm for the Lifshitz-invariant one. Boundary layers in the electric field distribution are also obtained here for the Lifshitz-invariant form [Fig. 15(d)]. However, in this case the bulk electric field is much larger than the boundary layer effect, and hence the electric potential distributions are much more alike [Fig. 15(b)].

In this section, we illustrate the richer physics in the regime of finite deformations by modeling a soft cantilever rod actuator.41 This device was considered by Bursian and Zaikovskii3 back in 1968 to experimentally demonstrate for the first time the flexoelectric effect.

Let us consider the direct flexoelectricity boundary value problem stated in Sec. II C 4 c in a compliant slender rod of length L=2μm and thickness H=100nm, see Fig. 16(a). The left tip of the rod is clamped, while all other boundaries are traction-free. The electric potential Φ at the top boundary is set to a certain value V, and the bottom boundary is grounded (Φ=0). The voltage drop induces a transverse electric field EY=V/H across the rod thickness which triggers (i) the flexoelectric effect, thereby generating a non-uniform axial strain that bends the rod, and (ii) the electrostrictive effect, yielding a uniform axial strain that elongates the rod. The induced curvature and axial stretch are constant throughout the neutral axis of the rod. Since the flexoelectric effect is reversible in sign, the cantilever will bend upwards or downwards depending on the direction of the applied electric field. However, since the electrostrictive effect is not, the beam length will always increase regardless of the electric field direction.

FIG. 16.

Actuation of soft cantilever rod with different flexoelectric tensors μ (nC/m). (a) Geometry and boundary conditions. (b) Axial strain ζ(V). (c) Curvature R1(V). (d) Deformed configuration and electric potential distribution in case A upon increasing voltage V (kV), indicated by the number at the free end. This figure was adapted from Codony et al., J. Mech. Phys. Solids 146, 104182 (2020). Copyright 2020 Elsevier.

FIG. 16.

Actuation of soft cantilever rod with different flexoelectric tensors μ (nC/m). (a) Geometry and boundary conditions. (b) Axial strain ζ(V). (c) Curvature R1(V). (d) Deformed configuration and electric potential distribution in case A upon increasing voltage V (kV), indicated by the number at the free end. This figure was adapted from Codony et al., J. Mech. Phys. Solids 146, 104182 (2020). Copyright 2020 Elsevier.

Close modal

To solve the resulting non-linear system of equations, the voltage V is applied incrementally in a sequence of load steps. We resort to a modified-step monolithic Newton–Raphson algorithm presented in Codony et al.,41 which is robust against the nonlinearities in the formulation, including the eventual presence of mechanical (geometrical) instabilities such as wrinkling, creasing, or buckling.

Figures 16(b)16(d) show the electromechanical response of the actuated rod modeled by an isotropic hyperelastic Neo-Hookean law augmented with strain gradient elasticity as follows:

ΨMech(E,0E)=λ2([log(J)]2+2EII,KEJJ,K)+μ(Tr(E)+2EIJ,KEIJ,K),
(95)

with Young’s modulus Y=1.0 GPa, Poisson’s ratio ν=0.37 (that is, Lamé’s moduli λ1.039 GPa and μ=0.685 GPa), mechanical length scale =30 nm, dielectric permittivity ε=0.11 nC/V m and different combinations of the flexoelectric constants. The results with negative curvature have been omitted, and the effect of varying the shear flexoelectric constant μS has been reported negligible for the studied cases.41 As expected, the axial strain of the rod [depicted in Fig. 16(b)] does not vary much with the different flexoelectric parameters, since it is mainly a consequence of electrostriction. The curvature [Fig. 16(c)], instead, varies significantly for the different combinations of flexoelectric parameters. The largest response is found with positive μT and negative μL, as shown in case A. For sufficiently large actuation, the rod is even able to roll up forming a closed circle [Fig. 16(d)].

We analyze next several devices and metamaterials exploiting the flexoelectric effect through the effective accumulation of the response to gradients in the mechanical and electric fields in the limit of infinitesimal deformation. The example in Sec. VI A shows a scalable flexoelectricity-based deformation sensor which works upon collective beam-bending, resulting in a potential difference that accumulates from beam to beam. In Secs. VI B and VI C, we simulate the RVE corresponding to different engineered, geometrically-polarized architected materials and show that they feature apparent piezoelectric behavior without piezoelectric constituents, allowing its use in sensing or actuation applications.

Previous works have shown that bending of thin structural elements is the most efficient way to mobilize flexoelectricity.9,75 We present here a setup achieving collective beam bending in such a way that the flexoelectrically-generated electric potential at the thin structural elements can be effectively accumulated throughout the structure and collected at the electrodes. The structure is composed of several thin beams connected by vertical elements on one end and to other beams through insulating connectors at the other end. The insulating connectors are placed as to achieve a non-centrosymmetric system and thus avoid internal cancelation of flexoelectrically-generated electric potential of opposite signs at the beams.83 This unit is appropriately repeated in series, producing a scalable device, here shown in 2D in Fig. 17(a). We analyze such a structure, with beam thicknesses of t=40 nm and lengths of =400 nm, and beam spacings of h=280 nm, under three point displacement as shown in Fig. 17(a), with δ=40 nm. The direct formulation for flexoelectricity in Sec. II C1 is used, with the interface conditions in Sec. II C42 at the dielectric-insulator interfaces. The material parameters are reported in Table I.

FIG. 17.

Collective-beam bending device. The applied displacement pattern induces strain gradients in the structure while the material arrangement breaks the overall centrosymmetry of the system, resulting in an accumulation of the electric potential. (a) Geometrical configuration. Flexoelectric material is depicted in orange, and the insulating material is depicted in gray. The arrows represent a vertical displacement of δ=40 nm. (b) Deformed shape and resulting electric potential distribution.

FIG. 17.

Collective-beam bending device. The applied displacement pattern induces strain gradients in the structure while the material arrangement breaks the overall centrosymmetry of the system, resulting in an accumulation of the electric potential. (a) Geometrical configuration. Flexoelectric material is depicted in orange, and the insulating material is depicted in gray. The arrows represent a vertical displacement of δ=40 nm. (b) Deformed shape and resulting electric potential distribution.

Close modal
TABLE I.

Material parameters in Sec. VI A.

Material ν mech ɛ μL μT μS 
 (GPa) — (nm) (nC/V m) (μC/m) (μC/m) (μC/m) 
Matrix 152 0.33 141 150 110 110 
Insulators 152 0.33 
Material ν mech ɛ μL μT μS 
 (GPa) — (nm) (nC/V m) (μC/m) (μC/m) (μC/m) 
Matrix 152 0.33 141 150 110 110 
Insulators 152 0.33 

Figure 17(b) shows the electric potential, ϕ, distribution. The accumulation of electric potential from the grounded (bottom-left) toward the top-right edge is apparent. This electric potential accumulation scales linearly with the number of connected beams.

According to Sharma et al.75 and Mocci et al.,83 geometrically-polarized inclusions embedded in a non-piezoelectric matrix produce an apparently piezoelectric response upon macroscopic homogeneous deformation. This response results from the effective accumulation of the flexoelectric response to the localized strain gradients. With this idea, one can endow any dielectric with apparent piezoelectricity.83 

We illustrate this idea by simulating the electromechanical response of flexoelectric metamaterial with triangular voids under homogeneous macroscopic vertical compression, Fig. 18. We consider the Lifshitz-invariant flexoelectric formulation, and barium strontium titanate (BST) in its paraelectric phase as base material, see Table II for material parameters. Paraelectric BST is a non-piezoelectric dielectric with high permittivity, and hence, a good flexoelectric.35 

FIG. 18.

Sensor device with geometrically polarized inclusions embedded in a non-piezoelectric matrix. The electric potential distribution is shown and plotted along the dashed line, both for the whole domain (thin line) and a single unit cell calculation (thick line). The electric potential has been normalized as ϕ^=ϕε/μ.

FIG. 18.

Sensor device with geometrically polarized inclusions embedded in a non-piezoelectric matrix. The electric potential distribution is shown and plotted along the dashed line, both for the whole domain (thin line) and a single unit cell calculation (thick line). The electric potential has been normalized as ϕ^=ϕε/μ.

Close modal
TABLE II.

Material parameters of BST in Secs. VI B and VI C.

 ν mech ɛ elec μL μT μS 
 (GPa— (nm) (nC/V m) (nm) (μC/m) (μC/m) (μC/m) 
 152 0.33 50 300 1.2 1.1 0.05 
 ν mech ɛ elec μL μT μS 
 (GPa— (nm) (nC/V m) (nm) (μC/m) (μC/m) (μC/m) 
 152 0.33 50 300 1.2 1.1 0.05 

To efficiently evaluate the performance of the periodic structure, high-order generalized periodic conditions in Eq. (76) are enforced in a single unit cell. Standard periodicity for the solution fields along the horizontal direction is applied ([[u]]x1=[[ϕ]]x1=0), whereas the jump in electric potential along the vertical direction ([[ϕ]]x2) is resolved as a result of a prescribed vertical jump in the displacement field ([[u]]x2=[[u]]x2¯).

Figure 18 (middle) shows the resulting electric potential distribution on the deformed sensor. The plot on the right represents the trend of the electric potential along the vertical dashed line by computing a portion of the whole domain (thin black line) and a single unit cell with a height of 3.15μm with generalized periodicity conditions (thick magenta line). Both simulations produce exactly the same result on a unit cell far away from the top and bottom boundaries. A non-vanishing net potential difference between the top and bottom layers of the device is achieved upon homogeneous deformation. The observed apparent piezoelectricity is the result of the accumulation of the flexoelectric response to strain gradients generated around the voids. This accumulation is in turn possible because the triangular voids produce a geometrical polarization of the microstructure. Such a material exhibits also inverse apparent piezoelectricity under the application of a homogeneous electrical bias.83 

Building on the ideas discussed in Secs. VI A and VI B, a new class of flexoelectricity-based metamaterials have been recently proposed by Mocci et al.83 These metamaterials are bending dominated, low area fraction lattices, with a non-centrosymmetric arrangement of small non-piezoelectric dielectric beams. Upon homogeneous macroscopic deformation or electric bias, the metamaterial behaves as an apparent piezoelectric due to the accumulation of the flexoelectric response of the lattice beam components. Mocci et al.83 show that it can reach performances comparable in some situations to those of very common piezoelectric materials, such as quartz or PZT.

Figure 19 shows a bending-dominated BST lattice metamaterial in actuation mode, i.e., deforming due to an applied macroscopic electric bias ([[ϕ]]x2=[[ϕ]]x2¯), modeled by the Lifshitz-invariant flexoelectric formulation. A RVE of the lattice is analyzed using generalized periodicity conditions on the top and bottom faces and standard periodicity in the horizontal direction. The lattice unit cell is has dimensions L1=4.8μm and L2=2.8μm and the beam thickness is t=160 nm. The resulting vertical displacement ([[u]]x2) and the deformed configuration are also depicted in a reduced portion of the whole domain. The metamaterial behaves again as an apparent piezoelectric despite its base material is not.

FIG. 19.

Bending-dominated lattice metamaterial device in the actuation mode. The electric field-induced displacement is depicted in a reduced portion of the lattice. The electric potential has been normalized as ϕ^=ϕε/μ, and the displacement field as u^=u/L2.

FIG. 19.

Bending-dominated lattice metamaterial device in the actuation mode. The electric field-induced displacement is depicted in a reduced portion of the lattice. The electric potential has been normalized as ϕ^=ϕε/μ, and the displacement field as u^=u/L2.

Close modal

We have reviewed the mathematical modeling of the flexoelectric effect in the regime of infinitesimal deformations, including different flexoelectric couplings, variational principles, and boundary value problems and presented its natural extension to the regime of finite deformations. We have also described in detail an efficient computational approach based on immersed-boundary hierarchical B-splines, enabling the computation on arbitrarily-shaped geometries with attractive numerical properties such as trivial mesh generation and optimal convergence rates. The provided numerical examples illustrate the mechanisms behind the flexoelectric effect and show the feasibility of engineered devices and materials designed for a variety of electromechanical applications. More specifically, we have shown that the Lifshitz-invariant formulation leads to boundary layers, and we have proposed material architectures capable of mobilizing the flexoelectric effect at a microstructural level and upscaling this effect to a mesoscale.

This paper exemplifies how accurate computational methods can enable new material designs that exploit flexoelectricity in engineering applications, for instance, to yield apparent piezoelectricity in metamaterials made of non-piezoelectric base materials,83 given the inherent complexity of setups mobilizing field gradients. With this background, we expect that materials designs can be optimized geometrically and topologically,44,130,131 and that ultimately the fundamental design principles of flexoelectric metamaterials can be distilled. There are some aspects that still hinder the quantitative evaluation of flexoelectricity. A clear physical interpretation of the high-order fields is still missing, which would allow us to solve the ambiguity regarding the high-order boundary conditions on physical grounds. Continuum models would benefit from a more precise connection with first principles. This would allow us not only to identify parameters for the continuum descriptions but also learn proper enrichments in order to capture additional physics, e.g., finite sample effects. These important effects, such as surface piezoelectricity and surface flexoelectricity, have not been included here.35 Finally, a wider and more precise experimental quantification of flexoelectricity in an ample catalog of materials would also be desirable.

Avenues for future research in the mathematical and computational modeling of flexoelectricity are vast. One of the most appealing ones is the complete understanding and exploitation of flexoelectricity in soft materials. Despite larger strain gradients (hence, larger flexoelectric response) are easily possible together with functional deformation modes, the fundamental mechanisms leading to flexoelectricity in polymers remain unclear.36 Other interesting research lines involve the understanding and exploitation of flexoelectricity in other important multifunctional materials, such as magnetoelectrics, liquid crystal elastomers, ferroelectrics, or semiconductors, with multiple potential applications.

This work was supported by the European Research Council (No. StG-679451 to I.A.), the Spanish Ministry of Economy and Competitiveness (No. RTI2018-101662-B-I00), and the Generalitat de Catalunya (ICREA Academia award for excellence in research to I.A. and Grant No. 2017-SGR-1278 to I.A.). CIMNE is a Severo Ochoa Centre of Excellence (2019-2023) under the grant CEX2018-000797-S funded by MCIN/AEI/10.13039/501100011033.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors have no conflicts to disclose.

The material enthalpy density in an infinitesimal deformation framework is characterized by the tensors of elasticity c, strain gradient elasticity h, dielectricity κ, flexoelectricity μ, and eventually gradient dielectricity M.

a. Elasticity tensor c

Isotropic elasticity is represented by the fourth-order tensor c, which depends on λ and μ in the following form:

cIJKL=λδIJδKL+2μδIKδJL.
(A1)
Remark 9:

The Lamé parameters λ and μ are related to the Young modulus Y and Poisson’s ratio ν as λ=Yν/(1+ν)(12ν) and μ=Y/2(1+ν).

b. Strain gradient elasticity tensor h

We consider an isotropic simplified strain gradient elasticity tensor,132 which depends on λ, μ and the mechanical length scale mech in the following form:

hIJKLMN=(λδIJδLM+2μδILδJM)mech2δKN.
(A2)

c. Dielectricity tensor κ

Isotropic dielectricity is represented by the second-order tensor κ, which depends on the electric permittivity ε as

κij=εδij.
(A3)

d. Gradient dielectricity tensor M

Isotropic gradient dielectricity is represented by the fourth-order tensor M. We take a simple form depending on the electric permittivity ε and the dielectric length scale elec as

Mijkl=εelec2δikδjl.
(A4)

e. Flexoelectricity tensor μ

The cubic flexoelectric tensor depends on the longitudinal μL, transversal μT, and shear μS parameters.63,133 When oriented along the Cartesian coordinates, it takes the following form:

μ<x1>LIJK={μL,forL=I=J=K,μT,forI=JK=L,μS,forL=IJ=KorL=JI=K,0otherwise.
(A5)

The flexoelectric tensor μ oriented in an arbitrarily rotated orthonormal basis is obtained by rotating μx1 as

μlijk=RlLRiIRjJRkKμx1LIJK,
(A6)

where R is the rotation matrix associated with the unit vectors of the rotated basis.

The isotropic flexoelectric tensor is a particular case of the cubic one with only two independent parameters, where μL=μT+2μS.

As discussed in Remark 8, depending on the geometry of Ω, the boundary conditions, the material parameters, and the mesh size, the numerical solution to the flexoelectricity problem features spurious oscillations that completely spoil the quality of the results. Figure 20 shows an example, with spurious oscillations arising in a Lifshitz-form flexoelectric cantilever beam with open circuit boundary conditions undergoing bending (i.e., a similar setting to Sec. IV A). The beam is L=3.2μm long by H=0.4μm thick. The dielectric material has Young’s modulus Y=100 GPa, Poisson’s ratio ν=0.37, dielectric permittivity ε=11nJ/V2m, and longitudinal and transversal flexoelectric coefficients μL=μT=1μJ/Vm. Other material coefficients are 0. A force of F=100μN is applied vertically on the top of the right end.

FIG. 20.

Spurious, mesh-dependent oscillatory behavior of the electric potential ϕ in a Lifshitz-invariant flexoelectric cantilever beam. The B-spline basis is quadratic (p=2) with mesh size h=0.0275μm. (a) Electric potential distribution within the cantilever. (b) Detail of the electric potential along the cross section at x=L/2.

FIG. 20.

Spurious, mesh-dependent oscillatory behavior of the electric potential ϕ in a Lifshitz-invariant flexoelectric cantilever beam. The B-spline basis is quadratic (p=2) with mesh size h=0.0275μm. (a) Electric potential distribution within the cantilever. (b) Detail of the electric potential along the cross section at x=L/2.

Close modal

In view of the results in Secs. IV A and IV B, and in agreement with the experience of this paper’s authors, the numerical instabilities are associated with the presence of boundary layers in the electric field, which may frequently appear when considering the Lifshitz-invariant flexoelectricity formulation.

Remark 10:

In the literature of computational flexoelectricity, the direct flexoelectricity form is much more popular than others, specially for the computational examples. Hence, it is difficult to say whether numerical instabilities are also present in other numerical schemes rather than smooth B-spline approximations.

A complete and exhaustive study of numerical instabilities and the proposal of stabilization strategies is a matter of current research of this paper’s authors and is out of the scope of this tutorial. However, it is instructive thinking of a simplified 1D-version of the Euler–Lagrange equations associated with the Lifshitz-invariant flexoelectricity form [cf. Eq. (65)]

{Luuu+LuϕϕYuIIYmech2uIVμϕIII=b,Lϕuu+LϕϕϕεϕII+εelec2ϕIV+μuIII=q,
(B1)

depending on just five material parameters, namely, the Young modulus Y and mechanical length scale mech, the dielectric constant ε and the dielectric length scale elec, and the flexoelectric coefficient μ.

The corresponding weak form can be written as

{auu(δu,u)+auϕ(δu,ϕ)=lu(δu),aϕϕ(δϕ,ϕ)+auϕ(u,δϕ)=lϕ(δϕ),
(B2)

with appropriate boundary conditions, where

auu(δu,u)=(δuI,YuI)Ω+(δuII,Ymech2uII)Ω,
(B3a)
aϕϕ(δϕ,ϕ)=(δϕI,εϕI)Ω+(δϕII,εelec2ϕII)Ω,
(B3b)
auϕ(δu,ϕ)=12μ((δuI,ϕII)Ω(δuII,ϕI)Ω),
(B3c)
lu(δu)=(δu,b)Ω,
(B3d)
lϕ(δϕ)=(δϕ,q)Ω.
(B3e)

To stabilize this simple 1D model, we resort to the Galerkin least squares (GLS) method134–136 due to its simple form and implementation. Following the GLS method, we make use of Eq. (B1) to define the following bilinear forms:

Amech({u,ϕ},{δu,δϕ})=(Luuδu+Luϕδϕ,τmechrmech(u,ϕ))Ω,
(B4a)
Aelec({u,ϕ},{δu,δϕ})=(Lϕuδu+Lϕϕδϕ,τelecrelec(u,ϕ))Ω,
(B4b)

with the residuals

rmech(u,ϕ)=Luuu+Luϕϕ+b,
(B5a)
relec(u,ϕ)=Lϕuu+Lϕϕϕq.
(B5b)

Note that a solution fulfilling Eq. (B1) yields vanishing residuals in Eq. (B5) and, therefore, the bilinear forms in Eq. (B4) vanish too. Hence, they can be added to the original weak form in Eq. (B2) while maintaining self-consistency

auu(δu,u)+auϕ(δu,ϕ)+Amech({u,ϕ},{δu,0})+Aelec({u,ϕ},{δu,0})=lu(δu),
(B6a)
aϕϕ(δϕ,ϕ)+auϕ(u,δϕ)+Amech({u,ϕ},{0,δϕ})+Aelec({u,ϕ},{0,δϕ})=lϕ(δϕ).
(B6b)

The stabilization parameters τmech,τelec provide control on the second derivatives of the state variables and can be tuned as a function of mesh size and material parameters. A reasonable choice for the stabilization parameters can be found in Codony62 as

τmech=αmechh2Y,τelec=αelech2ε,
(B7)

where h represents the mesh size and αmech,αelec are dimensionless scalars that can either be chosen as constant or dependent on the relation between material parameters and mesh size. By choosing large enough αmech and αelec, spurious oscillations in the numerical solution can be prevented.

A 2D implementation of the aforementioned GLS stabilization is illustrated in Fig. 21, showing the control on the spurious numerical oscillations as the stabilization parameters are tuned, and a robust approximation of {u,ϕ} for large enough stabilization parameters.

FIG. 21.

Residual-based weak form stabilization in a Lifshitz-invariant flexoelectric cantilever beam. The B-spline basis is quadratic (p=2) with a mesh size of h=0.0275μm. The stabilization parameters are chosen as αmech=αelec=α.

FIG. 21.

Residual-based weak form stabilization in a Lifshitz-invariant flexoelectric cantilever beam. The B-spline basis is quadratic (p=2) with a mesh size of h=0.0275μm. The stabilization parameters are chosen as αmech=αelec=α.

Close modal
1.
W.
Ma
and
L. E.
Cross
, “
Observation of the flexoelectric effect in relaxor Pb (Mg1/3Nb2/3)O3 ceramics
,”
Appl. Phys. Lett.
78
,
2920
2921
(
2001
).
2.
L.
Cross
, “
Flexoelectric effects: Charge separation in insulating solids subjected to elastic strain gradients
,”
J. Mater. Sci.
41
,
53
63
(
2006
).
3.
J.
Bursian
and
O.
Zaikovskii
, “
Changes in curvature of a ferroelectric film due to polarization
,”
Sov. Phys. Solid State
10
,
1121
1124
(
1968
).
4.
É. V.
Bursian
and
N. N.
Trunov
, “
Nonlocal piezoelectric effect
,”
Sov. Phys. Solid State
16
,
760
762
(
1974
).
5.
A.
Abdollahi
and
I.
Arias
, “
Constructive and destructive interplay between piezoelectricity and flexoelectricity in flexural sensors and actuators
,”
J. Appl. Mech.
82
,
121003
(
2015
).
6.
U. K.
Bhaskar
,
N.
Banerjee
,
A.
Abdollahi
,
Z.
Wang
,
D. G.
Schlom
,
G.
Rijnders
, and
G.
Catalan
, “
A flexoelectric microelectromechanical system on silicon
,”
Nat. Nanotechnol.
11
,
263
266
(
2016
).
7.
A.
Abdollahi
,
N.
Domingo
,
I.
Arias
, and
G.
Catalan
, “
Converse flexoelectricity yields large piezoresponse force microscopy signals in non-piezoelectric materials
,”
Nat. Commun.
10
,
1266
(
2019
).
8.
J. Y.
Fu
,
W.
Zhu
,
N.
Li
, and
L. E.
Cross
, “
Experimental studies of the converse flexoelectric effect induced by inhomogeneous electric field in a barium strontium titanate composition
,”
J. Appl. Phys.
100
,
024112
(
2006
).
9.
A.
Abdollahi
,
D.
Millán
,
C.
Peco
,
M.
Arroyo
, and
I.
Arias
, “
Revisiting pyramid compression to quantify flexoelectricity: A three-dimensional simulation study
,”
Phys. Rev. B
91
,
104103
(
2015
).
10.
S. M.
Kogan
, “
Piezoelectric effect during inhomogeneous deformation and acoustic scattering of carriers in crystals
,”
Sov. Phys. Solid State
5
,
2069
2070
(
1964
).
11.
V. S.
Mashkevich
and
K. B.
Tolpygo
, “
Electrical, optical and elastic properties of diamond type cristals. 1
,”
Sov. Phys. Solid State
5
,
435
439
(
1957
).
12.
K.
Tolpygo
, “
Long wavelength oscillations of diamond-type crystals including long range forces
,”
Sov. Phys. Solid State
4
,
1297
1305
(
1963
).
13.
A. K.
Tagantsev
, “
Piezoelectricity and flexoelectricity in crystalline dielectrics
,”
Phys. Rev. B
34
,
5883
5889
(
1986
).
14.
A. K.
Tagantsev
, “
Electric polarization in crystals and its response to thermal and elastic perturbations
,”
Phase Transit. A
35
,
119
203
(
1991
).
15.
R. D.
Mindlin
, “
Polarization gradient in elastic dielectrics
,”
Int. J. Solids Struct.
4
,
637
642
(
1968
).
16.
E.
Sahin
and
S.
Dost
, “
A strain-gradients theory of elastic dielectrics with spatial dispersion
,”
Int. J. Eng. Sci.
26
,
1231
1245
(
1988
).
17.
R.
Maranganti
,
N.
Sharma
, and
P.
Sharma
, “
Electromechanical coupling in nonpiezoelectric materials due to nanoscale nonlocal size effects: Green’s function solutions and embedded inclusions
,”
Phys. Rev. B
74
,
014110
(
2006
).
18.
R. D.
Mindlin
and
N. N.
Eshel
, “
On first strain-gradient theories in linear elasticity
,”
Int. J. Solids Struct.
4
,
109
124
(
1968
).
19.
H.
Askes
and
E. C.
Aifantis
, “
Gradient elasticity in statics and dynamics: An overview of formulations, length scale identification procedures, finite element implementations and new results
,”
Int. J. Solids Struct.
48
,
1962
1990
(
2011
).
20.
R. D.
Mindlin
and
H. F.
Tiersten
, “
Effects of couple-stresses in linear elasticity
,”
Arch. Ration. Mech. Anal.
11
,
415
448
(
1962
).
21.
A. R.
Hadjesfandiari
, “
Size-dependent piezoelectricity
,”
Int. J. Solids Struct.
50
,
2781
2791
(
2013
).
22.
R.
Poya
,
A. J.
Gil
,
R.
Ortigosa
, and
R.
Palma
, “
On a family of numerical models for couple stress based flexoelectricity for continua and beams
,”
J. Mech. Phys. Solids
125
,
613
652
(
2019
).
23.
L.
Anqing
,
Z.
Shenjie
,
Q.
Lu
, and
C.
Xi
, “
A flexoelectric theory with rotation gradient effects for elastic dielectrics
,”
Modell. Simul. Mater. Sci. Eng.
24
,
015009
(
2015
).
24.
A.
Li
,
S.
Zhou
,
L.
Qi
, and
X.
Chen
, “
A reformulated flexoelectric theory for isotropic dielectrics
,”
J. Phys. D: Appl. Phys.
48
,
465502
(
2015
).
25.
G.
Catalan
,
L.
Sinnamon
, and
J.
Gregg
, “
The effect of flexoelectricity on the dielectric properties of inhomogeneously strained ferroelectric thin films
,”
J. Condens. Matter Phys.
16
,
2253
(
2004
).
26.
E.
Eliseev
,
A.
Morozovska
,
M.
Glinchuk
, and
R.
Blinc
, “
Spontaneous flexoelectric/flexomagnetic effect in nanoferroics
,”
Phys. Rev. B
79
,
165433
(
2009
).
27.
E. A.
Eliseev
,
M. D.
Glinchuk
,
V.
Khist
,
V. V.
Skorokhod
,
R.
Blinc
, and
A. N.
Morozovska
, “
Linear magnetoelectric coupling and ferroelectricity induced by the flexomagnetic effect in ferroics
,”
Phys. Rev. B
84
,
174112
(
2011
).
28.
L.
Liu
, “
An energy formulation of continuum magneto-electro-elasticity with applications
,”
J. Mech. Phys. Solids
63
,
451
480
(
2014
).
29.
M.-M.
Yang
,
D. J.
Kim
, and
M.
Alexe
, “
Flexo-photovoltaic effect
,”
Science
360
,
904
907
(
2018
).
30.
L.
Shu
,
S.
Ke
,
L.
Fei
,
W.
Huang
,
Z.
Wang
,
J.
Gong
,
X.
Jiang
,
L.
Wang
,
F.
Li
,
S.
Lei
et al., “
Photoflexoelectric effect in halide perovskites
,”
Nat. Mater.
19
,
605
609
(
2020
).
31.
S.
Shen
and
S.
Hu
, “
A theory of flexoelectricity with surface effect for elastic dielectrics
,”
J. Mech. Phys. Solids
58
,
665
677
(
2010
).
32.
S.
Hu
and
S.
Shen
, “
Variational principles and governing equations in nano-dielectrics with the flexoelectric effect
,”
Sci. China: Phys., Mech. Astron.
53
,
1497
1504
(
2010
).
33.
P.
Yudin
and
A.
Tagantsev
, “
Fundamentals of flexoelectricity in solids
,”
Nanotechnology
24
,
432001
(
2013
).
34.
T. D.
Nguyen
,
S.
Mao
,
Y.-W.
Yeh
,
P. K.
Purohit
, and
M. C.
McAlpine
, “
Nanoscale flexoelectricity
,”
Adv. Mater.
25
,
946
974
(
2013
).
35.
P.
Zubko
,
G.
Catalan
, and
A. K.
Tagantsev
, “
Flexoelectric effect in solids
,”
Annu. Rev. Mat. Res.
24
,
387
421
(
2013
).
36.
S.
Krichen
and
P.
Sharma
, “
Flexoelectricity: A perspective on an unusual electromechanical coupling
,”
J. Appl. Mech.
83
,
030801
(
2016
).
37.
B.
Wang
,
Y.
Gu
,
S.
Zhang
, and
L.-Q.
Chen
, “
Flexoelectricity in solids: Progress, challenges, and perspectives
,”
Prog. Mater Sci.
106
,
100570
(
2019
).
38.
J.
Yvonnet
and
L.
Liu
, “
A numerical framework for modeling flexoelectricity and maxwell stress in soft dielectrics at finite strains
,”
Comput. Methods Appl. Mech. Eng.
313
,
450
482
(
2017
).
39.
T. Q.
Thai
,
T.
Rabczuk
, and
X.
Zhuang
, “
A large deformation isogeometric approach for flexoelectricity and soft materials
,”
Comput. Methods Appl. Mech. Eng.
341
,
718
739
(
2018
).
40.
A.
McBride
,
D.
Davydov
, and
P.
Steinmann
, “
Modelling the flexoelectric effect in solids: A micromorphic approach
,”
Comput. Methods Appl. Mech. Eng.
371
,
113320
(
2020
).
41.
D.
Codony
,
P.
Gupta
,
O.
Marco
, and
I.
Arias
, “
Modeling flexoelectricity in soft dielectrics at finite deformation
,”
J. Mech. Phys. Solids
146
,
104182
(
2020
).