## Abstract

Flexoelectricity, the linear coupling of strain gradient and electric polarization, is inherently a size-dependent phenomenon. The energy storage function for a flexoelectric material depends not only on polarization and strain, but also strain-gradient. Thus, conventional finite-element methods formulated solely on displacement are inadequate to treat flexoelectric solids since gradients raise the order of the governing differential equations. Here, we introduce a computational framework based on a mixed formulation developed previously by one of the present authors and a colleague. This formulation uses displacement and displacement-gradient as separate variables which are constrained in a ‘weighted integral sense’ to enforce their known relation. We derive a variational formulation for boundary-value problems for piezo- and/or flexoelectric solids. We validate this computational framework against available exact solutions. Our new computational method is applied to more complex problems, including a plate with an elliptical hole, stationary cracks, as well as tension and shear of solids with a repeating unit cell. Our results address several issues of theoretical interest, generate predictions of experimental merit and reveal interesting flexoelectric phenomena with potential for application.

## 1. Introduction

Continuum theories of electro-mechanical phenomena in solids have a long history and have been the subject of several texts including those of Landau *et al.* [1], Maugin & Eringen [2], Kovetz [3] among many others. The classical theory, including that of piezoelectricity, has been successful in analysing and predicting the electromechanical response of solids even in the nonlinear regime. Despite its broad applicability this classical theory does not possess an intrinsic length scale and does not account for gradient effects which are important in certain applications. Mindlin [4], Toupin [5] and Koiter [6] pioneered the study of gradient effects in elastic solids. They incorporated strain gradients into the elastic strain energy function and developed a consistent continuum theory of strain-gradient elasticity (SGE). Later Fleck *et al.* [7,8] extended the theory to strain-gradient plasticity. Various finite-element formulations based on these ideas are also documented in the literature, e.g. [9–13]. Based upon Toupin's variational principles [14], it is possible to generalize the above framework to include size-dependent electromechanical coupling phenomena [15,16], such as flexoelectricity.

Flexoelectricity refers to the linear coupling of strain-gradient and electric polarization. Like other gradient effects, it gives rise to non-local and size-dependent phenomena. Flexoelectricity was first proposed in theory half a century ago [17] and shortly after, discovered in experiments by [18,19]. However, it did not receive much attention within the field of mechanics of solids largely due to limited means of generating large strain gradients. Recently, there has been a revival of research interest in this area, mostly stimulated by advanced fabrication and characterization techniques in nanostructures. Since the gradient scales inversely with the size of the specimen, typically, strain gradients within nanostructures are much greater than their macroscopic counterparts. Experiments have convincingly illustrated the significance of flexoelectricity at the nanoscale and demonstrated its potential to open up novel and unique functionality that cannot be achieved through other means [20–25]. Concurrent theoretical studies have greatly advanced our understanding of the microscopic origins of flexoelectricity. Tagantsev [26,27] developed a framework based on point-charge models and attributed flexoelectricity to lattice effects. Maranganti & Sharma [28] built on this model and calculated flexoelectric constants through lattice dynamics. In addition, recent studies by Hong & Vanderbilt [29,30] revealed that electronic response is also an important aspect of flexoelectricity. We refer the reader to the reviews [31–33] for a detailed description of the state of the art in this field.

There also have been recent developments in the continuum modelling of flexoelectric solids [34–39]. The focus in these studies has been to establish a framework for the solution of useful boundary-value problems of flexoelectric solids that can lead to experiments and applications. For example, Mao & Purohit [39] have presented analytical solutions to several one- and two-dimensional problems and later extended their analysis to determine singular fields around point defects, dislocations and cracks [40]. Finite-element studies have also been conducted by Abdollahi *et al.* [41–44]. They studied several non-trivial geometries, e.g. beam and truncated pyramid structures, which have been extensively used for experimental measurements. Their studies have led to important insights. For example, the non-uniform strain-gradient distribution in a truncated-pyramid around sharp edges significantly influences the measurements of flexoelectric constants. To use the finite-element method with piecewise continuous shape functions to solve problems for flexoelectric solids, one must confront the difficulties arising from higher order differential equations. Abdollahi *et al.* avoided the use of these piecewise continuous functions by applying a mesh-free technique. For two-dimensional problems, they needed to discretize only three degrees of freedom, so their method is computationally efficient. By contrast, our approach still uses the shape function, so it is compatible with the framework of a majority of the current finite-element codes. Our method can be easily incorporated into software packages such as ABAQUS. Therefore, it can be used by non-expert engineers for the analysis of complex geometries.

This paper introduces a general framework for finite-element solutions of problems for an elastic dielectric with flexoelectricity and/or piezoelectricity. The generalized gradient theory developed by Mindlin [4] is used to model the gradient effect of elasticity. Piezoelectric as well as flexoelectric coupling are introduced into the formulation by adding polarization as a variable in the energy storage function. The energy storage function depends on the strain tensor, second gradient of displacement and polarization. To avoid using *C*^{1} finite elements in our numerical solution, a mixed formulation based on the work of Amanatidou & Aravas [13] is developed. In this formulation, displacement and displacement gradient are treated as separate degrees of freedom and their relationship is enforced in the variational form. This framework is entirely consistent with the continuum theory of flexoelectricity and is capable of capturing fine structures due to gradient effects. The finite-element code is validated against benchmark problems with known analytical solutions. Then it is employed to study three important classes of problems: plate with an elliptical hole, stationary crack and periodic meta-structures. In the stationary crack problem, for which an asymptotic solution has been developed by Mao & Purohit [40], the validity and region of dominance of the asymptotics is determined. The elliptical hole and periodic structure provide an alternative means of generating large strain gradients; the finite-element results show how these large gradients influence classical observations and generate crucial insights that can lead to better measurement in experiments as well as improved functionality in applications.

Standard notation is used throughout. Tensors, including vectors, are denoted by boldface symbols, whose orders are indicated by the context. Einstein summation convention is used here for repeated Latin indices of tensor components with respect to a fixed Cartesian coordinate system. A comma followed by a subscript, say *i*, denotes partial differentiation with respect to the corresponding spatial coordinate *x*_{i}, i.e. *f*_{,i}=∂*f*/∂*x*_{i}. Let **a** and **b** be vectors and **A** a second-order tensor; the following products are used in the text: **a**⋅**b**=*a*_{i}*b*_{i}, (**A**⋅**a**)_{i}=*A*_{ik}*a*_{k}, and (**a**⋅**A**)_{i}=*a*_{k}*A*_{ki}.

## 2. Constitutive model and boundary-value problem

Consider a homogeneous elastic dielectric body and introduce a fixed Cartesian coordinate system with an orthonormal basis (**e**_{1},**e**_{2},**e**_{3}) and spatial coordinates given as (*x*_{1},*x*_{2},*x*_{3}). The body occupies a region *Ω* in a fixed reference configuration with boundary ∂*Ω* and outward unit normal vector **n**.

In response to mechanical and electrical loads, the body deforms and polarizes. The mechanical response of the material is described by the displacement vector field **u**(*x*_{1},*x*_{2},*x*_{3}), and the electric response is characterized by a polarization vector field **P**(*x*_{1},*x*_{2},*x*_{3}), which is related to bounded or polarized charge in the dielectric.

Infinitesimal displacement gradients are assumed, hence the strain tensor is

Our constitutive model is based on an energy function per unit volume ** ε** and the second gradient of displacement

*σ*^{(0)}, the double-stress

**E**are

*ϕ*is the electric potential,

**b**the body force per volume,

*q*the free charge per volume and

*ϵ*

_{0}is the permittivity of free space. Equation (2.4) represents the conservation of charge. The corresponding boundary conditions are

^{n}=

**n**⋅

**∇**=

*n*

_{i}(∂/∂

*x*

_{i}) is the normal derivative,

**∇**

^{t}=

**∇**−

**n**∇

^{n}the ‘surface gradient’ on ∂

*Ω*, ∂

*Ω*

_{u}∪∂

*Ω*

_{Q}=∂

*Ω*

_{d}∪∂

*Ω*

_{R}=∂

*Ω*

_{ϕ}∪∂

*Ω*

_{ω}=∂

*Ω*, and

*C*

^{β}, and

**ℓ**=

**s**×

**n**, where

**s**is the unit vector tangent to

*C*

^{β}.

## 3. Variational formulation

Following the works of Amanatidou & Aravas [13] and Yang & Batra [45,46], we can show that the boundary-value problem defined in §2 can be formulated alternatively by the stationarity condition *δΠ*=0 of the functional
*α*^{t}=** α**−(

**⋅**

*α***n**)

**n**is the ‘tangential part’ of

**on ∂**

*α**Ω*, with

*δ*

**u**=

**0**on ∂

*Ω*

_{u}and

*δ*

**⋅**

*α***n**=

**0**on ∂

*Ω*

_{d}, and

*δϕ*=0 on ∂

*Ω*

_{ϕ}. Using the ‘surface divergence theorem’

*δΠ*=0 implies the field equations

In the above functional, the quantities *Ω* and on ∂*Ω*.

## 4. ‘Mixed’ finite-element formulation

Functional (3.1) forms the basis for a ‘mixed’ finite-element formulation, in which **u**, ** α**,

*ϕ*and

**P**are the nodal variables. The stationarity condition

*δΠ*leads to

*δ*

**u**=

**0**on ∂

*Ω*

_{u}and

*δ*

**⋅**

*α***n**=

**0**on ∂

*Ω*

_{d}, and

*δϕ*=0 on ∂

*Ω*

_{ϕ}.

The finite-element solutions are based on (4.1). We develop the 9-node isoparametric plane-strain element (I9-87) shown in figure 1. The quantities (*u*_{1},*u*_{2},*α*_{11},*α*_{22},*α*_{21},*α*_{12},*ϕ*) are used as degrees of freedom at all nodes; the quantities *u*_{1},*u*_{2},*α*_{11},*α*_{22},*α*_{21},*α*_{12},*ϕ*) and a bi-linear interpolation for

The element described above is implemented into the ABAQUS general purpose finite-element program [47]. This code provides a general interface so that a particular new element can be introduced as a ‘user subroutine’ (UEL).

The formulation described by the functional (3.1) is valid for materials with energy function of a general form, including those with nonlinear constitutive laws. Here we focus attention on linear materials with a general energy function *A*_{ijkpqr},*a*_{ij},*d*_{ijk},*f*_{ijkm}) are constitutive tensors. In the problems, we use isotropic materials with an energy function *μ*) are the usual Lamé parameters, ℓ is an internal ‘material length’, *a* is reciprocal susceptibility constant, which is related to the permittivity of the dielectric *ϵ* by *ϵ*=*ϵ*_{0}+*a*^{−1}. Constants *f*_{1},*f*_{2} are the two flexoelectric constants and we often refer to *f*=*f*_{1}+2*f*_{2} as the volumetric flexoelectric constant. Note that the third-order piezoelectric tensor *d*_{ijk} on the right-hand side of (4.3) vanishes for materials with centrosymmetry, e.g. isotropic or cubic materials. The corresponding constitutive equations are
*δ*_{ij} is the Kronecker delta. Note that, when **P**=**0** the energy function can be written also in the well-known form [13]

## 5. Applications

### (a) Code validation

The element I9-87 passes the patch test of bi-quadratic displacement field under pure gradient elasticity (all electric nodal degrees of freedom suppressed, i.e. *ϕ*=0,*P*_{i}=0) and bi-linear potential field for pure electrostatics (all displacement and stress nodal degrees of freedom suppressed, i.e. *a*).

The tube is loaded by internal and external pressures *p*_{i} and *p*_{o}, respectively, and a voltage difference *V* is applied across the inner and outer surfaces. The corresponding boundary conditions are
*A*,*B*,*C*,*D*,*G*,*H*) are constants determined from the boundary conditions, *f*=*f*_{1}+2*f*_{2}, *aϵ*=1+*aϵ*_{0}, *I*_{1}(*x*) and *K*_{1}(*x*) the first-order modified Bessel functions of the first and second kind, respectively, and

Calculations are carried out for the following non-dimensional parameters:
*E* is Young's modulus and *ν* Poisson's ratio with which we can recover the Lamé parameters. In the view of the axial symmetry, the problem is mathematically one-dimensional, since the solution depends only on the radial coordinate *r*. In the finite-element calculations, one-quarter of the cross section was analysed and appropriate symmetry conditions were enforced. Figure 2*b* shows the 40×20 finite-element mesh used in the calculations.

Figure 3*a* shows a comparison of the numerical and analytical solutions for the SGE (*f*_{1}=*f*_{2}=0) and the flexoelectric (coupled) problems. In both cases, the numerical solutions agree very well with the corresponding analytical solutions. For a mixed formulation, the ‘Lagrange multiplier’ fields are of interest due to potential instability. Therefore, a comparison of components of *b*. The finite-element solution exhibits good agreement and stability.

Figure 4 shows the variation of the electric potential *ϕ* and the polarization *P*_{r} as determined from the finite-element solution together with the analytical solution; again there is excellent agreement between the numerical and analytical solutions.

### (b) Elliptical hole in a plate

We consider the problem of a cylindrical elliptical hole in a plane strain tension field and a uniform electric field (figure 5). The major axis of the ellipse is in the horizontal *x*_{1}-direction and the tension field *x*_{2}-direction. The electric field is created by the opposite surface charges *r*_{a}/*r*_{b}=2. The surface of the hole is assumed to be traction- and charge-free.

The boundary conditions of the problem are

where *f*=*f*_{1}+2*f*_{2} and *aϵ*=1+*aϵ*_{0}. The boundary conditions listed above are consistent with uniform stress and electric fields at infinity.

A square plate with dimensions 2*w*×2*w*, with *w*=10*r*_{a}, is used in the finite-element calculations; because of symmetry, one half of the plate is analysed and the appropriate symmetry conditions are imposed (figure 5*b*). The side *w* is substantially larger than *r*_{a} and the solution of this finite-size plate problem is expected to be close to the infinite domain problem. The calculations are carried out for
*ε*_{22} and the polarization *P*_{2} along the *x*_{1}-axis ahead of the elliptical hole. A concentration of strain and polarization appears at the ‘tip’ of the hole over a distance approximately equal to the size *r*_{a} of the major semi-axis.

Figure 6*a* shows also the corresponding results of strain-gradient elasticity without any flexoelectric coupling, i.e. for *f*_{1}=*f*_{2}=0. Figure 6*b* shows the polarization *P*_{2} ahead of the hole as determined by pure electrostatics as well. For the values of the parameters used in the calculations, it appears that the flexoelectric effects have minimal influence on deformation field along *x*_{1}-axis but greater influence on the polarization field.

Figure 7 shows contour plots of the normal strain *ε*_{22} and polarization *P*_{2} in the plate. Owing to the flexoelectric coupling, the profiles are not symmetric with respect to *x*_{1}-axis, in spite of the centrosymmetric geometry. It is also interesting to note that this effect which breaks the symmetry of polarization field depends only on how the material is ‘poled’, not ‘material anisotropy’, because our constitutive equations are isotropic. If we flip the electric field, then the net polarization is rotated by 180°. This is similar to poling of certain piezoelectrics [44]. In our simplified material model this effect is reversible; in real materials, however, the stress and electric fields might cause elliptical (or other) defects to move or migrate (evolution of microstructure) in an irreversible manner and even create residual polarization.

### (c) Stationary crack

We consider the plane strain problem of an edge-cracked panel (ECP) loaded with a uniformly distributed load as shown in figure 8*a* or by a uniform far field electric load resulting in surface charge (figure 8*b*). The crack faces are assumed to be traction- and charge-free. This is an ideal model of insulating crack, also called the impermeable condition [50].

In the following, we use the finite-element solution to determine the coefficients (stress intensity factors) that enter the asymptotic solution developed by Mao & Purohit [40]. The panel that we studied here is a block with total width *w* and total height of 2*h*. The edge crack is placed at the left half of the specimen (starting from origin), with a length of *w*/2, as shown in figure 8.

Relevant non-dimensional parameters are

First, consider a Mode I insulating crack loaded by a mechanical load *T*. The boundary conditions for this problem are (figure 8*a*)

The asymptotic crack-tip fields in a flexoelectric solid were developed recently by Mao & Purohit [40]; these crack-tip fields are different from the corresponding fields in ‘linear elastic fracture mechanics’ (LEFM) and in ‘linear piezoelectric fracture mechanics’ (LPFM). In particular, it was shown that the leading term in the asymptotic expansion of the crack-tip displacement field is *r*^{3/2}, *r* being the radial distance from the crack-tip.

Using the notation of Mao & Purohit [40], we consider the crack-tip intensities
*Ω*_{3}(*r*,*θ*)=(*u*_{2,1}−*u*_{1,2})/2 is the out of plane component of the infinitesimal rotation vector ** Ω**. We write the solution in the form (see also [51])

*u*

_{2},

*Ω*

_{3}) on the crack face (

*θ*=

*π*) together with the predictions of the asymptotic solution. The leading term provides an accurate description of the displacement and rotation fields on

*θ*=

*π*in the range 0<

*r*<ℓ/10.

The asymptotic solution for the polarization field is [40]
*a* shows the variation of the polarization *P*_{r} ahead of the crack (*θ*=0) together with predictions of the asymptotic solution. The region of ‘C-dominance’ (in analogy to K-dominance in LEFM), which is the region where the asymptotic field dominates other terms, is about ℓ/10. Figure 10*b* shows the angular variation of *P*_{r} at radial distances from the tip *r*=ℓ/10,ℓ/15, ℓ/20 as determined numerically and as predicted by the asymptotic solution. It appears that the leading term of the asymptotic solution is very accurate ahead of the crack tip for values of *θ* in the range of 0 and about 120°; closer to the crack face, i.e. for

Another type of crack that uses the impermeable condition is the ‘Mode D’ crack shown in figure 8*b*. In a Mode D crack, there is no mechanical loading and an electric field is applied perpendicular to the direction of the crack faces so that charges of equal magnitude and opposite sign are induced on the top and bottom surfaces of the specimen. The corresponding boundary conditions are

The electric intensity factor for this type of crack which is given by
*Ω*_{3} are [40]

where *D* is a constant. The value of the constant *D* is determined from the numerical solution to be
*Ω*_{3} ahead of the crack together with the prediction of the two-term asymptotic solution (5.25).

We conclude this section by mentioning that the crack tip strains are non-singular, whereas the polarization field has an *r*^{−1/2} singularity at the crack tip. The region of dominance is quite small (≃ℓ/10) since the gradient characteristic length scale ℓ is usually in the range 10–100 nm [16]. However, the intensities of the asymptotic solution determine the energy release rate, which in turn determines the conditions for crack growth initiation.

### (d) Periodic structures

Recently, new material processing techniques have been used to produce solids with periodic structures to create meta-materials for improved or desired functionality. For instance, nanomeshes, a periodic array of squares with a circular hole inside can be fabricated to achieve higher thermoelectric responses of crystal silicon [52]. An example of such a structure is sketched in figure 12.

In order to determine the macroscopic electromechanical response of this periodic structure, we study the behaviour of a square unit cell ABCD using appropriate periodic boundary conditions as described in the following.

Let **u**(**x**) in the unit cell can be written in the form
**u***(**x**) is a periodic function with zero mean deformation gradient on the unit cell.

We denote the quantities on sides AB, BC, CD and DA with superscript l, b, r and u, respectively. Periodicity requires that
*L* be the length of the sides of the square unit cell. Then

We use the finite-element method to study the response of the square unit cell when subject to mechanical and electrical loads. Since the displacement gradients *α*_{ij}=*u*_{i,j} are treated as independent degrees of freedom in the finite-element formulation, similar periodicity conditions are imposed on ** α** in the numerical solution:

We consider first the case in which the macroscopic loads on the unit cell are a normal strain *x*_{2}-direction; the electric field is created by opposite charges *R* is the radius of the hole and *L* the length of the sides of the unit cell. This creates a meta-material with defect volume fraction of *π*/16≃19.6%.

Figure 13 shows the variation of the normal strain *ε*_{22} and the polarization *P*_{2} along the *x*_{1}-axis ahead of the hole. The corresponding solutions of SGE and electrostatics are also included in the same figure. The strain distribution *ε*_{22} appears to be insensitive to the values of the flexoelectric constants, whereas the polarization changes significantly when these constants are varied.

We consider next the problem in which the macroscopic loads on the square unit cell are a shear strain

Figure 14 shows the variation of the shear strain *ε*_{12} and the polarization *P*_{1} along the *x*_{2}-axis above the hole. The corresponding solutions of SGE and electrostatics are also included in the same figure. The variation of *ε*_{12} is affected by flexoelectricity and the relative effect is related to the magnitude of the flexoelectric constant. This effect is highly localized near the hole or meta-defect, whereas the overall profile of *ε*_{12} is relatively flat and a small gradient develops along the *x*_{2}-axis. It is very interesting to note that a polarization is produced in the *x*_{1}-direction due to the coupling of strain gradient and polarization, i.e. flexoelectricty rotates the polarization field towards the *x*_{1}-axis. Similar polarization rotation phenomena have also been reported in the literature [23,24]; they are realized in ferroelectric thin films and are believed to have applicability in memory devices.

The analyses of the periodic meta-structure and the elliptical hole problems suggest an alternative way of studying flexoelectricity. The classical solution of a circular hole in an infinite elastic body (under uniaxial tension) predicts a stress concentration factor of 3 and that strain/stress decays to the far field level *r*/*R*)^{−2} [53]. A good estimate for the strain gradient around the hole, where *r*∼*R*, can be
*η* is the concentration factor. Therefore, when the radii *R* of the holes are small, these periodic structures can generate considerably large strain gradients. For the same macroscopic deformation level, a reduction of the size of the hole produces larger strain gradients. In fact, periodic nano-scale or even atomistic scale holes have been observed to alter the electromechanical behaviours of certain two-dimensional materials [54]. However, holes of these scales are difficult to make. On the other hand, for meta-materials, the size of the hole can be in the range of hundreds of nanometres [52], which, by the above analysis, can also produce large gradients. For these structures, we can design the arrangement, size and spacing so as to meet different needs.

Moreover, this periodic structure could provide an alternative method to measure flexoelectric constants. So far, the most reliable means to measure them is through beam bending experiments. These, however, cannot determine all components of the flexoelectric tensor (even for the simplest isotropic case) [22]. For example, a recent study [30] predicted that some materials (such as, silicon) have finite volumetric flexoelectric constant *f*, but vanishing or very small bending flexoelectric constants. Beam bending experiments are not expected to be useful in determining the flexoelectric constants for these materials. Therefore, alternative measurement techniques are required. The truncated pyramid structure provides an alternative, but non-trivial deformation concentration around the pyramid edges makes it difficult to use [43]. The periodic structure studied here could be an ideal set up to overcome these difficulties. It gives a large gradient without singular fields due to sharp edges. The magnitude of the gradients can be easily controlled by altering loading or geometry (without exceeding the elastic limit). Both *f* and *f*_{1} (isotropic case) appear in the solutions, so we can determine them by appropriate measurements.

## 6. Conclusion

In this paper, we have formulated a variational form that is completely consistent with the continuum theory of flexoelectricity. The form uses a mixed formulation and circumvents the difficulties of modelling gradient effects in flexoelectric solids by introducing extra degrees of freedom. This variational form is general and can incorporate the piezoelectric effect as well. A new element is developed for adapting the variational form to finite-element calculation. The known analytic solution of a pressurized tube is employed as a benchmark problem for validation. Then the method is used to study three types of problems which are beyond current analytic capabilities. Asymptotic theories of cracks are confirmed and a more precise description of the fracture landscape is accomplished. Single hole in an infinite medium as well as periodic meta-structures illustrate the non-trivial coupling of electric loading and deformation. They also provide insights for alternative means of measuring and using flexoelectricity.

## Data accessibility

There is no data associated with this work.

## Authors' contributions

All authors contributed substantially to the work. In particular, the derivation of the formulation and the implementation of the code was done by S.M. and N.A., the particular problems were chosen and the paper was written by S.M., P.K.P. and N.A.

## Competing interests

We declare no competing interests.

## Funding

S.M. and P.K.P. acknowledge partial support for this work through the Army Research Office grant no. W911NF-11-1-0494.

## Acknowledgements

We thank Mr G. Tsantidis for providing code and useful discussions. Prof. K. Turner for providing access to ABAQUS v. 6.9 EF software package.

- Received December 24, 2015.
- Accepted May 12, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.