Mathematical and computational modeling of flexoelectricity

Codony, D.; Mocci, A.; Barceló-Mercader, J.; Arias, I.

doi:10.1063/5.0067852

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,⁷ under inhomogeneous applied electric field⁸ (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.⁹

Phenomenological models for flexoelectricity in crystalline dielectrics were first proposed by Kogan¹⁰ after the early studies by Mashkevich and Tolpygo¹¹ and Tolpygo.¹² First comprehensive theoretical works by Tagantsev^13,14 clarified the distinction between piezoelectricity and flexoelectricity. In the mechanics community, Mindlin¹⁵ 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.¹⁶ More recently, a simplified framework for isotropic dielectrics was proposed by Maranganti et al.¹⁷

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.³¹ General variational principles for flexoelectric materials can be found in Liu,²⁸ Shen and Hu,³¹ Hu and Shen.³² The reader is referred to Yudin and Tagantsev,³³ Nguyen et al.,³⁴ Zubko et al.,³⁵ Krichen and Sharma,³⁶ Wang et al.³⁷ 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 $R^{3}$ ⁠, 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.¹⁷ Yudin and Tagantsev,³³ 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:

\begin{aligned} ψ^{(0)} (ε, \nabla ε, P, \nabla P) & = \frac{1}{2} c_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} h_{i j k l m n} ε_{i j, k} ε_{l m, n} + \frac{1}{2} a_{k l} P_{k} P_{l} \\ + \frac{1}{2} b_{i j k l} P_{i, k} P_{j, l} - f_{l i j k}^{(1)} ε_{i j, k} P_{l} - f_{l i j k}^{(2)} ε_{i j} P_{l, k}, \end{aligned}

(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.¹⁵ 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} = d_{i j k} ε_{i j} P_{k}$ ⁠, 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⁵ and Bhaskar et al.⁴²). 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. can be rewritten as follows:

\begin{aligned} - f_{l i j k}^{(1)} ε_{i j, k} P_{l} - f_{l i j k}^{(2)} ε_{i j} P_{l, k} \\ = - f_{l i j k} ε_{i j, k} P_{l} - f_{l i j k}^{(2)} \frac{\partial (ε_{i j} P_{l})}{\partial x_{k}} \end{aligned}

(2a)

= f_{l i j k} ε_{i j} P_{l, k} - f_{l i j k}^{(1)} \frac{\partial (ε_{i j} P_{l})}{\partial x_{k}}

(2b)

= - \frac{1}{2} f_{l i j k} (ε_{i j, k} P_{l} - ε_{i j} P_{l, k}) - \frac{1}{2} (f_{l i j k}^{(1)} + f_{l i j k}^{(2)}) \frac{\partial (ε_{i j} P_{l})}{\partial x_{k}},

(2c)

with the (effective) flexocoupling tensor

f_{l i j k} = f_{l i j k}^{(1)} - f_{l i j k}^{(2)} .

(3)

The first terms in – 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 – are null-Lagrangians,⁴⁸ in the sense that their bulk integral can be written as a surface integral by means of the divergence theorem, e.g., for ,

\begin{aligned} \int_{Ω} \frac{1}{2} (f_{l i j k}^{(1)} + f_{l i j k}^{(2)}) \frac{\partial (ε_{i j} P_{l})}{\partial x_{k}} d Ω \\ = \int_{Γ} \frac{1}{2} (f_{l i j k}^{(1)} + f_{l i j k}^{(2)}) ε_{i j} P_{l} n_{k} d Γ . \end{aligned}

(4)

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

\begin{aligned} ψ^{(Dir)} (ε, \nabla ε, P, \nabla P) & = \frac{1}{2} c_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} h_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ + \frac{1}{2} a_{k l} P_{k} P_{l} + \frac{1}{2} b_{i j k l} P_{i, j} P_{k, l} - f_{l i j k} ε_{i j, k} P_{l}, \end{aligned}

(5)

\begin{aligned} ψ^{(Con)} (ε, \nabla ε, P, \nabla P) & = \frac{1}{2} c_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} h_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ + \frac{1}{2} a_{k l} P_{k} P_{l} + \frac{1}{2} b_{i j k l} P_{i, j} P_{k, l} + f_{l i j k} ε_{i j} P_{l, k}, \end{aligned}

(6)

and

\begin{aligned} ψ^{(Lif)} (ε, \nabla ε, P, \nabla P) & = \frac{1}{2} c_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} h_{i j k l m n} ε_{i j, k} ε_{l m, n} + \frac{1}{2} a_{k l} P_{k} P_{l} \\ + \frac{1}{2} b_{i j k l} P_{i, j} P_{k, l} - \frac{1}{2} f_{l i j k} (ε_{i j, k} P_{l} - ε_{i j} P_{l, k}) . \end{aligned}

(7)

As derived later in Eq. , the bulk constitutive equations for the physical stress

σ

and physical electric field

E

are

E_{l} = \frac{\partial ψ}{\partial P_{l}} - \frac{\partial}{\partial x_{k}} (\frac{\partial ψ}{\partial P_{l, k}}) = a_{k l} P_{k} - b_{i j l k} P_{i, j k} - f_{l i j k} ε_{i j, k},

(8a)

σ_{i j} = \frac{\partial ψ}{\partial ε_{i j}} - \frac{\partial}{\partial x_{k}} (\frac{\partial ψ}{\partial ε_{i j, k}}) = c_{i j k l} ε_{k l} - h_{i j k l m n} ε_{l m, n k} + f_{l i j k} P_{l, 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.⁴⁹ 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)} (ε, \nabla ε, P) = {ψ^{(Dir)} (ε, \nabla ε, P, \nabla P) |}_{b = 0} .

(9)

In turn, the strain gradient elasticity tensor

h

can be neglected in

ψ^{(Con)}

⁠, resulting in

ψ^{(Con)} (ε, P, \nabla P) = {ψ^{(Con)} (ε, \nabla ε, P, \nabla 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 beams^50–54 and Timoshenko beams⁵⁵ 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, $ε = \frac{1}{2} (\nabla u + {(\nabla 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.²⁸ Such free energy takes the form^17,28

Π [u, P] = \int_{Ω} (ψ (u, P) + \frac{1}{2} ε_{0} {‖ E ‖}^{2}) d Ω - W^{ext},

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

\frac{1}{2} ε_{0} {‖ E ‖}^{2}

is the corresponding electrostatic energy density, with

ε_{0}

being the vacuum permittivity, and

W^{ext}

is the external work.

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

(u^{*}, P^{*}) = \arg min_{u} min_{P} Π [u, P],

(12a)

subject to

\nabla \times E = 0 in Ω,

(12b)

\nabla \cdot D = q in Ω,

(12c)

where Eqs. and 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. and 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

\nabla \cdot σ = - f^{ext} in Ω,

(13a)

\nabla \cdot (ε_{0} E + P) = q in Ω,

(13b)

\nabla \times E = 0 in Ω,

(13c)

where

f^{ext}

represents the external body forces per unit volume, and

σ (u, P) = \frac{d ψ (u, P)}{d ε} = [\frac{\partial ψ}{\partial ε} - \nabla \cdot \frac{\partial ψ}{\partial \nabla ε}] (u, P),

(14a)

E (u, P) = \frac{d ψ (u, P)}{d P} = [\frac{\partial ψ}{\partial P} - \nabla \cdot \frac{\partial ψ}{\partial \nabla P}] (u, P) .

(14b)

The system of equations in is composed by the differential Eqs. and , corresponding to a (fourth-order) coupled elliptic problem and by Eq. which represents an additional constraint, requiring the irrotationality of the resulting electric field

E

⁠. This constraint hinders the solution of the system in Eq. to find

u^{*}

and

P^{*}

⁠, either by means of analytical closed forms or numerical approximations.²⁸ 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 $\bar{Π} [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. implies the existence of an electric potential

ϕ

such that

E = - \nabla ϕ .

(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

\bar{ψ} (u, ϕ)

⁠. Analogously to Sec. II A, one can define several enthalpy densities for flexoelectric materials, namely,

\begin{aligned} {\bar{ψ}}^{(0)} (ε, \nabla ε, E, \nabla E) = \frac{1}{2} {\bar{c}}_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} {\bar{h}}_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ - \frac{1}{2} κ_{k l} E_{k} E_{l} - \frac{1}{2} M_{i j k l} E_{i, j} E_{k, l} - μ_{l i j k}^{(1)} ε_{i j, k} E_{l} - μ_{l i j k}^{(2)} ε_{i j} E_{l, k}, \end{aligned}

(16)

\begin{aligned} {\bar{ψ}}^{(Dir)} (ε, \nabla ε, E, \nabla E) = \frac{1}{2} {\bar{c}}_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} {\bar{h}}_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ - \frac{1}{2} κ_{k l} E_{k} E_{l} - \frac{1}{2} M_{i j k l} E_{i, j} E_{k, l} - μ_{l i j k} ε_{i j, k} E_{l}, \end{aligned}

(17)

\begin{aligned} {\bar{ψ}}^{(Con)} (ε, \nabla ε, E, \nabla E) = \frac{1}{2} {\bar{c}}_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} {\bar{h}}_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ - \frac{1}{2} κ_{k l} E_{k} E_{l} - \frac{1}{2} M_{i j k l} E_{i, j} E_{k, l} + μ_{l i j k} ε_{i j} E_{l, k}, \end{aligned}

(18)

and

\begin{aligned} {\bar{ψ}}^{(Lif)} (ε, \nabla ε, E, \nabla E) = \frac{1}{2} {\bar{c}}_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} {\bar{h}}_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ - \frac{1}{2} κ_{k l} E_{k} E_{l} - \frac{1}{2} M_{i j k l} E_{i, j} E_{k, l} - \frac{1}{2} μ_{l i j k} (ε_{i j, k} E_{l} - ε_{i j} E_{l, k}), \end{aligned}

(19)

where

\bar{c}

is the elasticity tensor,

κ

is the electric permittivity tensor,

\bar{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

{\bar{ψ}}^{(Dir)} (ε, \nabla ε, E) = {{\bar{ψ}}^{(Dir)} (ε, \nabla ε, E, \nabla E) |}_{M = 0}

(20)

and

{\bar{ψ}}^{(Con)} (ε, E, \nabla E) = {{\bar{ψ}}^{(Con)} (ε, \nabla ε, E, \nabla E) |}_{\bar{h} = 0},

(21)

where the gradient dielectricity tensor

M

or the gradient elasticity tensor

\bar{h}

is, respectively, vanished.

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

\bar{Π} [u, ϕ] = \int_{Ω} \bar{ψ} (u, ϕ) d Ω - W^{ext},

(22)

where

\bar{ψ}

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

(u^{*}, ϕ^{*}) = \arg min_{u} max_{ϕ} \bar{Π} [u, ϕ],

(23)

which is now an unconstrained optimization (⁠

min

–

max

⁠) problem given that Eq. is automatically fulfilled and that Eq. follows directly from the solution of the variational principle as one of the two associated Euler–Lagrange equations⁵⁹

\nabla \cdot σ = - f^{ext} in Ω,

(24a)

\nabla \cdot D = q in Ω,

(24b)

where

σ (u, ϕ) = \frac{d \bar{ψ} (u, ϕ)}{d ε} = [\frac{\partial \bar{ψ}}{\partial ε} - \nabla \cdot \frac{\partial \bar{ψ}}{\partial \nabla ε}] (u, ϕ),

(25a)

D (u, ϕ) = - \frac{d \bar{ψ} (u, ϕ)}{d E} = [- \frac{\partial \bar{ψ}}{\partial E} + \nabla \cdot \frac{\partial \bar{ψ}}{\partial \nabla E}] (u, ϕ) .

(25b)

Equations and describe a (fourth-order) coupled elliptic problem equivalent to that in Eqs. and , 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} \nabla ϕ^{*} .

(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 $\bar{Π} (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

\bar{ψ}

is obtained by substracting the electrostatic energy density to the internal enthalpy density,^28,60,61

{\bar{ψ}}^{Int}

⁠, which is found via the following partial Legendre transform:

\begin{aligned} \bar{ψ} (u, ϕ) & = {\bar{ψ}}^{Int} (u, ϕ) - \frac{1}{2} ε_{0} {‖ E ‖}^{2} \\ = min_{P} (ψ (u, P) - P \cdot E) - \frac{1}{2} ε_{0} {‖ E ‖}^{2} . \end{aligned}

(27)

In the case of simplified direct internal energy density

ψ^{(Dir)} (ε, \nabla ε, P)

in Eq. , i.e., without gradient polarization term (⁠

b = 0

⁠), Eq. reads

\begin{aligned} {\bar{ψ}}^{(Dir)} (u, ϕ) & = min_{P} (\frac{1}{2} c_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} h_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ + \frac{1}{2} a_{m l} P_{m} P_{l} - f_{l i j k} ε_{i j, k} P_{l} - P_{l} E_{l}) - \frac{1}{2} ε_{0} {‖ E ‖}^{2} . \end{aligned}

(28)

The stationarity condition of the minimization results in

E_{l} = \frac{d ψ}{d P_{l}} = a_{m l} P_{m} - f_{l i j k} ε_{i j, k},

(29)

which is explicitly invertible to

P_{m} = a_{l m}^{- 1} (E_{l} + f_{l i j k} ε_{i j, k}) = a_{l m}^{- 1} (E_{l} + E_{l}^{Flexo}),

(30)

where

E_{l}^{Flexo} = f_{l i j k} ε_{i j, k}

is well-known as the flexoelectric field. By inserting Eq. into , the free enthalpy density is obtained as

\begin{aligned} {\bar{ψ}}^{(Dir)} (u, ϕ) & = \frac{1}{2} c_{i j k l} ε_{i j} ε_{k l} - \frac{1}{2} (a_{i j}^{- 1} + ε_{0} δ_{i j}) E_{i} E_{j} \\ + \frac{1}{2} (h_{i j k l m n} - a_{l m}^{- 1} f_{m i j k} f_{b l m n}) ε_{i j, k} ε_{l m, n} \\ - (a_{l m}^{- 1} f_{m i j k}) E_{l} ε_{i j, k} . \end{aligned}

(31)

The relation between the tensors

c, a, f, h

in the direct internal energy density and their counterparts

\bar{c}, κ, μ, \bar{h}

in the free enthalpy density is revealed by Eq. as follows:

{\bar{c}}_{i j k l} = c_{i j k l},

(32a)

κ_{i j} = a_{i j}^{- 1} + ε_{0} δ_{i j},

(32b)

μ_{l i j k} = a_{l m}^{- 1} f_{m i j k},

(32c)

{\bar{h}}_{i j k l m n} = h_{i j k l m n} - a_{a b}^{- 1} f_{a i j k} f_{b l m n} .

(32d)

Further, by assuming the standard expression for isotropic reciprocal dielectric susceptibility tensors

a = {(χ_{e} ε_{0})}^{- 1} I = {(ε - ε_{0})}^{- 1} I

⁠, Eq. simplifies to

{\bar{c}}_{i j k l} = c_{i j k l},

(33a)

κ_{i j} = ε δ_{i j},

(33b)

μ_{l i j k} = χ_{e} ε_{0} f_{l i j k},

(33c)

{\bar{h}}_{i j k l m n} = h_{i j k l m n} - χ_{e} ε_{0} f_{a i j k} f_{a l m n},

(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 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. , hold also for other density forms without gradient polarization/dielectricity terms (⁠

b, M = 0

⁠), differing from and in null Lagrangians only.⁶² However, considering the gradient polarization/dielectricity terms (⁠

b, M \neq 0

⁠) breaks this equivalence, since the resulting stationarity condition of the minimization, analogous to Eq. , is rewritten as

E_{l} = \frac{d ψ}{d P_{l}} = a_{m l} P_{m} - f_{l i j k} ε_{i j, k} - b_{i j l k} P_{i, j k},

(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

E^{local}

and nonlocal

E^{nl}

fields as follows:

E = E^{local} + E^{nl},

(35a)

E_{l}^{local} = a_{m l} P_{m} - f_{l i j k} ε_{i j, k},

(35b)

E_{l}^{nl} = - b_{i j l k} P_{i, j k}

(35c)

and invert Eq. , analogously to Eq. ,

P_{m} = a_{l m}^{- 1} (E_{l}^{local} + f_{l i j k} ε_{i j, k}) .

(36)

With these definitions, both

E (ε, E^{local})

and

P (ε, E^{local})

can be expressed in terms of

ε

and

E^{local}

⁠, and introduced in Eq. to yield

\begin{aligned} {\bar{ψ}}^{(Dir)} (u, E^{local}) = \frac{1}{2} {\bar{c}}_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} {\bar{h}}_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ - \frac{1}{2} κ_{k l} E_{k}^{local} E_{l}^{local} - μ_{l i j k} ε_{i j, k} E_{l}^{local} \\ - \frac{1}{2} (a_{i p}^{- 1} a_{k q}^{- 1} b_{p j q l}) E_{i, j}^{local} E_{k, l}^{local}, \end{aligned}

(37)

where higher order terms—i.e., terms depending on

\nabla \nabla ε

⁠,

\nabla \nabla E^{local}

or higher derivatives—have been neglected. Upon inspection of Eqs. and , the following relation is revealed:

M_{i j k l} = a_{i p}^{- 1} a_{k q}^{- 1} b_{p j q l} .

(38)

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

\nabla \times E = \nabla \times E^{local} + \nabla \times E^{nl} = 0,

(39)

hence

\nabla \times E^{local} \neq 0

⁠, preventing writing

E^{local}

in terms of a scalar potential field

ϕ

⁠. Therefore, the associated variational formulation remains constrained by Maxwell–Faraday’s law in Eq. . This fact is sometimes overlooked in the literature, where there is no distinction made between

E^{local}

and

E

⁠, expressed in terms of

ϕ

via Eq. , 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 $\bar{Π} [u, ϕ]$ in Sec. II B 2 and derive the corresponding boundary value problems for the Direct ${\bar{ψ}}^{(Dir)} (ε, \nabla ε, E)$ and Lifshitz ${\bar{ψ}}^{(Lif)} (ε, \nabla ε, E, \nabla E)$ free enthalpy densities in Eqs. (19) and (20).

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

FIG. 1.

Sketch of the geometry of Ω ∈ R 3. (a) Detail of ∂ Ω subdivided in smooth portions ∂ Ω i and ∂ Ω j, with their corresponding normal vectors n i and n j, (b) detail of ∂ Ω i, with the triplet { m i , s i , n i } defined on ∂ ∂ Ω i, and (c) detail of ∂ Ω j, with the triplet { m j , s j , n j } 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.

View large Download slide

Sketch of the geometry of $Ω \in R^{3}$ ⁠. (a) Detail of $\partial Ω$ subdivided in smooth portions $\partial Ω_{i}$ and $\partial Ω_{j}$ ⁠, with their corresponding normal vectors $n^{i}$ and $n^{j}$ ⁠, (b) detail of $\partial Ω_{i}$ ⁠, with the triplet ${m^{i}, s^{i}, n^{i}}$ defined on $\partial \partial Ω_{i}$ ⁠, and (c) detail of $\partial Ω_{j}$ ⁠, with the triplet ${m^{j}, s^{j}, n^{j}}$ defined on $\partial \partial Ω_{j}$ ⁠. The orientation of the tangent vector $s$ is arbitrary and not relevant in this context. This figure was published in Codony et al.⁶³ Copyright 2019 Elsevier.

FIG. 2.

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

View large Download slide

Sketch of the geometry of $Ω \in R^{2}$ ⁠. (a) $Ω$ and the boundary $\partial Ω$ subdivided in several smooth portions $\partial Ω_{f}$ ⁠, e.g., $\partial Ω_{i}$ and $\partial Ω_{j}$ with their corresponding normal vectors $n^{i}$ and $n^{j}$ ⁠. The intersection between $\partial Ω_{i}$ and $\partial Ω_{j}$ is the corner $C^{i j}$ ⁠. (b) Detail of $\partial Ω_{i}$ ⁠, with the pair ${m^{i}, n^{i}$ } defined on $C^{i j}$ ⁠, and (c) detail of $\partial Ω_{j}$ ⁠, with the pair ${m^{j}, n^{j}}$ defined on $C^{i j}$ ⁠. This figure was published in Codony et al.⁶³ Copyright 2019 Elsevier.

The boundary of the $f$ th portion of $\partial Ω$ is a closed curve (⁠ $D = 3$ ⁠) or a pair of points (⁠ $D = 2$ ⁠) denoted as $\partial \partial Ω_{f}$ ⁠. The union $C = ⋃_{f} \partial \partial Ω_{f}$ represents the regions at which $\partial Ω$ is not smooth, and $C^{i j} = \partial \partial Ω_{i} \cap \partial \partial Ω_{j}$ ⁠, $(i \neq j)$ ⁠, represents the part of $C$ that is adjacent to $\partial Ω_{i}$ and $\partial Ω_{j}$ ⁠. In case $D = 3$ ⁠, at each point $x \in \partial \partial Ω_{f}$ ⁠, we define $m^{f}$ as the unit co-normal vector pointing outward of $\partial Ω_{f}$ ⁠, which is orthogonal to the normal vector $n^{f}$ and to the tangent vector $s^{f}$ of the curve $\partial \partial Ω_{f}$ [see Figs. 1(b) and 1(c)]. In case $D = 2$ ⁠, $m^{f}$ is defined on $\partial \partial Ω_{f}$ as the outward unit tangent vector to $\partial Ω_{f}$ [see Figs. 2(b) and 2(c)].

The jump operator $[[\cdot]]$ acting on a given quantity $a$ is defined on $C^{i j}$ as the sum of that quantity evaluated at both sides of $C^{i j}$ ⁠, namely, $[[a]] = a^{i} + a^{j}$ ⁠, where $a^{k}$ is the value of $a$ from $\partial Ω_{k}$ [see Fig. 1(a)].

1. Direct flexoelectricity

From Eqs. and , it follows that

\begin{aligned} \bar{Π} [u, ϕ] & = \int_{Ω} (\frac{1}{2} c_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} {\bar{h}}_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ - \frac{1}{2} κ_{k l} E_{k} E_{l} - μ_{l i j k} ε_{i j, k} E_{l}) d Ω \\ - \int_{Ω} W^{Ω} d Ω - \int_{\partial Ω} W^{\partial Ω} d Γ - \int_{C} W^{C} d s, \end{aligned}

(40)

where

W^{Ω}

⁠,

W^{\partial Ω}

⁠, and

W^{C}

represent the work densities exerted by external sources in the domain

Ω

⁠, on its boundary

\partial Ω

and at the non-smooth region

C

of its boundary, respectively. In this case, the admissible external sources of work are⁶³

W^{Ω} (u, ϕ) = {f^{ext}}_{i} u_{i} - q ϕ,

(41a)

W^{\partial Ω} (u, ϕ) = t_{i} u_{i} + r_{i} \partial^{n} u_{i} - w ϕ,

(41b)

W^{C} (u, ϕ) = j_{i} u_{i},

(41c)

where

f^{ext}

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

\partial^{n} \equiv \partial / \partial n

denotes the normal derivative.

Both the boundary of the domain

\partial Ω

and its non-smooth region

C

are split into several disjoint Dirichlet and Neumann sets as

\partial Ω = \partial Ω_{u} \cup \partial Ω_{t} = \partial Ω_{v} \cup \partial Ω_{r} = \partial Ω_{ϕ} \cup \partial Ω_{w},

(42a)

C = C_{u} \cup C_{j},

(42b)

where

\partial Ω_{u}

⁠,

C_{u}

⁠,

\partial Ω_{v}

⁠, and

\partial Ω_{ϕ}

are the Dirichlet regions where either the displacement field, its normal derivative or the electric potential are prescribed,

u = \bar{u} on \partial Ω_{u} and C_{u},

(43a)

\partial^{n} (u) = \bar{v} on \partial Ω_{v},

(43b)

ϕ = \bar{ϕ} on \partial Ω_{ϕ} .

(43c)

The Neumann regions

\partial Ω_{t}

⁠,

C_{j}

⁠,

\partial Ω_{r}

⁠, and

\partial Ω_{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

\partial Ω_{f}

belongs to

\partial Ω_{u}

⁠, its boundaries

\partial \partial Ω_{f}

must belong to

C_{u}

as well, namely,

C_{u} = C \cap \bar{\partial Ω_{u}}

⁠, and hence

C_{j} = C ∖ C_{u}

⁠.

Introducing the definitions and into Eq. yields

\begin{aligned} \bar{Π} [u, ϕ] & = \int_{Ω} (\frac{1}{2} c_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} {\bar{h}}_{i j k l m n} ε_{i j, k} ε_{l m, n} - \frac{1}{2} κ_{k l} E_{k} E_{l} \\ - μ_{l i j k} ε_{i j, k} E_{l} - f_{i}^{ext} u_{i} + q ϕ) d Ω - \int_{\partial Ω_{t}} t_{i} u_{i} d Γ \\ - \int_{\partial Ω_{r}} r_{i} \partial^{n} u_{i} d Γ + \int_{\partial Ω_{w}} w ϕ d Γ - \int_{C_{j}} j_{i} u_{i} d s . \end{aligned}

(44)

According to Eq. , the equilibrium states

(u^{*}, ϕ^{*})

correspond to the saddle points in the enthalpy potential meeting the corresponding Dirichlet boundary conditions , implying that all the admissible variations with respect to the displacement and potential fields must vanish at

(u^{*}, ϕ^{*})

⁠,

\begin{aligned} δ \bar{Π} [u^{*}, ϕ^{*}] & = δ_{u} \bar{Π} [u^{*}, ϕ^{*}] + δ_{ϕ} \bar{Π} [u^{*}, ϕ^{*}] = 0, \\ \forall (δ u, δ ϕ) \in U_{0} \otimes P_{0} . \end{aligned}

(45)

The functional spaces of the admissible variations

U_{0}

⁠,

P_{0}

must meet homogeneous Dirichlet conditions, i.e.,

\begin{aligned} U_{0} & = {δ u \in {[H^{2} (Ω)]}^{D} | δ u = 0 on \partial Ω_{u} and C_{u}, \\ \partial^{n} δ u = 0 on \partial Ω_{v}}, \end{aligned}

(46a)

P_{0} = {δ ϕ \in H^{1} (Ω) | δ ϕ = 0 on \partial Ω_{ϕ}} .

(46b)

Introducing Eq. into yields the weak form of the problem.

Find $(u, ϕ) \in U_{D} \otimes P_{D}$ such that $\forall (δ u, δ ϕ) \in U_{0} \otimes P_{0}$ ⁠,

\begin{aligned} \int_{Ω} ({\hat{σ}}_{i j} δ ε_{i j} + {\tilde{σ}}_{i j k} δ ε_{i j, k} - {\hat{D}}_{l} δ E_{l}) d Ω \\ = \int_{Ω} ({f^{ext}}_{i} δ u_{i} - q δ ϕ) d Ω + \int_{\partial Ω_{t}} t_{i} δ u_{i} d Γ + \int_{\partial Ω_{r}} r_{i} \partial^{n} δ u_{i} d Γ \\ - \int_{\partial Ω_{w}} w δ ϕ d Γ + \int_{C_{j}} j_{i} δ u_{i} d s, \end{aligned}

(47)

with functional spaces

U_{D}

and

P_{D}

fulfilling the Dirichlet boundary conditions in Eq. , that is,

\begin{aligned} U_{D} & = {u \in {[H^{2} (Ω)]}^{D} | u = \bar{u} on \partial Ω_{u} and C_{u}, \\ \partial^{n} u = \bar{v} on \partial Ω_{v}}, \end{aligned}

(48a)

P_{D} = {ϕ \in H^{1} (Ω) | ϕ = \bar{ϕ} on \partial Ω_{ϕ}} .

(48b)

The Cauchy stress

\hat{σ}

⁠, the double stress

\tilde{σ}

⁠, and the electric displacement

\hat{D}

fields in Eq. are the conjugate quantities to the strain

ε

⁠, the strain gradient

\nabla ε

⁠, and the electric field

E

⁠, respectively, as follows:

{\hat{σ}}_{i j} = \frac{\partial {\bar{ψ}}^{(Dir)} (ε, \nabla ε, E)}{\partial ε_{i j}} = c_{i j k l} ε_{k l},

(49a)

{\tilde{σ}}_{i j k} = \frac{\partial {\bar{ψ}}^{(Dir)} (ε, \nabla ε, E)}{\partial ε_{i j, k}} = {\bar{h}}_{i j k l m n} ε_{l m, n} - μ_{l i j k} E_{l},

(49b)

{\hat{D}}_{l} = - \frac{\partial {\bar{ψ}}^{(Dir)} (ε, \nabla ε, E)}{\partial E_{l}} = κ_{l m} E_{m} + μ_{l i j k} ε_{i j, 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, $U_{D}$ and $U_{0}$ 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 $\partial Ω_{v} = \emptyset$ ⁠, which implies that only second order Neumann boundary conditions are allowed, i.e., $\partial Ω = \partial Ω_{r}$ ⁠.^58,64–67 This choice is further justified by the unclear physical interpretation of second order Dirichlet boundary conditions.⁴⁰

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

Ω

⁠, as follows:

\begin{aligned} \int_{Ω} ({\hat{σ}}_{i j} δ ε_{i j} + {\tilde{σ}}_{i j k} δ ε_{i j, k} - {\hat{D}}_{l} δ E_{l}) d Ω \\ = \int_{Ω} (- {\hat{σ}}_{i j, j} δ u_{i} - {\tilde{σ}}_{i j k, k} δ ε_{i j} - {\hat{D}}_{l, l} δ ϕ) d Ω \\ + \int_{\partial Ω} ({\hat{σ}}_{i j} n_{j} δ u_{i} + {\hat{D}}_{l} n_{l} δ ϕ) d Γ + \int_{\partial Ω} {\tilde{σ}}_{i j k} n_{k} δ u_{i, j} d Γ \\ = - \int_{Ω} [{({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k})}_{, j} δ u_{i} + {\hat{D}}_{l, l} δ ϕ] d Ω \\ + \int_{\partial Ω} [({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k}) n_{j} δ u_{i} + {\hat{D}}_{l} n_{l} δ ϕ] d Γ \\ + \int_{\partial Ω} {\tilde{σ}}_{i j k} n_{k} (\nabla_{j}^{S} + n_{j} \partial^{n}) δ u_{i} d Γ \\ = - \int_{Ω} [{({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k})}_{, j} δ u_{i} + {\hat{D}}_{l, l} δ ϕ] d Ω \\ + \int_{\partial Ω} [({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k}) n_{j} δ u_{i} + {\tilde{σ}}_{i j k} n_{j} n_{k} \partial^{n} δ u_{i} + {\hat{D}}_{l} n_{l} δ ϕ] d Γ \\ + \int_{\partial Ω} {\tilde{σ}}_{i j k} n_{k} \nabla_{j}^{S} (δ u_{i}) d Γ, \end{aligned}

(50)

where the gradient of the variations of the displacement field

δ u_{i, j}

in the last integral on

\partial Ω

is expressed in terms of its tangential

\nabla_{j}^{S} δ u_{i}

and normal

n_{j} \partial^{n} δ u_{i}

components, with

\nabla^{S}

being the surface gradient operator. This separation is required since the tangential components of

\nabla δ u

on

\partial Ω

are dependent on

δ u

⁠, whereas the normal component is not. The last term in Eq. is expressed in terms of

δ u

via integration by parts and the surface divergence theorem on

\partial Ω

⁠,

\begin{aligned} \int_{\partial Ω} {\tilde{σ}}_{i j k} n_{k} \nabla_{j}^{S} (δ u_{i}) d Γ & = - \int_{\partial Ω} \nabla_{j}^{S} ({\tilde{σ}}_{i j k} n_{k}) δ u_{i} d Γ \\ + \int_{\partial Ω} \nabla_{j}^{S} ({\tilde{σ}}_{i j k} n_{k} δ u_{i}) d Γ \\ = - \int_{\partial Ω} \nabla_{j}^{S} ({\tilde{σ}}_{i j k} n_{k}) δ u_{i} d Γ \\ + \int_{\partial Ω} \nabla_{l}^{S} (n_{l}) {\tilde{σ}}_{i j k} n_{j} n_{k} δ u_{i} d Γ \\ + \int_{C} [[{\tilde{σ}}_{i j k} n_{k} m_{k}]] δ u_{i} d s . \end{aligned}

(51)

Combining Eqs. , , , and , we obtain

\begin{aligned} \int_{Ω} [[{({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k})}_{, j} + {f^{ext}}_{i}] δ u_{i} + [{\hat{D}}_{l, l} - q] δ ϕ] d Ω \\ = \int_{\partial Ω_{t}} [({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k} + \nabla_{l}^{S} (n_{l}) {\tilde{σ}}_{i j k} n_{k}) n_{j} - \nabla_{j}^{S} ({\tilde{σ}}_{i j k} n_{k}) - t_{i}] δ u_{i} d Γ \\ + \int_{\partial Ω_{r}} [{\tilde{σ}}_{i j k} n_{j} n_{k} - r_{i}] \partial^{n} δ u_{i} d Γ + \int_{\partial Ω_{w}} [{\hat{D}}_{l} n_{l} + w] δ ϕ d Γ \\ + \int_{C_{j}} [[[{\tilde{σ}}_{i j k} n_{k} m_{k}]] - j_{i}] δ u_{i} d s, \end{aligned}

(52)

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

{\begin{matrix} {({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k})}_{, j} + {f^{ext}}_{i} = 0 & in Ω, \\ {\hat{D}}_{l, l} - q = 0 & in Ω, \end{matrix}

(53)

Dirichlet boundary conditions and Neumann boundary conditions

\begin{aligned} ({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k} + \nabla_{l}^{S} (n_{l}) {\tilde{σ}}_{i j k} n_{k}) n_{j} - \nabla_{j}^{S} ({\tilde{σ}}_{i j k} n_{k}) & = t_{i} on \partial Ω_{t}, \end{aligned}

(54a)

{\tilde{σ}}_{i j k} n_{j} n_{k} = r_{i} on \partial Ω_{r},

(54b)

- {\hat{D}}_{l} n_{l} = w on \partial Ω_{w},

(54c)

[[{\tilde{σ}}_{i j k} m_{j} n_{k}]] = j_{i} on C_{j} .

(54d)

In view of Eq. , the definition of the physical stress

σ

arises naturally as

σ_{i j} = {\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k} = c_{i j k l} ε_{k l} - {\bar{h}}_{i j k l m n} ε_{l m, n k} + μ_{l i j k} E_{l, k},

(55)

whereas the physical electric displacement

D

is simply

\hat{D}

⁠.

Remark 2:

The traction

t

can alternatively be written as⁶³

t_{i} = ({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k} - \nabla_{k}^{S} {\tilde{σ}}_{i k j}) n_{j} + {\tilde{σ}}_{i j k} {\tilde{N}}_{j k},

(56)

where

\tilde{N} : = S - Tr (S) n \otimes n

is the second-order geometry tensor written in terms of the shape operator

S

of the surface

\partial Ω

⁠. The first term involves stress measures dotted with the normal vector

n

(a first-order measure of

\partial Ω

⁠, i.e., the orientation), whereas the second term involves the double stress dotted with the second-order geometry tensor (a second-order measure of

\partial Ω

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. , plus the high order dielectric quantities

℘

(electric charge density per unit length) and

r

(double charge density, i.e., charge moment per unit area),⁶² analogous to

j

and

r

from mechanics

W^{Ω} (u, ϕ) = {f^{ext}}_{i} u_{i} - q ϕ,

(57a)

W^{\partial Ω} (u, ϕ) = t_{i} u_{i} + r_{i} \partial^{n} u_{i} - w ϕ - r \partial^{n} ϕ,

(57b)

W^{C} (u, ϕ) = j_{i} u_{i} - ℘ ϕ .

(57c)

Accordingly, the boundary of the domain

\partial Ω

is split into several disjoint Dirichlet and Neumann boundaries as

\begin{aligned} \partial Ω & = \partial Ω_{u} \cup \partial Ω_{t} = \partial Ω_{v} \cup \partial Ω_{r} \\ = \partial Ω_{ϕ} \cup \partial Ω_{w} = \partial Ω_{φ} \cup \partial Ω_{r}, \end{aligned}

(58a)

C = C_{u} \cup C_{j} = C_{ϕ} \cup C_{℘},

(58b)

where the high-order nature of the dielectrics leads to the definition of the Dirichlet

\partial Ω_{φ}

and Neumann

\partial Ω_{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 \cap \bar{\partial Ω_{ϕ}}

and

C_{℘} = C ∖ C_{ϕ}

⁠. The corresponding Dirichlet conditions of the problem are mathematically written as

u = \bar{u} on \partial Ω_{u} and C_{u},

(59a)

\partial^{n} (u) = \bar{v} on \partial Ω_{v},

(59b)

ϕ = \bar{ϕ} on \partial Ω_{ϕ} and C_{ϕ},

(59c)

\partial^{n} (ϕ) = \bar{φ} on \partial Ω_{φ} .

(59d)

With the definitions in and and the enthalpy density in , the enthalpy potential in reads

\begin{aligned} \bar{Π} [u, ϕ] & = \int_{Ω} (\frac{1}{2} c_{i j k l} ε_{i j} ε_{k l} + \frac{1}{2} {\bar{h}}_{i j k l m n} ε_{i j, k} ε_{l m, n} \\ - \frac{1}{2} κ_{k l} E_{k} E_{l} - \frac{1}{2} μ_{l i j k} (ε_{i j, k} E_{l} - ε_{i j} E_{l, k}) \\ - \frac{1}{2} M_{i j k l} E_{i, j} E_{k, l} - f_{i}^{ext} u_{i} + q ϕ) d Ω \\ - \int_{\partial Ω_{t}} t_{i} u_{i} d Γ - \int_{\partial Ω_{r}} r_{i} \partial^{n} u_{i} d Γ - \int_{C_{j}} j_{i} u_{i} d s \\ + \int_{\partial Ω_{w}} w ϕ d Γ + \int_{\partial Ω_{r}} r \partial^{n} ϕ d Γ + \int_{C_{℘}} ℘ ϕ d s . \end{aligned}

(60)

The variational principle in Eq. implies that feasible equilibrium states

u^{*} \in U_{D}

⁠,

ϕ^{*} \in P_{D}

must meet Dirichlet boundary conditions , and their corresponding admissible variations

δ u \in U_{0}

⁠,

δ ϕ \in P_{0}

must vanish as stated in Eq. , where

\begin{aligned} U_{D} & = {u \in {[H^{2} (Ω)]}^{D} | u = \bar{u} on \partial Ω_{u} and C_{u}, \\ \partial^{n} u = \bar{v} on \partial Ω_{v}}, \\ P_{D} & = {ϕ \in H^{2} (Ω) | ϕ = 0 on \partial Ω_{ϕ} and C_{ϕ}, \end{aligned}

(61a)

\partial^{n} ϕ = \bar{φ} on \partial Ω_{φ}},

(61b)

and

\begin{aligned} U_{0} & = {δ u \in {[H^{2} (Ω)]}^{D} | δ u = 0 on \partial Ω_{u} and C_{u}, \\ \partial^{n} δ u = 0 on \partial Ω_{v}}, \\ P_{0} & = {δ ϕ \in H^{2} (Ω) | δ ϕ = 0 on \partial Ω_{ϕ} and C_{ϕ}, \end{aligned}

(62a)

\partial^{n} δ ϕ = 0 on \partial Ω_{φ}} .

(62b)

Introducing Eq. into yields the weak form of the Lifshitz-invariant flexoelectricity problem

Find $(u, ϕ) \in U_{D} \otimes P_{D}$ such that, $\forall (δ u, δ ϕ) \in U_{0} \otimes P_{0}$ ⁠,

\begin{aligned} \int_{Ω} ({\hat{σ}}_{i j} δ ε_{i j} + {\tilde{σ}}_{i j k} δ ε_{i j, k} - {\hat{D}}_{l} δ E_{l} - {\tilde{D}}_{l m} δ E_{l, m}) d Ω \\ = \int_{Ω} ({f^{ext}}_{i} δ u_{i} - q δ ϕ) d Ω \\ + \int_{\partial Ω_{t}} t_{i} δ u_{i} d Γ + \int_{\partial Ω_{r}} r_{i} \partial^{n} δ u_{i} d Γ + \int_{C_{j}} j_{i} δ u_{i} d s \\ - \int_{\partial Ω_{w}} w δ ϕ d Γ - \int_{\partial Ω_{r}} r \partial^{n} δ ϕ d Γ - \int_{C_{℘}} ℘ δ ϕ d s . \end{aligned}

(63)

The Cauchy stress

\hat{σ}

⁠, the double stress

\tilde{σ}

⁠, the electric displacement

\hat{D}

⁠, and the double electric displacement

\tilde{D}

in Eq. are the conjugate quantities to the strain

ε

⁠, the strain gradient

\nabla ε

⁠, the electric field

E

⁠, and the electric field gradient

\nabla E

⁠, respectively, as follows:

{\hat{σ}}_{i j} = \frac{\partial {\bar{ψ}}^{(Lif)}}{\partial ε_{i j}} = c_{i j k l} ε_{k l} + \frac{1}{2} μ_{l i j k} E_{l, k},

(64a)

{\tilde{σ}}_{i j k} = \frac{\partial {\bar{ψ}}^{(Lif)}}{\partial ε_{i j, k}} = {\bar{h}}_{i j k l m n} ε_{l m, n} - \frac{1}{2} μ_{l i j k} E_{l},

(64b)

{\hat{D}}_{l} = - \frac{\partial {\bar{ψ}}^{(Lif)}}{\partial E_{l}} = κ_{l m} E_{m} + \frac{1}{2} μ_{l i j k} ε_{i j, k},

(64c)

{\tilde{D}}_{l k} = - \frac{\partial {\bar{ψ}}^{(Lif)}}{\partial E_{l, k}} = M_{i j l k} E_{i, j} - \frac{1}{2} μ_{l i j k} ε_{i j} .

(64d)

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

{\begin{matrix} {({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k})}_{, j} + {f^{ext}}_{i} = 0 in Ω, \\ {({\hat{D}}_{l} - {\tilde{D}}_{l k, k})}_{, l} - q = 0 in Ω, \end{matrix}

(65)

complemented with Dirichlet conditions and Neumann conditions

\begin{aligned} ({\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k} + \nabla_{l}^{S} (n_{l}) {\tilde{σ}}_{i j k} n_{k}) n_{j} \\ - \nabla_{j}^{S} ({\tilde{σ}}_{i j k} n_{k}) & = t_{i} on \partial Ω_{t}, \end{aligned}

(66a)

[[{\tilde{σ}}_{i j k} m_{j} n_{k}]] = j_{i} on C_{j},

(66b)

\begin{aligned} {\tilde{σ}}_{i j k} n_{j} n_{k} & = r_{i} on \partial Ω_{r}, \\ - ({\hat{D}}_{l} - {\tilde{D}}_{l k, k} + \nabla_{i}^{S} (n_{i}) {\tilde{D}}_{l k} n_{k}) n_{l} \\ + \nabla_{l}^{S} ({\tilde{D}}_{l k} n_{k}) = w on \partial Ω_{w}, \end{aligned}

(66c)

- [[{\tilde{D}}_{j k} m_{j} n_{k}]] = ℘ on C_{℘},

(66d)

- {\tilde{D}}_{j k} n_{j} n_{k} = r on \partial Ω_{r} .

(66e)

The physical stress

σ

and physical electric displacement

D

arise from Eq. as

σ_{i j} = {\hat{σ}}_{i j} - {\tilde{σ}}_{i j k, k} = c_{i j k l} ε_{k l} - {\bar{h}}_{i j k l m n} ε_{l m, n k} + μ_{l i j k} E_{l, k},

(67a)

D_{l} = {\hat{D}}_{l} - {\tilde{D}}_{l k, k} = κ_{l m} E_{m} - M_{i j l k} E_{i, j k} + μ_{l i j k} ε_{i j, 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 $\partial Ω$ 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 $C^{1}$ continuity requirement.⁵⁹

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 method⁶⁸ 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⁶⁹) as compared to other alternatives such as the Lagrange multipliers or penalty methods.

Following Nitsche’s approach, an additional term

{\bar{Π}}^{Nitsche} [u, ϕ]

acting on the Dirichlet boundaries is introduced in the enthalpy functional in Eq. or Eq. to weakly enforce Dirichlet boundary conditions Eq. or Eq. , respectively. In the case of Lifshitz-invariant formulation, it has the following form:

\begin{aligned} {\bar{Π}}^{Nitsche} [u, ϕ] \\ = \int_{\partial Ω_{u}} (\frac{1}{2} β_{u} {(u_{i} - {\bar{u}}_{i})}^{2} - (u_{i} - {\bar{u}}_{i}) t_{i} (u, ϕ)) d Γ \\ + \int_{\partial Ω_{v}} (\frac{1}{2} β_{v} {(\partial^{n} u_{i} - {\bar{v}}_{i})}^{2} - (\partial^{n} u_{i} - {\bar{v}}_{i}) r_{i} (u, ϕ)) d Γ \\ + \int_{C_{u}} (\frac{1}{2} β_{C_{u}} {(u_{i} - {\bar{u}}_{i})}^{2} - (u_{i} - {\bar{u}}_{i}) j_{i} (u, ϕ)) d s \\ - \int_{\partial Ω_{ϕ}} (\frac{1}{2} β_{ϕ} {(ϕ - \bar{ϕ})}^{2} - (ϕ - \bar{ϕ}) w (u, ϕ)) d Γ \\ - \int_{\partial Ω_{φ}} (\frac{1}{2} β_{φ} {(\partial^{n} ϕ - \bar{φ})}^{2} - (\partial^{n} ϕ - \bar{φ}) r (u, ϕ)) d Γ \\ - \int_{C_{u}} (\frac{1}{2} β_{C_{ϕ}} {(ϕ - \bar{ϕ})}^{2} - (ϕ - \bar{ϕ}) ℘ (u, ϕ)) d s, \end{aligned}

(68)

with the numerical parameters

β_{u}, β_{v}, β_{C_{u}}, β_{ϕ}, β_{φ}, β_{C_{ϕ}} \in R^{+}

and Neumann terms explicitly depending on the state variables as defined in Eq. . For the direct formulation, the last two lines in Eq. must be omitted according to the Dirichlet conditions . 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

ϕ

⁠.⁶³

Remark 3:

Nitsche’s formulations are self-consistent for any value of the penalty parameters

β_{u}, β_{v}, β_{C_{u}}, β_{ϕ}, β_{φ}, β_{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

min

–

max

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.⁶³ 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

ζ \in R^{+}

as follows:

\begin{matrix} β_{u} = \frac{Y}{h} ζ, β_{v} = (ℓ_{mech}^{2} + ℓ_{μ}^{2}) β_{u}, β_{C_{u}} = \frac{β_{v}}{h}, \\ β_{ϕ} = \frac{ε}{h} ζ, β_{φ} = (ℓ_{elec}^{2} + ℓ_{μ}^{2}) β_{ϕ}, β_{C_{ϕ}} = \frac{β_{φ}}{h}, \end{matrix}

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

ℓ_{μ} \sim μ / \sqrt{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

ζ = 10

–

100

typically provides stable results.⁶³

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

\begin{aligned} U_{D} = U_{0} & = {[H^{2} (Ω)]}^{D} with L^{2} - integrable third \\ derivatives on the boundary \partial Ω_{u}, \end{aligned}

(70a)

P_{D} = P_{0} = H^{1} (Ω),

(70b)

where

U_{D}

and

U_{0}

account for the integrals involving

t

in Eq. , requiring the computation of third order derivatives. In practice,

U_{D} = U_{0} = {[H^{3} (Ω)]}^{D}

are typically considered.⁶³ Analogously, the Lifshitz-invariant form requires

\begin{aligned} U_{D} = U_{0} & = {[H^{2} (Ω)]}^{D} with L^{2} - integrable third \\ derivatives on the boundary \partial Ω_{u}, \end{aligned}

(71a)

\begin{aligned} P_{D} = P_{0} & = H^{2} (Ω) with L^{2} - integrable third \\ derivatives on the boundary \partial Ω_{ϕ} . \end{aligned}

(71b)

Remark 4:

Many authors in computational flexoelectricity neglect the edge Dirichlet conditions (43a) on $C_{u}$ ^9,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 $\partial Ω_{u}$ (⁠ $\partial Ω_{ϕ}$ ⁠) automatically implies the strong imposition on the adjacent edges in $C_{u}$ (⁠ $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.⁶³

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

\partial Ω_{Φ}^{i}

and its edge

\partial \partial Ω_{Φ}^{i}

is mathematically written as

ϕ = Φ^{i} on \partial Ω_{Φ}^{i} and \partial \partial Ω_{Φ}^{i}, \forall i = 1, \dots, N_{Φ},

(72)

where

Φ^{1}, \dots, Φ^{N_{Φ}} \in R

are additional state variables and

N_{Φ}

denotes the number of sensing electrodes in the system. Since each

Φ^{i}

is not prescribed on

\partial Ω_{Φ}^{i}

nor

\partial \partial Ω_{Φ}^{i}

⁠, their energy conjugates—i.e., the total surface charge density on the electrode⁶³ and on their edges—must vanish, that is,

\int_{\partial Ω_{Φ}^{i}} w d Ω = 0 on \partial Ω_{Φ}^{i}, \forall i = 1, \dots, N_{Φ},

(73a)

\int_{\partial \partial Ω_{Φ}^{i}} ℘ d s = 0 on \partial \partial Ω_{Φ}^{i}, \forall i = 1, \dots, N_{Φ} .

(73b)

In the direct flexoelectricity model, high-order electrostatics are not present, hence Eqs. and must hold on

\partial Ω_{Φ}^{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 multipliers⁷² or a modified version of Nitsche’s method.⁶³

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.⁸¹ 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.

View large Download slide

(a) Composite domain $Ω = Ω_{1} \cup Ω_{2} \cup Ω_{3} .$ (b) Detail of the interface $I$ between $Ω_{1}$ and $Ω_{2}$ ⁠. (c) Detail of the point $S$ shared by all three subdomains.

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

u^{L} = u^{R} on I,

(74a)

u^{i} = u^{j} on S,

(74b)

\partial u^{L} / \partial n^{L} = \partial u^{R} / \partial n^{L} on I,

(74c)

ϕ^{L} = ϕ^{R} on I,

(74d)

ϕ^{i} = ϕ^{j} on S,

(74e)

\partial ϕ^{L} / \partial n^{L} = \partial ϕ^{R} / \partial n^{L} on I,

(74f)

and high-order equilibrium conditions read

t^{L} + t^{R} = 0 on I,

(75a)

\sum_{k} j^{k} = 0 on S,

(75b)

r^{L} - r^{R} = 0 on I,

(75c)

w^{L} + w^{R} = 0 on I,

(75d)

\sum_{k} ℘^{k} = 0 on S,

(75e)

r^{L} - r^{R} = 0 on I,

(75f)

where

L

and

R

represent either sides across the interface

I

⁠,

i \neq j

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. , , , and must not be considered.

The conditions in Eqs. (74) and (75) are implicitly fulfilled in $C^{1}$ -conforming body-fitted discretizations. Otherwise, they must be explicitly enforced, e.g., weakly via Nitsche’s method.⁸¹

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 $(L_{1} e_{1}, L_{2} e_{2}, L_{3} e_{3})$ in $R^{3}$ ⁠. Without loss of generality, the flexoelectricity boundary value problem can be reduced to a single RVE defined as $Ω \in (0, L_{1}) \times (0, L_{2}) \times (0, L_{3})$ with the proper generalized periodic conditions (i) between a given periodic boundary $A$ of the RVE (i.e., $x_{1} = 0$ ⁠, $x_{2} = 0$ ⁠, or $x_{3} = 0$ ⁠) and its corresponding translational image $A^{'}$ (i.e., $x_{1} = L_{1}$ ⁠, $x_{2} = L_{2}$ ⁠, or $x_{3} = L_{3}$ ⁠), 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.

View large Download slide

(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.

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

u^{A} = u^{A^{'}} - {[[u]]}^{A} on A and S,

(76a)

\partial u^{A} / \partial n^{A} = \partial u^{A^{'}} / \partial n^{A} on A,

(76b)

ϕ^{A} = ϕ^{A^{'}} - {[[ϕ]]}^{A} on A and S,

(76c)

\partial ϕ^{A} / \partial n^{A} = \partial ϕ^{A^{'}} / \partial n^{A} on A,

(76d)

t^{A} + t^{A^{'}} = 0 on A,

(76e)

j^{S} + j^{S^{'}} = 0 on S,

(76f)

r^{A} - r^{A^{'}} = 0 on A,

(76g)

w^{A} + w^{A^{'}} = 0 on A,

(76h)

℘^{S} + ℘^{S^{'}} = 0 on S,

(76i)

r^{A} - r^{A^{'}} = 0 on A,

(76j)

where

{[[\cdot]]}^{A}

are either known (enforced) or unknown (resolved) constant quantities denoting the difference of the state variable within the brackets

(\cdot)

across either sides of the RVE along the direction orthogonal to the plane

A

⁠. For the direct flexoelectricity model, Eqs. , and are disregarded and Eq. 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.⁸³ 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.⁸² Alternatively, one can resort to enforcing generalized periodicity conditions (76) weakly.⁸¹

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 $Ω_{0}$ in the reference (or undeformed) configuration and by $Ω$ in the current (or deformed) configuration. The deformation map $χ : Ω_{0} \to Ω$ maps every material point $X \in Ω_{0}$ to the spatial point $x = χ (X) \in Ω$ ⁠. Relevant deformation measures in this context are the deformation gradient $F = \nabla_{0} χ$ ⁠, the Jacobian determinant $J = det (F)$ ⁠, and the right Cauchy–Green deformation tensor $C = F^{T} \cdot F$ ⁠. In the reference configuration, the strain is represented by the Green–Lagrangian strain tensor $E = (C - I) / 2$ ⁠, and the electric field in by $E = - \nabla_{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.,⁴¹ resulting in

\begin{aligned} {\bar{Ψ}}^{Enth} (E, \nabla_{0} E, E) \\ = Ψ^{Mech} (E, \nabla_{0} E) - \frac{1}{2} E \cdot \bar{κ} \cdot E - E \cdot \bar{μ} ⋮ \nabla_{0} E, \end{aligned}

(77)

where

Ψ^{Mech}

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

\bar{κ} (E) = J C^{- 1} ε,

(78a)

\bar{μ} (E) = J C^{- 1} \cdot μ

(78b)

are the effective dielectricity and flexoelectricity tensors,⁴¹ 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. , as⁴¹

\begin{aligned} \bar{Π} [χ, Φ] & = \int_{Ω_{0}} ({\bar{Ψ}}^{Enth} (E, \nabla_{0} E, E) - {F^{ext}}_{i} χ_{i} + Q Φ) d Ω_{0} \\ - \int_{\partial Ω_{0}^{T}} T_{i} χ_{i} d Γ_{0} - \int_{\partial Ω_{0}^{R}} R_{i} \partial_{0}^{N} χ_{i} d Γ_{0} \\ - \int_{C_{0}^{J}} J_{i} χ_{i} d s_{0} + \int_{\partial Ω_{0}^{W}} W Φ d Γ_{0}, \end{aligned}

(79)

where

F^{ext}

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

\bar{Π}

satisfying

{χ^{*}, Φ^{*}} = \arg min_{χ \in X} max_{Φ \in P} \bar{Π} [χ, Φ],

(80)

where

X

and

P

are the functional spaces for

χ

and

Φ

with sufficient regularity fulfilling Dirichlet boundary conditions. Vanishing the first variations of

\bar{Π} [χ, Φ]

leads to the weak form of the boundary value problem, whose corresponding Euler–Lagrange equations read as follows:⁴¹

{(F_{i I} S_{I J})}_{, J} + {F^{ext}}_{i} = 0,

(81a)

D_{K, K} - Q = 0.

(81b)

The physical second Piola–Kirchhoff stress

S

and the electric displacement

D

are defined in terms of

χ

and

Φ

as

S_{I J} (χ, Φ) = {S^{Mech}}_{I J} (χ, Φ) + {S^{Maxwell}}_{I J} (χ, Φ),

(82a)

D_{L} (χ, Φ) = J C_{K L}^{- 1} (ε E_{K} + μ_{K I J M} E_{I J, M}),

(82b)

where

\begin{aligned} {S^{Mech}}_{I J} & = \frac{\partial Ψ^{Mech}}{\partial E_{I J}} - \frac{\partial}{\partial X_{K}} (\frac{\partial Ψ^{Mech}}{\partial E_{I J, K}}) \\ + J C_{L M}^{- 1} E_{M, K} μ_{L I J K} \end{aligned}

(83)

and

D

correspond to the extensions of Eqs. and to finite deformations, respectively, and the additional term

{S^{Maxwell}}_{I J} = \frac{J}{2} (C_{M I}^{- 1} C_{L J}^{- 1} + C_{L I}^{- 1} C_{M J}^{- 1} - C_{M L}^{- 1} C_{I J}^{- 1}) E_{M} ε E_{L}

(84)

is the total second Piola–Maxwell stress tensor,⁴¹ arising as a direct consequence of the strain-dependence of the effective dielectricity tensor in Eq. . 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.⁴¹

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,⁵⁵ or Cosserat rod⁴¹ models. Otherwise, it is necessary to resort to computational flexoelectricity.

The major challenge is to handle the $C^{1}$ 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 $C^{1}$ Argyris triangular element approximation,³⁸ the immersed boundary B-spline method,⁶³ or the C0-interior penalty method,⁷² among others. Another family of numerical methods are those circumventing the $C^{1}$ 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.⁵⁹ 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 $C^{p - 1}$ -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

ξ \in Ξ = [k_{0}, k_{m}]

in terms of the knot vector

k = [k_{0}, k_{1}, k_{2}, \dots, k_{m}]

⁠, where

k_{i} < k_{i + 1}

⁠. The

i

th function of this basis is defined recursively as⁹⁷

\begin{aligned} B_{i}^{0} (ξ) & = {\begin{matrix} 1 & k_{i - 1} \leq ξ < k_{i}, \\ 0 & otherwise, \end{matrix} \\ B_{i}^{p} (ξ) & = \frac{ξ - k_{i - 1}}{k_{i + p - 1} - k_{i - 1}} B_{i}^{p - 1} (ξ) + \frac{k_{i + p} - ξ}{k_{i + p} - k_{i}} B_{i + 1}^{p - 1} (ξ), \\ i = 1, \dots, n . \end{aligned}

(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 $ξ = k_{0}$ to $ξ = k_{p + 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 $k_{i} = 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 $B_{i}^{p} (ξ) = B_{1}^{p} (ξ - i + 1)$ (see Fig. 5).

FIG. 5.

View large Download slide

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

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 $k_{0}$ and $k_{m}$ in parametric space. To do so, the basis is modified by knot multiplicity, i.e., $k = [k_{0}, \dots, k_{0}, k_{1}, k_{2}, \dots, k_{m}, \dots, k_{m}]$ ⁠, where $k_{0}$ and $k_{m}$ are repeated $p + 1$ times. In doing so, the continuity of the B-spline basis at these knots is decreased to $C^{- 1}$ (i.e., discontinuous), yielding a boundary-interpolant (or open) basis suitable to enforce essential boundary conditions strongly (see Fig. 6).

FIG. 6.

View large Download slide

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]$ ⁠.

In a

D

-dimensional space, the

i

th B-spline function

B_{i}^{p} (ξ)

of a

D

-variate B-spline basis (where

i

is the

D

-variate index

[i_{1}, \dots, i_{D}]

⁠) is defined as the tensor product of

D

univariate B-spline functions as

\begin{aligned} B_{i}^{p} (ξ) = B_{[i_{1}, \dots, i_{D}]}^{p} ([ξ_{1}, \dots, ξ_{D}]) = \prod_{d = 1}^{D} B_{i_{d}}^{p} (ξ_{d}), \\ i_{d} = 1, \dots, n_{d}, \end{aligned}

(86)

which is defined on the

D

-dimensional parametric space

ξ \in Ξ = [k_{1, 0}, k_{1, m_{1}}] \otimes \dots \otimes [k_{d, 0}, k_{d, m_{d}}]

⁠.

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

\begin{matrix} φ : Ξ \to Ω \\ ξ \mapsto x = φ (ξ), \end{matrix}

(87)

which maps a given point

ξ \in Ξ

in the parametric space to a given point

x \in Ω

in the physical space. The basis functions

N (x) \in Ω

are defined as

N = [B^{p} ° φ^{- 1}]

⁠, in such a way that

{[u (x)]}_{d} ≃ {[u^{h} (x)]}_{d} = N_{i} (x) {a^{u}}_{i d} = B_{i}^{p} (ξ) {a^{u}}_{i d},

(88a)

ϕ (x) ≃ ϕ^{h} (x) = N_{i} (x) {a^{ϕ}}_{i} = B_{i}^{p} (ξ) {a^{ϕ}}_{i},

(88b)

for

d = 1, \dots, D

⁠, where

{a^{u}, a^{ϕ}}

are the degrees of freedom (DOF) of the approximations

u^{h} (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.,

Ω = x_{0} + [0, L_{1}] \otimes \dots \otimes [0, L_{D}]

⁠, the map

φ

simply corresponds to an affine transformation composed by a non-uniform scaling and a translation,

x = φ (ξ) = x_{0} + J \cdot ξ, J = [\begin{matrix} h_{1} & 0 & 0 \\ 0 & ⋱ & 0 \\ 0 & 0 & h_{D} \end{matrix}],

(89)

where

h_{d}

is the size of each cell in physical space along the dimension

d

⁠, determined as

h_{d} = L_{d} / m_{d}

⁠.

This simple expression, with a constant and diagonal Jacobian matrix

J

⁠, facilitates the computation of (high-order) gradients of the basis functions

N (x) \in Ω

⁠, since the tensor product structure of the basis is preserved, e.g.,

{\nabla^{x}}_{j} N_{i} (x) = h_{j}^{- 1} (\frac{\partial}{\partial ξ_{j}} B_{i_{j}}^{p} (ξ_{j})) \prod_{d \neq j}^{D} B_{i_{d}}^{p} (ξ_{d}),

(90a)

\begin{aligned} {\nabla^{x}}_{j} [{\nabla^{x}}_{j} N_{i}] (x) = h_{j}^{- 2} (\frac{\partial^{2}}{\partial^{2} ξ_{j}} B_{i_{j}}^{p} (ξ_{j})) \prod_{d \neq j}^{D} B_{i_{d}}^{p} (ξ_{d}), \\ {\nabla^{x}}_{k} [{\nabla^{x}}_{j} N_{i}] (x) \\ = h_{j}^{- 1} h_{k}^{- 1} (\frac{\partial}{\partial ξ_{j}} B_{i_{j}}^{p} (ξ_{j})) (\frac{\partial}{\partial ξ_{k}} B_{i_{k}}^{p} (ξ_{k})) \prod_{d \neq j, d \neq k}^{D} B_{i_{d}}^{p} (ξ_{d}), \end{aligned}

(90b)

for j \neq k,

(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 $L_{2}$ 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 $a^{u}$ and $a^{ϕ}$ ⁠, and hence, the approximated fields $u^{h} (x)$ and $ϕ^{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

Ω

is embedded into a larger domain

Ω_{□}

with rectangular/cuboidal shape, i.e.,

Ω \subseteq Ω_{□}

⁠. The geometrical map in Eq. is redefined as

\begin{matrix} φ : Ξ \to Ω_{□} \\ ξ \mapsto x = φ (ξ), \end{matrix}

(91)

which is independent on

Ω

⁠, and hence, arbitrary geometries are allowed. The basis functions

N (x) \in Ω_{□}

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. – hold.

FIG. 7.

Sketch of the immersed boundary method. (a) Physical domain Ω ∈ R 2 (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.

View large Download slide

Sketch of the immersed boundary method. (a) Physical domain $Ω \in R^{2}$ (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.

Remark 6:

Since the approximation space is unfitted to the geometry, it is not interpolant on $\partial Ω$ ⁠, 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 $\partial Ω$ 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} \subseteq Ω$ ⁠, 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ć,¹⁰⁰ Fries,¹⁰¹ Kudela et al.,¹⁰² Legrain et al.¹⁰³). 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 $\partial Ω$ ⁠, 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 $Ω_{□}^{c}$ since only their intersection with $Ω$ ⁠, i.e., $Ω_{□}^{c} \cap Ω$ ⁠, is to be considered. The standard approach consists on dividing $Ω_{□}^{c} \cap Ω$ 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,¹⁰⁶ 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 $2^{D}$ 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.

View large Download slide

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

FIG. 9.

View large Download slide

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

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 $\partial Ω$ is considered. In the case of implicit boundary representations such as level sets, an explicit approximation of $\partial Ω$ 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} \cap Ω | / | Ω_{□}^{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} = min_{c} (η_{c})$ in meshes of fixed size $h$ at a rate of $η_{\min}^{- (2 p + 1 - 2 / D)}$ ⁠,⁶² 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).$

View large Download slide

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$ ⁠).

Several strategies have been proposed to alleviate ill-conditioning of trimmed cells, such as the ghost penalty method,¹¹² the artificial stiffness approach,^110,111 the extended B-spline method,^113–116 or special preconditioning techniques specifically designed for immersed boundary methods,¹¹⁷ 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.¹¹³ 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.

View large Download slide

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 ${{\tilde{B}}_{3}^{2} (ξ)$ ⁠, ${\tilde{B}}_{4}^{2} (ξ)$ ⁠, ${\tilde{B}}_{5}^{2} (ξ)}$ ⁠, as a result of the extrapolation of the basis ${B_{3}^{2} (ξ)$ ⁠, $B_{4}^{2} (ξ)$ ⁠, $B_{5}^{2} (ξ)}$ in $Ω_{□}^{4}$ toward $Ω_{□}^{5}$ ⁠. The critical basis $B_{6}^{2} (ξ)$ is removed.

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,¹¹⁸ 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 B i 3 , 1 ( ξ 1 ), i = 3 , … , 6 in the coarsest level of refinement ( ℓ = 1) have been replaced by their children B j 3 , 2 ( ξ 2 ), j = 5 , … , 15 in the next ( ℓ = 2) refinement level.

View large Download slide

Univariate uniform cubic (⁠ $p = 3$ ⁠) hierarchical B-spline basis. The basis functions $B_{i}^{3, 1} (ξ^{1})$ ⁠, $i = 3, \dots, 6$ in the coarsest level of refinement (⁠ $ℓ = 1$ ⁠) have been replaced by their children $B_{j}^{3, 2} (ξ^{2})$ ⁠, $j = 5, \dots, 15$ in the next (⁠ $ℓ = 2$ ⁠) refinement level.

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¹¹⁹

\begin{aligned} B_{i}^{p, ℓ} (ξ^{ℓ}) & = \sum_{j = 1}^{p + 2} {[s^{p}]}_{j} B_{2 (i - 1) + j}^{p, ℓ} (2 ξ^{ℓ}) \\ = \sum_{j = 2 i - 1}^{2 i + p} {[s^{p}]}_{j - 2 (i - 1)} B_{j}^{p, ℓ + 1} (ξ^{ℓ + 1}), \end{aligned}

(92)

where

ℓ \in N

stands for the refinement level,

ξ^{ℓ + 1} = 2 ξ^{ℓ}

and the scaling coefficients

{[s^{p}]}_{j} = \frac{1}{2^{p}} (\begin{matrix} p + 1 \\ j - 1 \end{matrix}) = \frac{2^{- p} (p + 1)!}{(j - 1)! (p + 2 - j)!},

(93)

for

j = {1, \dots, p + 2}

⁠. This process maintains the properties of the spanned functional space and, in particular, its regularity.

FIG. 13.

View large Download slide

Two-scale relation of univariate B-splines. Top: original (parent) B-spline $B_{1}^{p, 1} (ξ^{1})$ at level $ℓ = 1$ ⁠. Bottom: the children B-spline basis functions $B_{j}^{p, 2} (ξ^{2})$ at level $ℓ = 2$ scaled by the factors ${[s^{p}]}_{j}$ ⁠, for $j = {1, \dots, p + 2}$ ⁠, (a) linear (⁠ $p = 1$ ⁠), $s^{1} = \frac{1}{2} [1, 2, 1]$ ⁠, (b) quadratic (⁠ $p = 2$ ⁠), $s^{2} = \frac{1}{4} [1, 3, 3, 1]$ ⁠, (c) cubic (⁠ $p = 3$ ⁠), $s^{3} = \frac{1}{8} [1, 4, 6, 4, 1]$ ⁠.

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 $K_{u u}$ 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 $10^{9}$ (or higher) and $10^{- 9}$ (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

f_{L}

⁠,

f_{S}

⁠, and

f_{P}

for each of them, the normalized aforementioned magnitudes are defined as

u_{0} = u / f_{L}, σ_{0} = σ / f_{S}, P_{0} = P / f_{P} .

(94)

All the other quantities in the problem are divided by normalization factors that are derived as a combination of

f_{L}

⁠,

f_{S}

⁠, and

f_{P}

⁠, e.g.,

the electric potential $ϕ_{0} = ϕ / (f_{L} f_{S} f_{P}^{- 1})$ ⁠,
Young’s modulus $Y_{0} = Y / f_{S}$ ⁠,
the flexoelectric coefficients $μ_{0} = μ / (f_{P} f_{L})$ ⁠,
the electric permittivity $ε_{0} = ε / (f_{P}^{2} f_{S}^{- 1})$ ⁠,

and so on. Typically, suitable values for the normalizing factors are $f_{L} \sim L_{geom}$ ⁠, $f_{S} \sim Y$ ⁠, and $f_{P} \sim \sqrt{Y ε}$ ⁠, where $L_{geom}$ 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.

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 μ 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 / V m$ ⁠. 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.

View large Download slide

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$ ⁠.

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⁴¹ 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.³⁵ 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.

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³ to experimentally demonstrate for the first time the flexoelectric effect, which had been predicted theoretically by Mashkevich and Tolpygo.¹¹ Computational studies are also present in Abdollahi et al.⁵⁸ and Zhuang et al.⁵⁹

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.

View large Download slide

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$ ⁠.

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

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.⁴¹ This device was considered by Bursian and Zaikovskii³ 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 = 100 nm$ ⁠, 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 $E_{Y} = - 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.

View large Download slide

Actuation of soft cantilever rod with different flexoelectric tensors $μ$ (nC/m). (a) Geometry and boundary conditions. (b) Axial strain $ζ (V)$ ⁠. (c) Curvature $R^{- 1} (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.

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.,⁴¹ 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:

\begin{aligned} Ψ^{Mech} (E, \nabla_{0} E) & = \frac{λ}{2} ({[\log (J)]}^{2} + ℓ^{2} E_{I I, K} E_{J J, K}) \\ + μ (Tr (E) + ℓ^{2} E_{I J, K} E_{I J, K}), \end{aligned}

(95)

with Young’s modulus

Y = 1.0

GPa, Poisson’s ratio

ν = 0.37

(that is, Lamé’s moduli

λ \approx 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.⁴¹ 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)].

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.⁸³ 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 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.

FIG. 17.

View large Download slide

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.

TABLE I.

Material parameters in Sec. VI A.

Material

Y

ν

ℓ_mech

ɛ

μ_L

μ_T

μ_S

(GPa)

—

(nm)

(nC/V m)

(μC/m)

Matrix

152

0.33

1

141

150

110

Insulators

152

0.33

1

0

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.

B. Geometrically-polarized periodic inclusions

According to Sharma et al.⁷⁵ and Mocci et al.,⁸³ 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.⁸³

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.³⁵

FIG. 18.

View large Download slide

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 $\hat{ϕ} = ϕ ε / μ$ ⁠.

TABLE II.

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

Y

ν

ℓ_mech

ɛ

ℓ_elec

μ_L

μ_T

μ_S

(GPa)

—

(nm)

(nC/V m)

(nm)

(μ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]]}^{x_{1}} = {[[ϕ]]}^{x_{1}} = 0$ ⁠), whereas the jump in electric potential along the vertical direction (⁠ ${[[ϕ]]}^{x_{2}}$ ⁠) is resolved as a result of a prescribed vertical jump in the displacement field (⁠ ${[[u]]}^{x_{2}} = \bar{{[[u]]}^{x_{2}}}$ ⁠).

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.⁸³

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.⁸³ 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.⁸³ 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 (⁠ ${[[ϕ]]}^{x_{2}} = \bar{{[[ϕ]]}^{x_{2}}}$ ⁠), 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 $L_{1} = 4.8 μ$ m and $L_{2} = 2.8 μ$ m and the beam thickness is $t = 160$ nm. The resulting vertical displacement (⁠ ${[[u]]}^{x_{2}}$ ⁠) 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.

View large Download slide

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 $\hat{ϕ} = ϕ ε / μ$ ⁠, and the displacement field as $\hat{u} = u / L_{2} .$

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,⁸³ 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.³⁵ 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.³⁶ 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.

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 $κ$ ⁠, 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:

c_{I J K L} = λ δ_{I J} δ_{K L} + 2 μ δ_{I K} δ_{J L} .

(A1)

Remark 9:

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

b. Strain gradient elasticity tensor h

We consider an isotropic simplified strain gradient elasticity tensor,¹³² which depends on

λ

⁠,

μ

and the mechanical length scale

ℓ_{mech}

in the following form:

h_{I J K L M N} = (λ δ_{I J} δ_{L M} + 2 μ δ_{I L} δ_{J M}) ℓ_{mech}^{2} δ_{K N} .

(A2)

c. Dielectricity tensor κ

Isotropic dielectricity is represented by the second-order tensor

κ

⁠, which depends on the electric permittivity

ε

as

κ_{i j} = ε δ_{i j} .

(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

M_{i j k l} = ε ℓ_{elec}^{2} δ_{i k} δ_{j l} .

(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:

{μ^{< x_{1} >}}_{L I J K} = {\begin{matrix} μ_{L}, & for L = I = J = K, \\ μ_{T}, & for I = J \neq K = L, \\ μ_{S}, & for L = I \neq J = K or L = J \neq I = K, \\ 0 & otherwise . \end{matrix}

(A5)

The flexoelectric tensor

μ

oriented in an arbitrarily rotated orthonormal basis is obtained by rotating

μ^{⟨ x_{1} ⟩}

as

μ_{l i j k} = R_{l L} R_{i I} R_{j J} R_{k K} {μ^{⟨ x_{1} ⟩}}_{L I J K},

(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}$ ⁠.

APPENDIX B: RESIDUAL-BASED WEAK FORM STABILIZATION

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 $ε = 11 nJ / V^{2} m$ ⁠, 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.

View large Download slide

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$ ⁠.

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

{\begin{matrix} L_{u u} u + L_{u ϕ} ϕ \equiv Y u^{I I} - Y ℓ_{mech}^{2} u^{I V} - μ ϕ^{I I I} = - b, \\ L_{ϕ u} u + L_{ϕ ϕ} ϕ \equiv - ε ϕ^{I I} + ε ℓ_{elec}^{2} ϕ^{I V} + μ u^{I I I} = q, \end{matrix}

(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

{\begin{matrix} a_{u u} (δ u, u) + a_{u ϕ} (δ u, ϕ) = l_{u} (δ u), \\ a_{ϕ ϕ} (δ ϕ, ϕ) + a_{u ϕ} (u, δ ϕ) = l_{ϕ} (δ ϕ), \end{matrix}

(B2)

with appropriate boundary conditions, where

a_{u u} (δ u, u) = {(δ u^{I}, Y u^{I})}_{Ω} + {(δ u^{I I}, Y ℓ_{mech}^{2} u^{I I})}_{Ω},

(B3a)

a_{ϕ ϕ} (δ ϕ, ϕ) = {(δ ϕ^{I}, - ε ϕ^{I})}_{Ω} + {(δ ϕ^{I I}, - ε ℓ_{elec}^{2} ϕ^{I I})}_{Ω},

(B3b)

a_{u ϕ} (δ u, ϕ) = - \frac{1}{2} μ ({(δ u^{I}, ϕ^{I I})}_{Ω} - {(δ u^{I I}, ϕ^{I})}_{Ω}),

(B3c)

l_{u} (δ u) = {(δ u, b)}_{Ω},

(B3d)

l_{ϕ} (δ ϕ) = {(δ ϕ, - q)}_{Ω} .

(B3e)

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. to define the following bilinear forms:

\begin{aligned} A_{mech} ({u, ϕ}, {δ u, δ ϕ}) \\ = {(L_{u u} δ u + L_{u ϕ} δ ϕ, τ_{mech} r_{mech} (u, ϕ))}_{Ω}, \end{aligned}

(B4a)

\begin{aligned} A_{elec} ({u, ϕ}, {δ u, δ ϕ}) \\ = {(L_{ϕ u} δ u + L_{ϕ ϕ} δ ϕ, - τ_{elec} r_{elec} (u, ϕ))}_{Ω}, \end{aligned}

(B4b)

with the residuals

r_{mech} (u, ϕ) = L_{u u} u + L_{u ϕ} ϕ + b,

(B5a)

r_{elec} (u, ϕ) = L_{ϕ u} u + L_{ϕ ϕ} ϕ - q .

(B5b)

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

\begin{aligned} a_{u u} (δ u, u) + a_{u ϕ} (δ u, ϕ) + A_{mech} ({u, ϕ}, {δ u, 0}) \\ + A_{elec} ({u, ϕ}, {δ u, 0}) = l_{u} (δ u), \end{aligned}

(B6a)

\begin{aligned} a_{ϕ ϕ} (δ ϕ, ϕ) + a_{u ϕ} (u, δ ϕ) + A_{mech} ({u, ϕ}, {0, δ ϕ}) \\ + A_{elec} ({u, ϕ}, {0, δ ϕ}) = l_{ϕ} (δ ϕ) . \end{aligned}

(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 Codony⁶² as

τ_{mech} = α_{mech} \frac{h^{2}}{Y}, τ_{elec} = α_{elec} \frac{h^{2}}{ε},

(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.

View large Download slide

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} = α$ ⁠.

REFERENCES

1.

W.

Ma

and

L. E.

Cross

, “

Observation of the flexoelectric effect in relaxor Pb (Mg

_{1 / 3}

Nb

_{2 / 3}

⁠)O

_{3}

ceramics

,”

Appl. Phys. Lett.

78

,

2920

–

2921

(

2001

).

https://doi.org/10.1063/1.1356444

Google Scholar

Crossref

2.

L.

Cross

, “

Flexoelectric effects: Charge separation in insulating solids subjected to elastic strain gradients

,”

J. Mater. Sci.

41

,

53

–

63

(

2006

).

https://doi.org/10.1007/s10853-005-5916-6

Google Scholar

Crossref

3.

J.

Bursian

and

O.

Zaikovskii

, “

Changes in curvature of a ferroelectric film due to polarization

,”

Sov. Phys. Solid State

10

,

1121

–

1124

(

1968

).

Google Scholar

4.

É. V.

Bursian

and

N. N.

Trunov

, “

Nonlocal piezoelectric effect

,”

Sov. Phys. Solid State

16

,

760

–

762

(

1974

).

Google Scholar

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

).

https://doi.org/10.1115/1.4031333

Google Scholar

Crossref

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

).

https://doi.org/10.1038/nnano.2015.260

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1038/s41467-019-09266-y

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1063/1.2219990

Google Scholar

Crossref

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

).

https://doi.org/10.1103/PhysRevB.91.104103

Google Scholar

Crossref

10.

S. M.

Kogan

, “

Piezoelectric effect during inhomogeneous deformation and acoustic scattering of carriers in crystals

,”

Sov. Phys. Solid State

5

,

2069

–

2070

(

1964

).

Google Scholar

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

).

Google Scholar

12.

K.

Tolpygo

, “

Long wavelength oscillations of diamond-type crystals including long range forces

,”

Sov. Phys. Solid State

4

,

1297

–

1305

(

1963

).

Google Scholar

13.

A. K.

Tagantsev

, “

Piezoelectricity and flexoelectricity in crystalline dielectrics

,”

Phys. Rev. B

34

,

5883

–

5889

(

1986

).

https://doi.org/10.1103/PhysRevB.34.5883

Google Scholar

Crossref

14.

A. K.

Tagantsev

, “

Electric polarization in crystals and its response to thermal and elastic perturbations

,”

Phase Transit. A

35

,

119

–

203

(

1991

).

https://doi.org/10.1080/01411599108213201

Google Scholar

Crossref

15.

R. D.

Mindlin

, “

Polarization gradient in elastic dielectrics

,”

Int. J. Solids Struct.

4

,

637

–

642

(

1968

).

https://doi.org/10.1016/0020-7683(68)90079-6

Google Scholar

Crossref

16.

E.

Sahin

and

S.

Dost

, “

A strain-gradients theory of elastic dielectrics with spatial dispersion

,”

Int. J. Eng. Sci.

26

,

1231

–

1245

(

1988

).

https://doi.org/10.1016/0020-7225(88)90043-2

Google Scholar

Crossref

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

).

https://doi.org/10.1103/PhysRevB.74.014110

Google Scholar

Crossref

18.

R. D.

Mindlin

and

N. N.

Eshel

, “

On first strain-gradient theories in linear elasticity

,”

Int. J. Solids Struct.

4

,

109

–

124

(

1968

).

https://doi.org/10.1016/0020-7683(68)90036-X

Google Scholar

Crossref

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

).

https://doi.org/10.1016/j.ijsolstr.2011.03.006

Google Scholar

Crossref

20.

R. D.

Mindlin

and

H. F.

Tiersten

, “

Effects of couple-stresses in linear elasticity

,”

Arch. Ration. Mech. Anal.

11

,

415

–

448

(

1962

).

https://doi.org/10.1007/BF00253946

Google Scholar

Crossref

21.

A. R.

Hadjesfandiari

, “

Size-dependent piezoelectricity

,”

Int. J. Solids Struct.

50

,

2781

–

2791

(

2013

).

https://doi.org/10.1016/j.ijsolstr.2013.04.020

Google Scholar

Crossref

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

).

https://doi.org/10.1016/j.jmps.2019.01.013

Google Scholar

Crossref

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

).

https://doi.org/10.1088/0965-0393/24/1/015009

Google Scholar

Crossref

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

).

https://doi.org/10.1088/0022-3727/48/46/465502

Google Scholar

Crossref

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

).

https://doi.org/10.1088/0953-8984/16/13/006

Google Scholar

Crossref

26.

E.

Eliseev

,

A.

Morozovska

,

M.

Glinchuk

, and

R.

Blinc

, “

Spontaneous flexoelectric/flexomagnetic effect in nanoferroics

,”

Phys. Rev. B

79

,

165433

(

2009

).

https://doi.org/10.1103/PhysRevB.79.165433

Google Scholar

Crossref

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

).

https://doi.org/10.1103/PhysRevB.84.174112

Google Scholar

Crossref

28.

L.

Liu

, “

An energy formulation of continuum magneto-electro-elasticity with applications

,”

J. Mech. Phys. Solids

63

,

451

–

480

(

2014

).

https://doi.org/10.1016/j.jmps.2013.08.001

Google Scholar

Crossref

29.

M.-M.

Yang

,

D. J.

Kim

, and

M.

Alexe

, “

Flexo-photovoltaic effect

,”

Science

360

,

904

–

907

(

2018

).

https://doi.org/10.1126/science.aan3256

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1038/s41563-020-0659-y

Google Scholar

Crossref

PubMed

31.

S.

Shen

and

S.

Hu

, “

A theory of flexoelectricity with surface effect for elastic dielectrics

,”

J. Mech. Phys. Solids

58

,

665

–

677

(

2010

).

https://doi.org/10.1016/j.jmps.2010.03.001

Google Scholar

Crossref

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

).

https://doi.org/10.1007/s11433-010-4039-5

Google Scholar

Crossref

33.

P.

Yudin

and

A.

Tagantsev

, “

Fundamentals of flexoelectricity in solids

,”

Nanotechnology

24

,

432001

(

2013

).

https://doi.org/10.1088/0957-4484/24/43/432001

Google Scholar

Crossref

PubMed

34.

T. D.

Nguyen

,

S.

Mao

,

Y.-W.

Yeh

,

P. K.

Purohit

, and

M. C.

McAlpine

, “

Nanoscale flexoelectricity

,”

Adv. Mater.

25

,

946

–

974

(

2013

).

https://doi.org/10.1002/adma.201203852

Google Scholar

Crossref

PubMed

35.

P.

Zubko

,

G.

Catalan

, and

A. K.

Tagantsev

, “

Flexoelectric effect in solids

,”

Annu. Rev. Mat. Res.

24

,

387

–

421

(

2013

).

https://doi.org/10.1146/annurev-matsci-071312-121634

Google Scholar

Crossref

36.

S.

Krichen

and

P.

Sharma

, “

Flexoelectricity: A perspective on an unusual electromechanical coupling

,”

J. Appl. Mech.

83

,

030801

(

2016

).

https://doi.org/10.1115/1.4032378

Google Scholar

Crossref

37.

B.

Wang

,

Y.

Gu

,

S.

Zhang

, and

L.-Q.

Chen

, “

Flexoelectricity in solids: Progress, challenges, and perspectives

,”

Prog. Mater Sci.

106

,

100570

(

2019

).

https://doi.org/10.1016/j.pmatsci.2019.05.003

Google Scholar

Crossref

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

).

https://doi.org/10.1016/j.cma.2016.09.007

Google Scholar

Crossref

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

).

https://doi.org/10.1016/j.cma.2018.05.019

Google Scholar

Crossref

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

).

https://doi.org/10.1016/j.cma.2020.113320

Google Scholar

Crossref

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

).

https://doi.org/10.1016/j.jmps.2020.104182

Google Scholar

Crossref

42.

U. K.

Bhaskar

,

N.

Banerjee

,

A.

Abdollahi

,

E.

Solanas

,

G.

Rijnders

, and

G.

Catalan

, “

Flexoelectric MEMS: Towards an electromechanical strain diode

,”

Nanoscale

8

,

1293

–

1298

(

2016

).

https://doi.org/10.1039/C5NR06514C

Google Scholar

Crossref

PubMed

43.

N.

Mawassy

,

H.

Reda

,

J.-F.

Ganghoffer

,

V. A.

Eremeyev

, and

H.

Lakiss

, “

A variational approach of homogenization of piezoelectric composites towards piezoelectric and flexoelectric effective media

,”

Int. J. Eng. Sci.

158

,

103410

(

2021

).

https://doi.org/10.1016/j.ijengsci.2020.103410

Google Scholar

Crossref

44.

X.

Chen

,

J.

Yvonnet

,

H. S.

Park

, and

S.

Yao

, “

Enhanced converse flexoelectricity in piezoelectric composites by coupling topology optimization with homogenization

,”

J. Appl. Phys.

129

,

245104

(

2021

).

https://doi.org/10.1063/5.0051062

Google Scholar

Crossref

45.

E.

Lifshitz

and

L.

Landau

, Statistical Physics: Course of Theoretical Physics (Butterworth-Heinemann, London, 1984), Vol. 5.

46.

N.

Sharma

,

C.

Landis

, and

P.

Sharma

, “

Piezoelectric thin-film superlattices without using piezoelectric materials

,”

J. Appl. Phys.

108

,

024304

(

2010

).

https://doi.org/10.1063/1.3443404

Google Scholar

Crossref

47.

L. D.

Landau

and

E. M.

Lifshitz

,

Course of Theoretical Physics

(

Elsevier

,

2013

).

Google Scholar

48.

L.

Evans

, Partial Differential Equations, Graduate Studies in Mathematics (American Mathematical Society, 2010).

49.

I.-D.

Ghiba

,

P.

Neff

,

A.

Madeo

, and

I.

Münch

, “

A variant of the linear isotropic indeterminate couple-stress model with symmetric local force-stress, symmetric nonlocal force-stress, symmetric couple-stresses and orthogonal boundary conditions

,”

Math. Mech. Solids

22

,

1221

–

1266

(

2017

).

https://doi.org/10.1177/1081286515625535

Google Scholar

Crossref

50.

S.

Baroudi

,

F.

Najar

, and

A.

Jemai

, “

Static and dynamic analytical coupled field analysis of piezoelectric flexoelectric nanobeams: A strain gradient theory approach

,”

Int. J. Solids Struct.

135

,

110

–

124

(

2018

).

https://doi.org/10.1016/j.ijsolstr.2017.11.014

Google Scholar

Crossref

51.

X.

Liang

,

S.

Hu

, and

S.

Shen

, “

Effects of surface and flexoelectricity on a piezoelectric nanobeam

,”

Smart Mater. Struct.

23

,

035020

(

2014

).

https://doi.org/10.1088/0964-1726/23/3/035020

Google Scholar

Crossref

52.

Q.

Deng

,

M.

Kammoun

,

A.

Erturk

, and

P.

Sharma

, “

Nanoscale flexoelectric energy harvesting

,”

Int. J. Solids Struct.

51

,

3218

–

3225

(

2014

).

https://doi.org/10.1016/j.ijsolstr.2014.05.018

Google Scholar

Crossref

53.

Z.

Yan

and

L.

Jiang

, “

Flexoelectric effect on the electroelastic responses of bending piezoelectric nanobeams

,”

J. Appl. Phys.

113

,

194102

(

2013

).

https://doi.org/10.1063/1.4804949

Google Scholar

Crossref

54.

Z.

Yan

, “

Modeling of a nanoscale flexoelectric energy harvester with surface effects

,”

Phys. E

88

,

125

–

132

(

2017

).

https://doi.org/10.1016/j.physe.2017.01.001

Google Scholar

Crossref

55.

R.

Zhang

,

X.

Liang

, and

S.

Shen

, “

A timoshenko dielectric beam model with flexoelectric effect

,”

Meccanica

51

,

1181

–

1188

(

2016

).

https://doi.org/10.1007/s11012-015-0290-1

Google Scholar

Crossref

56.

R. D.

Mindlin

, “

Micro-structure in linear elasticity

,”

Arch. Ration. Mech. Anal.

16

,

51

–

78

(

1964

).

https://doi.org/10.1007/BF00248490

Google Scholar

Crossref

57.

S.

Mao

and

P. K.

Purohit

, “

Insights into flexoelectric solids from strain-gradient elasticity

,”

J. Appl. Mech.

81

,

1

–

10

(

2014

).

https://doi.org/10.1115/1.4027451

Google Scholar

Crossref

58.

A.

Abdollahi

,

C.

Peco

,

D.

Millán

,

M.

Arroyo

, and

I.

Arias

, “

Computational evaluation of the flexoelectric effect in dielectric solids

,”

J. Appl. Phys.

116

,

093502

(

2014

).

https://doi.org/10.1063/1.4893974

Google Scholar

Crossref

59.

X.

Zhuang

,

B. H.

Nguyen

,

S. S.

Nanthakumar

,

T. Q.

Tran

,

N.

Alajlan

, and

T.

Rabczuk

, “

Computational modeling of flexoelectricity: A review

,”

Energies

13

,

1326

(

2020

).

https://doi.org/10.3390/en13061326

Google Scholar

Crossref

60.

L.

Dorfmann

and

R. W.

Ogden

,

Nonlinear Theory of Electroelastic and Magnetoelastic Interactions

(

Springer

,

2014

), Vol. 1.

Mathematical and computational modeling of flexoelectricity

I. INTRODUCTION

II. CONTINUUM MODELING

A. Direct and converse flexoelectricity

B. Variational principles

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

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

3. Equivalence between different variational principles via partial Legendre transform

C. Boundary value problems

1. Direct flexoelectricity

2. Lifshitz-invariant flexoelectricity

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

4. Beyond standard boundary conditions

a. Sensing electrode condition

b. Interface condition

c. Generalized periodicity condition

D. Flexoelectricity in soft materials: Beyond infinitesimal deformations

III. NUMERICAL MODELING

A. Uniform open B-spline basis

B. Immersed boundary B-spline approximation

4. Cell classification

5. Enforcement of Dirichlet boundary conditions

6. Numerical integration

7. Cut-cell stabilization

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

D. Nondimensionalization of the equations

IV. COMPUTATIONAL EXPERIMENTS OF CANTILEVER BEAM BENDING WITH DIFFERENT FORMULATIONS

A. Cantilever beam bending

B. Cantilever beam actuator

V. COMPUTATIONAL ANALYSIS OF A SOFT CANTILEVER ROD ACTUATOR

VI. COMPUTATIONAL ANALYSIS OF FLEXOELECTRIC DEVICES AND METAMATERIALS

A. Collective-beam bending scalable flexoelectric device

B. Geometrically-polarized periodic inclusions

C. Bending-dominated lattice metamaterial

VII. CONCLUSIONS

ACKNOWLEDGMENTS

DATA AVAILABILITY

APPENDIX A: MATERIAL CHARACTERIZATION

a. Elasticity tensor c

b. Strain gradient elasticity tensor h

c. Dielectricity tensor κ

d. Gradient dielectricity tensor M

e. Flexoelectricity tensor μ

APPENDIX B: RESIDUAL-BASED WEAK FORM STABILIZATION

REFERENCES

2. Optimization of the free enthalpy $\bar{Π} [u, ϕ]$