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.

## I. INTRODUCTION

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 field^{8} (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 Kogan^{10} after the early studies by Mashkevich and Tolpygo^{11} and Tolpygo.^{12} First comprehensive theoretical works by Tagantsev^{13,14} clarified the distinction between piezoelectricity and flexoelectricity. In the mechanics community, Mindlin^{15} 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 fields^{27,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.

## II. CONTINUUM MODELING

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.

### A. Direct and converse flexoelectricity

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 $\epsilon $, the electric polarization field $P$, and their corresponding spatial gradients in the following form:

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, $\psi piezo=dijk\epsilon 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 Arias^{5} 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:

with the (effective) flexocoupling tensor

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),

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

and

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

regardless of the choice of internal energy density, i.e., with $\psi $ representing either the *original*$\psi (0)$, *direct*$\psi (Dir)$, *converse*$\psi (Con)$, or *Lifshitz-invariant*$\psi (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 $\psi (Dir)$ and converse $\psi (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 $\psi (Dir)$, the polarization gradient tensor $b$ is typically neglected, yielding

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

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 beams^{50–54} and Timoshenko beams^{55} 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.

### B. Variational principles

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, $\epsilon =12(\u2207u+(\u2207u)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*]

*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 form^{17,28}

where $\Omega $ represents the flexoelectric material with internal energy density $\psi $ (either $\psi (0)$, $\psi (Dir)$, $\psi (Con)$, or $\psi (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\epsilon 0\Vert E\Vert 2$ is the corresponding electrostatic energy density, with $\epsilon 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

subject to

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 $\Omega $ 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}

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

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\u2217$ and $P\u2217$, 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 $\Pi \xaf[u,\varphi ]$

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 $\varphi $ such that

Hence, by considering $\varphi $ 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 $\psi (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 $\psi \xaf(u,\varphi )$. Analogously to Sec. II A, one can define several enthalpy densities for flexoelectric materials, namely,

and

where $c\xaf$ is the elasticity tensor, $\kappa $ is the electric permittivity tensor, $h\xaf$ is the strain gradient elasticity tensor, $M$ is the gradient dielectricity tensor, and $\mu =\mu (1)\u2212\mu (2)$ is the flexoelectricity tensor. The direct and converse forms can be simplified as

and

where the gradient dielectricity tensor $M$ or the gradient elasticity tensor $h\xaf$ is, respectively, vanished.

The associated free enthalpy of the system reads^{52,58,59}

where $\psi \xaf$ represents one of the aforementioned free enthalpy densities. The corresponding variational principle is stated as

which is now an unconstrained optimization ($min$–$max$) 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 equations^{59}

where

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\u2217$ is easily retrieved in terms of the obtained solution ($u\u2217$, $\varphi \u2217$) as

#### 3. Equivalence between different variational principles via partial Legendre transform

The variational models based on free energy $\Pi (u,P)$ minimization and the ones based on free enthalpy $\Pi \xaf(u,\varphi )$ 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 $\psi $, the corresponding free enthalpy density $\psi \xaf$ is obtained by substracting the electrostatic energy density to the *internal* enthalpy density,^{28,60,61} $\psi \xafInt$, which is found via the following partial Legendre transform:

In the case of *simplified direct* internal energy density $\psi (Dir)(\epsilon ,\u2207\epsilon ,P)$ in Eq. (9), i.e., without gradient polarization term ($b=0$), Eq. (27) reads

The stationarity condition of the minimization results in

which is explicitly invertible to

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

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

Further, by assuming the standard expression for isotropic reciprocal dielectric susceptibility tensors $a=(\chi e\epsilon 0)\u22121I=(\epsilon \u2212\epsilon 0)\u22121I$, Eq. (32) simplifies to

which yields the standard definition of $\kappa $ for isotropic dielectrics, as a function of its electric permittivity $\epsilon $, and reveals a well-known feature of flexoelectricity: its linear growth with the dielectric susceptibility $\chi 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,M\u22600$) breaks this equivalence, since the resulting stationarity condition of the minimization, analogous to Eq. (29), is rewritten as

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:

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

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

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

hence $\u2207\xd7Elocal\u22600$, preventing writing $Elocal$ in terms of a scalar potential field $\varphi $. 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 $\varphi $ via Eq. (15), regardless whether the enthalpy density features gradient dielectricity or not.

### C. Boundary value problems

In this section, we consider the $min$–$max$ variational principle on $\Pi \xaf[u,\varphi ]$ in Sec. II B 2 and derive the corresponding boundary value problems for the Direct $\psi \xaf(Dir)(\epsilon ,\u2207\epsilon ,E)$ and Lifshitz $\psi \xaf(Lif)(\epsilon ,\u2207\epsilon ,E,\u2207E)$ free enthalpy densities in Eqs. (19) and (20).

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

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

The *jump* operator $[[\u22c5]]$ 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 $\u2202\Omega k$ [see Fig. 1(a)].

#### 1. Direct flexoelectricity

where $W\Omega $, $W\u2202\Omega $, and $WC$ represent the work densities exerted by external sources in the domain $\Omega $, on its boundary $\u2202\Omega $ and at the non-smooth region $C$ of its boundary, respectively. In this case, the admissible external sources of work are^{63}

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 $\u2202n\u2261\u2202/\u2202n$ denotes the normal derivative.

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

where $\u2202\Omega u$, $Cu$, $\u2202\Omega v$, and $\u2202\Omega \varphi $ are the Dirichlet regions where either the displacement field, its normal derivative or the electric potential are prescribed,

The Neumann regions $\u2202\Omega t$, $Cj$, $\u2202\Omega r$, and $\u2202\Omega 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 $\u2202\Omega f$ belongs to $\u2202\Omega u$, its boundaries $\u2202\u2202\Omega f$ must belong to $Cu$ as well, namely, $Cu=C\u2229\u2202\Omega u\xaf$, and hence $Cj=C\u2216Cu$.

According to Eq. (23), the equilibrium states $(u\u2217,\varphi \u2217)$ 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\u2217,\varphi \u2217)$,

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

*Find $(u,\varphi )\u2208UD\u2297PD$ such that $\u2200(\delta u,\delta \varphi )\u2208U0\u2297P0$,*

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

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

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 $\u2202\Omega v=\u2205$, which implies that only second order *Neumann* boundary conditions are allowed, i.e., $\u2202\Omega =\u2202\Omega 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 $\Omega $, as follows:

where the gradient of the variations of the displacement field $\delta ui,j$ in the last integral on $\u2202\Omega $ is expressed in terms of its tangential $\u2207jS\delta ui$ and normal $nj\u2202n\delta ui$ components, with $\u2207S$ being the surface gradient operator. This separation is required since the tangential components of $\u2207\delta u$ on $\u2202\Omega $ are dependent on $\delta u$, whereas the normal component is not. The last term in Eq. (50) is expressed in terms of $\delta u$ via integration by parts and the surface divergence theorem on $\u2202\Omega $,

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

Dirichlet boundary conditions (43) and Neumann boundary conditions

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

whereas the physical electric displacement $D$ is simply $D^$.

^{63}

*second-order geometry tensor*written in terms of the shape operator $S$ of the surface $\u2202\Omega $. The first term involves stress measures dotted with the normal vector $n$ (a first-order measure of $\u2202\Omega $, i.e., the orientation), whereas the second term involves the double stress dotted with the second-order geometry tensor (a second-order measure of $\u2202\Omega $ 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 $\mathcal{P}$ (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

Accordingly, the boundary of the domain $\u2202\Omega $ is split into several disjoint Dirichlet and Neumann boundaries as

where the high-order nature of the dielectrics leads to the definition of the Dirichlet $\u2202\Omega \phi $ and Neumann $\u2202\Omega 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\varphi $ and Neumann $C\mathcal{P}$ edge sets, respectively, where either the electric potential $\varphi $ or the electric charge density per unit length $\mathcal{P}$ are prescribed. In this case, analogously to the mechanical partition of $C$, the sets $C\varphi =C\u2229\u2202\Omega \varphi \xaf$ and $C\mathcal{P}=C\u2216C\varphi $. The corresponding Dirichlet conditions of the problem are mathematically written as

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

The variational principle in Eq. (23) implies that feasible equilibrium states $u\u2217\u2208UD$, $\varphi \u2217\u2208PD$ must meet Dirichlet boundary conditions (59), and their corresponding admissible variations $\delta u\u2208U0$, $\delta \varphi \u2208P0$ must vanish as stated in Eq. (45), where

and

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

*Find $(u,\varphi )\u2208UD\u2297PD$ such that, $\u2200(\delta u,\delta \varphi )\u2208U0\u2297P0$,*

The Cauchy stress $\sigma ^$, the double stress $\sigma ~$, the electric displacement $D^$, and the double electric displacement $D~$ in Eq. (63) are the conjugate quantities to the strain $\epsilon $, the strain gradient $\u2207\epsilon $, the electric field $E$, and the electric field gradient $\u2207E$, respectively, as follows:

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

complemented with Dirichlet conditions (59) and Neumann conditions

The physical stress $\sigma $ and physical electric displacement $D$ arise from Eq. (65) as

#### 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 $\u2202\Omega $ 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\u2217,\varphi \u2217)$ 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 method^{68} 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 freedom^{69}) as compared to other alternatives such as the Lagrange multipliers or penalty methods.

Following Nitsche’s approach, an additional term $\Pi \xafNitsche[u,\varphi ]$ 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:

with the numerical parameters $\beta u,\beta v,\beta Cu,\beta \varphi ,\beta \phi ,\beta C\varphi \u2208R+$ 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\u2217,\varphi \u2217)$ being, respectively, actual minima and maxima of the energy functional with respect to $u$ and $\varphi $.^{63}

*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 $\zeta \u2208R+$ as follows:

^{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,

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(\Omega )]D$ are typically considered.^{63} Analogously, the Lifshitz-invariant form requires

Many authors in computational flexoelectricity neglect the edge Dirichlet conditions (43a) on $Cu$^{9,57,58,71} [or analogously (59a) on $C\varphi $]. 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 $\u2202\Omega u$ ($\u2202\Omega \varphi $) automatically implies the strong imposition on the adjacent edges in $Cu$ ($C\varphi $) 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 $\varphi $ [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 $\u2202\Omega \Phi i$ and its edge $\u2202\u2202\Omega \Phi i$ is mathematically written as

where $\Phi 1,\u2026,\Phi N\Phi \u2208R$ are additional state variables and $N\Phi $ denotes the number of sensing electrodes in the system. Since each $\Phi i$ is not prescribed on $\u2202\Omega \Phi i$ nor $\u2202\u2202\Omega \Phi i$, their energy conjugates—i.e., the total surface charge density on the electrode^{63} and on their edges—must vanish, that is,

In the direct flexoelectricity model, high-order electrostatics are not present, hence Eqs. (72) and (73) must hold on $\u2202\Omega \Phi 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 $\varphi $ to all the nodes on the electrode, and Eq. (73) is then automatically fulfilled. Otherwise, they must be enforced weakly, e.g., via Lagrange multipliers^{72} 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 stacks^{46,75,76} or structured materials with geometrically polarized cavities^{2,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).

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

and high-order equilibrium conditions read

where $L$ and $R$ represent either sides across the interface $I$, $i\u2260j$ 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.

##### 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 $\Omega \u2208(0,L1)\xd7(0,L2)\xd7(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\u2032$ (i.e., $x1=L1$, $x2=L2$, or $x3=L3$), as well as (ii) between the intersections of $A$ and $A\u2032$ with the actual boundary of the infinite lattice, denoted as $S$ and $S\u2032$, respectively (see Fig. 4).

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

where $[[\u22c5]]A$ are either known (enforced) or unknown (resolved) constant quantities denoting the difference of the state variable within the brackets $(\u22c5)$ 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}

### D. Flexoelectricity in soft materials: Beyond infinitesimal deformations

The interest on computational flexoelectricity in soft materials has increased in recent years^{22,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 $\Omega 0$ in the reference (or undeformed) configuration and by $\Omega $ in the current (or deformed) configuration. The deformation map $\chi :\Omega 0\u2192\Omega $ maps every material point $X\u2208\Omega 0$ to the spatial point $x=\chi (X)\u2208\Omega $. Relevant deformation measures in this context are the deformation gradient $F=\u22070\chi $, the Jacobian determinant $J=det(F)$, and the right Cauchy–Green deformation tensor $C=FT\u22c5F$. In the reference configuration, the strain is represented by the Green–Lagrangian strain tensor $E=(C\u2212I)/2$, and the electric field in by $E=\u2212\u22070\Phi $.

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

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

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), as^{41}

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 ${\chi \u2217,\Phi \u2217}$ are located at the saddle points of $\Pi \xaf$ satisfying

where $X$ and $P$ are the functional spaces for $\chi $ and $\Phi $ with sufficient regularity fulfilling Dirichlet boundary conditions. Vanishing the first variations of $\Pi \xaf[\chi ,\Phi ]$ leads to the weak form of the boundary value problem, whose corresponding Euler–Lagrange equations read as follows:^{41}

The physical second Piola–Kirchhoff stress $S$ and the electric displacement $D$ are defined in terms of $\chi $ and $\Phi $ as

where

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

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}

## III. NUMERICAL MODELING

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 rod^{41} 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.

### A. Uniform open B-spline basis

B-spline functions^{97–99} are smooth, positive-valued piece-wise polynomials with compact support. Being $p$ the polynomial degree, they are by construction $Cp\u22121$-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 $\xi \u2208\Xi =[k0,km]$ in terms of the knot vector $k=[k0,k1,k2,\u2026,km]$, where $ki<ki+1$. The $i$th function of this basis is defined recursively as^{97}

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 $\xi =k0$ to $\xi =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 $i$th B-spline function can be expressed as a translation of the first one as $Bip(\xi )=B1p(\xi \u2212i+1)$ (see Fig. 5).

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,\u2026,k0,k1,k2,\u2026,km,\u2026,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 $C\u22121$ (i.e., discontinuous), yielding a boundary-interpolant (or *open*) basis suitable to enforce essential boundary conditions strongly (see Fig. 6).

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

which is defined on the $D$-dimensional parametric space $\xi \u2208\Xi =[k1,0,k1,m1]\u2297\u2026\u2297[kd,0,kd,md]$.

In order to use B-spline basis (defined in the parametric space $\Xi $) to approximate the state variables $u$ and $\varphi $ (defined in the physical space $\Omega $), let us define the geometrical map

which maps a given point $\xi \u2208\Xi $ in the parametric space to a given point $x\u2208\Omega $ in the physical space. The basis functions $N(x)\u2208\Omega $ are defined as $N=[Bp\xb0\phi \u22121]$, in such a way that

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

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)\u2208\Omega $, since the tensor product structure of the basis is preserved, e.g.,

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\varphi $, and hence, the approximated fields $uh(x)$ and $\varphi h(x)$ are obtained.

### B. Immersed boundary B-spline approximation

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 $\Omega $ is embedded into a larger domain $\Omega \u25a1$ with rectangular/cuboidal shape, i.e., $\Omega \u2286\Omega \u25a1$. The geometrical map in Eq. (87) is redefined as

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

Since the approximation space is unfitted to the geometry, it is not interpolant on $\u2202\Omega $, 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 $\u2202\Omega $ is allowed to intersect the cells $\Omega \u25a1c$ of the embedding mesh arbitrarily, leading to a unfitted discretization of $\Omega $. The cells with a nonempty intersection with $\Omega $ are labeled as active, whereas outer (or inactive) cells are neglected. Active cells can be either inner [$\Omega \u25a1c\u2286\Omega $, see Fig. 7(b)] or cut [$\Omega \u25a1c\u2288\Omega $, see Fig. 7(c)].

#### 4. Cell classification

Labeling a cell depending on its intersection with $\Omega $ is usually accomplished by checking whether all vertices of each cell (and possibly more points within the cell) lie within the domain $\Omega $ (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 $\Omega $, it is not interpolant on its boundary $\u2202\Omega $, hence, Dirichlet boundary conditions must be weakly enforced, e.g., by means of Nitsche’s method (cf. Sec. II C 3).

#### 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 $\Omega \u25a1c$ since only their intersection with $\Omega $, i.e., $\Omega \u25a1c\u2229\Omega $, is to be considered. The standard approach consists on dividing $\Omega \u25a1c\u2229\Omega $ 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 schemes^{110,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 $\Omega $ are considered.

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 $\u2202\Omega $ is considered. In the case of implicit boundary representations such as level sets, an explicit approximation of $\u2202\Omega $ 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 $\Omega $. 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 $\eta c=|\Omega \u25a1c\u2229\Omega |/|\Omega \u25a1c|\u226a1$. In particular, for the boundary value problems of flexoelectricity considered in this paper, the condition number scales with the minimum volume fraction $\eta min=minc(\eta c)$ in meshes of fixed size $h$ at a rate of $\eta min\u2212(2p+1\u22122/D)$,^{62} which implies that ill-conditioning is more severe for high-order basis (see Fig. 10).

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 $\Omega $) 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 $\eta 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}

### C. Beyond uniform meshes: Local *h*-refinement via hierarchical B-splines

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.

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 relation*^{119}

where $\u2113\u2208N$ stands for the refinement level, $\xi \u2113+1=2\xi \u2113$ and the scaling coefficients

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

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.

### D. Nondimensionalization of the equations

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\varphi \varphi $ 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 $\epsilon $ of the material. In a typical dielectric, they are in the order of $109$ (or higher) and $10\u22129$ (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

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 $\varphi 0=\varphi /(fLfSfP\u22121)$,

Young’s modulus $Y0=Y/fS$,

the flexoelectric coefficients $\mu 0=\mu /(fPfL)$,

the electric permittivity $\epsilon 0=\epsilon /(fP2fS\u22121)$,

and so on. Typically, suitable values for the normalizing factors are $fL\u223cLgeom$, $fS\u223cY$, and $fP\u223cY\epsilon $, where $Lgeom$ is a characteristic length of the geometry.

The normalized flexoelectric coefficients $\mu 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.

## IV. COMPUTATIONAL EXPERIMENTS OF CANTILEVER BEAM BENDING WITH DIFFERENT FORMULATIONS

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 numerically^{58,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.

### A. Cantilever beam bending

We consider a 2D (plane strain) cantilever beam of length $L=8\mu m$ and thickness $H=0.4\mu m$. The material properties are simple enough to isolate the transversal flexoelectric effect, i.e., a Young modulus $Y=100$ GPa, electric permittivity $\epsilon =11$ nC/V m, and transversal flexoelectric coefficient $\mu T=1\mu 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=\u22121\mu N/\mu 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.

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 beam^{41} is different in each case. Comparing the maximum deflection of the direct-flexoelectric beam ($0.30\mu m$) and the Lifshitz-invariant-flexoelectric beam ($0.24\mu m$) with respect to a standard elastic one ($0.32\mu 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.

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.

### B. Cantilever beam actuator

In this experiment, we explore the transversal flexoelectric effect triggered by electrical actuation. This device was first used by Bursian and Zaikovskii^{3} 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.

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\mu m$ for the direct case and only $0.12\mu 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)].

## V. COMPUTATIONAL ANALYSIS OF A SOFT CANTILEVER ROD ACTUATOR

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 Zaikovskii^{3} 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\mu 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 $\Phi $ at the top boundary is set to a certain value $V$, and the bottom boundary is grounded ($\Phi =0$). The voltage drop induces a transverse electric field $EY=\u2212V/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.

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:

with Young’s modulus $Y=1.0$ GPa, Poisson’s ratio $\nu =0.37$ (that is, Lamé’s moduli $\lambda \u22481.039$ GPa and $\mu =0.685$ GPa), mechanical length scale $\u2113=30$ nm, dielectric permittivity $\epsilon =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 $\mu 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 $\mu T$ and negative $\mu 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)].

## VI. COMPUTATIONAL ANALYSIS OF FLEXOELECTRIC DEVICES AND METAMATERIALS

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.

### A. Collective-beam bending scalable flexoelectric device

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 $\u2113=400$ nm, and beam spacings of $h=280$ nm, under three point displacement as shown in Fig. 17(a), with $\delta =40$ nm. The direct formulation for flexoelectricity in Sec. II C 1 is used, with the interface conditions in Sec. II C 4 2 at the dielectric-insulator interfaces. The material parameters are reported in Table I.

Material | Y | ν | ℓ_{mech} | ɛ | μ_{L} | μ_{T} | μ_{S} |

(GPa) | — | (nm) | (nC/V m) | (μC/m) | (μC/m) | (μC/m) | |

Matrix | 152 | 0.33 | 1 | 141 | 150 | 110 | 110 |

Insulators | 152 | 0.33 | 1 | 0 | 0 | 0 | 0 |

Material | Y | ν | ℓ_{mech} | ɛ | μ_{L} | μ_{T} | μ_{S} |

(GPa) | — | (nm) | (nC/V m) | (μC/m) | (μC/m) | (μC/m) | |

Matrix | 152 | 0.33 | 1 | 141 | 150 | 110 | 110 |

Insulators | 152 | 0.33 | 1 | 0 | 0 | 0 | 0 |

Figure 17(b) shows the electric potential, $\varphi $, 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.

### B. Geometrically-polarized periodic inclusions

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}

Y | ν | ℓ_{mech} | ɛ | ℓ_{elec} | μ_{L} | μ_{T} | μ_{S} | |

(GPa) | — | (nm) | (nC/V m) | (nm) | (μC/m) | (μC/m) | (μC/m) | |

152 | 0.33 | 50 | 8 | 300 | 1.2 | 1.1 | 0.05 |

Y | ν | ℓ_{mech} | ɛ | ℓ_{elec} | μ_{L} | μ_{T} | μ_{S} | |

(GPa) | — | (nm) | (nC/V m) | (nm) | (μC/m) | (μC/m) | (μC/m) | |

152 | 0.33 | 50 | 8 | 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=[[\varphi ]]x1=0$), whereas the jump in electric potential along the vertical direction ($[[\varphi ]]x2$) is resolved as a result of a prescribed vertical jump in the displacement field ($[[u]]x2=[[u]]x2\xaf$).

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\mu $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}

### C. Bending-dominated lattice metamaterial

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 ($[[\varphi ]]x2=[[\varphi ]]x2\xaf$), 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\mu $m and $L2=2.8\mu $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.

## VII. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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

## AUTHOR DECLARATIONS

### Conflicts of Interest

The authors have no conflicts to disclose.

### APPENDIX A: MATERIAL CHARACTERIZATION

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

#### a. Elasticity tensor **c**

Isotropic elasticity is represented by the fourth-order tensor $c$, which depends on $\lambda $ and $\mu $ in the following form:

The Lamé parameters $\lambda $ and $\mu $ are related to the Young modulus $Y$ and Poisson’s ratio $\nu $ as $\lambda =Y\nu /(1+\nu )(1\u22122\nu )$ and $\mu =Y/2(1+\nu )$.

#### b. Strain gradient elasticity tensor **h**

We consider an isotropic simplified strain gradient elasticity tensor,^{132} which depends on $\lambda $, $\mu $ and the mechanical length scale $\u2113mech$ in the following form:

#### c. Dielectricity tensor *κ*

*κ*Isotropic dielectricity is represented by the second-order tensor $\kappa $, which depends on the electric permittivity $\epsilon $ as

#### 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 $\epsilon $ and the dielectric length scale $\u2113elec$ as

#### e. Flexoelectricity tensor *μ*

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

The flexoelectric tensor $\mu $ oriented in an arbitrarily rotated orthonormal basis is obtained by rotating $\mu \u27e8x1\u27e9$ as

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 $\mu L=\mu T+2\mu S$.

### APPENDIX B: RESIDUAL-BASED WEAK FORM STABILIZATION

As discussed in Remark 8, depending on the geometry of $\Omega $, 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\mu m$ long by $H=0.4\mu m$ thick. The dielectric material has Young’s modulus $Y=100$ GPa, Poisson’s ratio $\nu =0.37$, dielectric permittivity $\epsilon =11nJ/V2m$, and longitudinal and transversal flexoelectric coefficients $\mu L=\mu T=1\mu J/Vm$. Other material coefficients are $0$. A force of $F=100\mu N$ is applied vertically on the top of the right end.

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.

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)]

depending on just five material parameters, namely, the Young modulus $Y$ and mechanical length scale $\u2113mech$, the dielectric constant $\epsilon $ and the dielectric length scale $\u2113elec$, and the flexoelectric coefficient $\mu $.

The corresponding weak form can be written as

with appropriate boundary conditions, where

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

with the residuals

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

The stabilization parameters $\tau mech,\tau 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 Codony^{62} as

where $h$ represents the mesh size and $\alpha mech,\alpha 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 $\alpha mech$ and $\alpha 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,\varphi}$ for large enough stabilization parameters.

## REFERENCES

*et al.*, “