An Immersed Boundary Hierarchical B-spline method for flexoelectricity

TL;DR

This paper introduces a novel immersed boundary hierarchical B-spline computational framework for flexoelectricity and piezoelectricity problems, enabling accurate simulations on complex geometries with adaptive refinement.

Contribution

It develops an unfitted mesh method using hierarchical B-splines for high-order PDEs in flexoelectricity, incorporating NURBS boundary representation and Nitsche's method for boundary conditions.

Findings

Achieves optimal error convergence rates with high-order B-splines.

Successfully verified against analytical solutions and benchmarks.

Handles complex geometries with adaptive local refinement.

Abstract

This paper develops a computational framework with unfitted meshes to solve linear piezoelectricity and flexoelectricity electromechanical boundary value problems including strain gradient elasticity at infinitesimal strains. The high-order nature of the coupled PDE system is addressed by a sufficiently smooth hierarchical B-spline approximation on a background Cartesian mesh. The domain of interest is embedded into the background mesh and discretized in an unfitted fashion. The immersed boundary approach allows us to use B-splines on arbitrary domain shapes, regardless of their geometrical complexity, and could be directly extended, for instance, to shape and topology optimization. The domain boundary is represented by NURBS, and exactly integrated by means of the NEFEM mapping. Local adaptivity is achieved by hierarchical refinement of B-spline basis, which are efficiently evaluated and integrated thanks to their piecewise polynomial definition. Nitsche's formulation is derived to weakly enforce essential boundary conditions, accounting also for the non-local conditions on the non-smooth portions of the domain boundary (i.e. edges in 3D or corners in 2D) arising from Mindlin's strain gradient elasticity theory. Boundary conditions modeling sensing electrodes are formulated and enforced following the same approach. Optimal error convergence rates are reported using high-order B-spline approximations. The method is verified against available analytical solutions and well-known benchmarks from the literature.

Figures

Equations

$$p_l=d_{lij}{\epsilon }_{ij},$$

$$p_l=f_{lijk}\frac{\partial {\epsilon }_{ij}}{\partial x_k},$$

$${\sigma }_{ij}=f_{lijk}\frac{\partial E_l}{\partial x_k}.$$

$${\left[{\nabla }^{\text{sym}}(\mathbf{a})\right]}_{ij}:=\frac{1}{2}{\left[\nabla (\mathbf{a})+\nabla {(\mathbf{a})}^T\right]}_{ij}=\frac{1}{2}(a_{i,j}+a_{j,i}).$$

$${\left[{\nabla }^S(\mathbf{A})\right]}_{ijk}={\nabla }_k^S{A}_{ij}:=$$

$${\left[{\nabla }^S\cdot (\mathbf{A})\right]}_i={\nabla }_j^S{A}_{ij}:=$$

$$[\mathbf{S}{]}_{ij}:=-{\left[{\nabla }^S(\mathbf{n})\right]}_{ij}=-n_{i,l}{P}_{lj}.$$

$$H:=\frac{1}{2}Tr(\mathbf{S})=-\frac{1}{2}{\nabla }^S\cdot \mathbf{n}=-\frac{1}{2}{n}_{i,j}{P}_{ij}.$$

$${\left[\tilde{\mathbf{N}}\right]}_{ij}:=S_{ij}-2H{n}_i{n}_j.$$

$${\delta_{f_{i}}F}[f_{1},\dots,f_{n};\psi_{i}]\coloneqq\frac{d}{d\epsilon}F[f_{1},\dots,f_{i-1},f_{i}+\epsilon\psi_{i},f_{i+1},\dots,f_{n}]\Big{|}_{\epsilon=0},$$

$$\delta F[f_1,\dots ,f_n;{\psi }_1,\dots ,{\psi }_n]:=\sum _{i=1}^n{\delta }_{f_i}F[f_1,\dots ,f_n;{\psi }_i].$$

$${\delta^{2}_{f_{i}}F}[f_{1},\dots,f_{n};\psi_{i}]\coloneqq{\delta_{f_{i}}\Big{(}{\delta_{f_{i}}F}[f_{1},\dots,f_{n};\psi_{i}]\Big{)}}[f_{1},\dots,f_{n};\psi_{i}],$$

$$[\mathbf{\epsilon }(\mathbf{u}){]}_{ij}=[\mathbf{\epsilon }(\mathbf{u}){]}_{ji}:={\left[{\nabla }^{\text{sym}}(\mathbf{u})\right]}_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i}),$$

$$[\mathbf{E}(\varphi ){]}_l:=-[\nabla \varphi {]}_l=-{\varphi }_{,l}.$$

$$\begin{gathered} {\mathcal{H}}^{\Omega }[\mathbf{u},\varphi ]={\mathcal{H}}^{\Omega }[\mathbf{\epsilon },\nabla \mathbf{\epsilon },\mathbf{E}]:=\frac{1}{2}{\epsilon }_{ij}\text{varmathbb}{C}_{ijkl}{\epsilon }_{kl}+\frac{1}{2}{\epsilon }_{ij,k}{h}_{ijklmn}{\epsilon }_{lm,n} \\ -\frac{1}{2}{E}_l{\kappa }_{lm}{E}_m-E_l{e}_{lij}{\epsilon }_{ij}-E_l{\mu }_{lijk}{\epsilon }_{ij,k}. \end{gathered}$$

$${\mathcal{W}}^{\Omega }[\mathbf{u},\varphi ]:=-b_i{u}_i+q\varphi .$$

$${\mathcal{W}}_0^{\partial \Omega }[\mathbf{u},\varphi ]:=-t_i{u}_i-r_i{\partial }^n{u}_i+w\varphi .$$

$${\mathcal{W}}_0^C[\mathbf{u},\varphi ]:=-j_i{u}_i.$$

$$\Pi_{0}[{{\boldsymbol{u}}},\phi]=\int_{\Omega}\Big{(}\mathcal{H}^{\Omega}[{{\boldsymbol{u}}},\phi]+\mathcal{W}^{\Omega}[{{\boldsymbol{u}}},\phi]\Big{)}\text{\,}\mathrm{d}\Omega+\int_{\partial\Omega}\mathcal{W}^{\partial\Omega}_{0}[{{\boldsymbol{u}}},\phi]\text{\,}\mathrm{d}\Gamma+\int_{C}\mathcal{W}^{C}_{0}[{{\boldsymbol{u}}},\phi]\text{\,}\mathrm{d}\mathrm{s}.$$

$$\mathbf{u}-\bar{\mathbf{u}}$$

$${\partial }^n(\mathbf{u})-\overline{\mathbf{v}}$$

$$\varphi -\bar{\varphi }$$

$$\mathbf{u}-\bar{\mathbf{u}}$$

$${\Pi }_0[\mathbf{u},\varphi ]={\Pi }^{\Omega }[\mathbf{u},\varphi ]+{\Pi }_0^{\text{Dirichlet}}[\mathbf{u},\varphi ]+{\Pi }^{\text{Neumann}}[\mathbf{u},\varphi ],$$

$$\Pi^{\Omega}[{{\boldsymbol{u}}},\phi]=\int_{\Omega}\Big{(}\mathcal{H}^{\Omega}[{{\boldsymbol{u}}},\phi]+\mathcal{W}^{\Omega}[{{\boldsymbol{u}}},\phi]\Big{)}\text{\,}\mathrm{d}\Omega,$$

$$\begin{gathered} {\Pi }_0^{\text{Dirichlet}}[\mathbf{u},\varphi ]={\int }_{\partial {\Omega }_u}-u_i{t}_i(\mathbf{u},\varphi )\mathrm{d}\Gamma +{\int }_{\partial {\Omega }_v}-{\partial }^n{u}_i{r}_i(\mathbf{u},\varphi )\mathrm{d}\Gamma \\ +{\int }_{\partial {\Omega }_{\varphi }}\varphi w(\mathbf{u},\varphi )\mathrm{d}\Gamma +{\int }_{C_u}-u_i{j}_i(\mathbf{u},\varphi )\mathrm{ds} \end{gathered}$$

$${\Pi }^{\text{Neumann}}[\mathbf{u},\varphi ]={\int }_{\partial {\Omega }_t}-u_i{\bar{t}}_i\mathrm{d}\Gamma +{\int }_{\partial {\Omega }_r}-{\partial }^n{u}_i{\bar{r}}_i\mathrm{d}\Gamma +{\int }_{\partial {\Omega }_w}\varphi \bar{w}\mathrm{d}\Gamma +{\int }_{C_j}-u_i{\bar{j}}_i\mathrm{ds}.$$

$${\mathcal{U}}_0$$

$${\mathcal{P}}_0$$

$$({\mathbf{u}}^{\text{*}},{\varphi }^{\text{*}})=\arg \underset{\mathbf{u}\in {\mathcal{U}}_0}{\min }\underset{\varphi \in {\mathcal{P}}_0}{\max }{\Pi }_0[\mathbf{u},\varphi ].$$

$$\delta {\Pi }_0={\delta }_{\mathbf{u}}{\Pi }_0+{\delta }_{\varphi }{\Pi }_0=0;$$

$${\delta }_{\mathbf{u}}^2{\Pi }_0>0,{\delta }_{\varphi }^2{\Pi }_0<0,$$

$$\{\begin{array}{ll}{\left({\hat{\sigma }}_{ij}(\mathbf{u},\varphi )-{\tilde{\sigma }}_{ijk,k}(\mathbf{u},\varphi )\right)}_{,j}+b_i=0_i& \text{in }\Omega ,\\ \text{hfill}{\hat{D}}_{l,l}(\mathbf{u},\varphi )-q=0& \text{in }\Omega ;\end{array}$$

Full text

*Keywords: * Flexoelectricity, Piezoelectricity, Strain gradient elasticity , Immersed boundary B-spline approximation , High-order PDE , Nitsche’s method

1 Introduction

Electroactive materials are able to transform mechanical energy into electrical energy (and viceversa), which can be used for sensing, actuating or energy harvesting applications. A wide range of modern technologies are based on the electromechanical properties of these materials, such as cameras, printers or motors.

Different electromechanical couplings can be found depending on the material. The most common coupling is piezoelectricity, by which the strain ${\boldsymbol{\varepsilon}}$ and polarization ${\boldsymbol{p}}$ are linearly coupled:

[TABLE]

where $\mathbf{d}$ is the third-rank tensor of piezoelectricity. This is the case of piezoelectric ceramics, which are polarized by deformation, and conversely deform when an electrical field is applied. Some piezoelectrics exhibit further electromechanical couplings, such as pyroelectricity (temperature-dependent polarization) or ferroelectricity (reversible spontaneous polarization). Soft materials such as piezoelectric polymers or dielectric elastomers exhibit also electrostriction, a nonlinear electromechanical coupling between the strain state and the square of the polarization field.

This variety of electromechanical couplings has been largely studied, is quite well understood and is suitable to model electromechanical couplings in materials at a macroscale. However, micro- and nanoscale electromechanics cannot be described by just considering traditional models, because additional effects become relevant at small scales, prominently flexoelectricity.

Flexoelectricity is a two-way linear coupling between electric polarization and strain gradient. The (direct) flexoelectric effect is understood as the material polarization due to inhomogeneous deformation (e.g. bending) and is mathematically expressed as

[TABLE]

where $\mathbf{f}$ is the fourth-rank tensor of flexoelectricity. There also exists a thermodynamically conjugate converse flexoelectric effect that consists on the generation of stress ${\boldsymbol{\sigma}}$ due to the application of an inhomogeneous electric field ${\boldsymbol{E}}$, i.e.

[TABLE]

Compared to piezoelectricity, flexoelectricity has two distinctive features. On the one hand, it is universal, meaning it is present in any dielectric material. For a crystalline material to be piezoelectric, its crystalline structure is required to be non-centrosymmetric in order to allow for a net polarization as a result of a uniform deformation (Fig. 4). Otherwise, the relative position of positive and negative ions remains unchanged after deformation and no net polarization is expected (Fig. 4). However, flexoelectricity generically breaks the inversion symmetry of the material, regardless of the internal crystalline structure, and a net polarization is observed after a non-uniform mechanical stimulus such as bending (Fig. 4). On the other hand, the flexoelectric material constants are typically small, and therefore sufficiently large strain gradients are required in order to trigger a sizable flexoelectric effect. Since strain-gradients scale inversely to spatial dimension, they are considerably large in the micro- and nanoscale. Therefore, flexoelectricity is by nature a size dependent effect.

Flexoelectricity in crystalline dielectrics was first studied by Mashkevich [1], Tolpygo [2] and Kogan [3], who proposed the first phenomenological model. In 1968, Bursian et. al. [4] performed the first experiment showing evidence of flexoelectricity in ferroelectric films and, in fact, it is not until 1981 that the phenomenon is named flexoelectricity [5]. The first comprehensive theoretical works by Tagantsev [6, 7] clarified the distinction between piezoelectricity and flexoelectricity. However, since its effect is negligible at the macro-scale, it received little attention. In recent years, there has been a renewed interest in the scientific community, motivated by the need to downscale electromechanical transduction, and enabled by recent developments in nanotechnology [8]. Maranganti proposed the first mathematical framework for the flexoelectric governing equations [9]. Following this work, flexoelectricity has been studied analytically for simple reduced models under restrictive assumptions, such as cantilever beams [10] and thin films [11], to name a few. In this last work by Sharma, atomistic calculations are also performed in order to verify the analytical results. Cross developed a formulation to measure experimentally the longitudinal flexoelectric effect [12]. Other authors consider further physics, such as the flexoelectric effect in ferroelectrics [13, 14], the coupling with magnetic fields [15] and the contributions of surface effects [16]. The general variational principles for flexoelectric materials can be found in [17, 16, 15]. The reader is referred to [18, 8, 19, 20] for recent reviews of flexoelectricity in solids.

Within the continuum flexoelectric theory, the symmetry of the flexoelectric tensor is well understood [21, 22], although its full characterization is still lacking for most materials [19]. The equations are a coupled system of 4th-order partial differential equations, which renders analytical solutions difficult to obtain and precludes the use of conventional $C^{0}$ finite elements. Several numerical alternatives have been proposed in the literature, based on smooth approximations with at least $C^{1}$ continuity [23, 24, 25, 26, 27, 28, 29] or on mixed formulations [30, 31]. The first self-consistent numerical solution of the linear flexoelectric problem was provided by Abdollahi et. al. [23, 24, 25] using a mesh-free approach in 3D. The degrees of freedom correspond only to displacements and electric potential, discretized with a $C^{\infty}$-continuous approximation to address the high-order nature of the equations. This method was successfully applied to study the effect of flexoelectricity on the fracture of piezoelectric materials [25], and on the design of bimorph microsensors and microactuators [32]. Later, an alternative 2D continuum approach was proposed by the group of Aravas et. al. [30], extending the mixed FEM formulation originally developed in [33] for strain-gradient elasticity to flexoelectricity. Displacement and displacement gradient fields are treated as separate degrees of freedom in order to circumvent the $C^{1}$-continuity requirement. This approach was also used by Deng [31]. Another alternative is the isogeometric approach, which has been used to perform topology optimization on 2D flexoelectric cantilever beams [26, 27, 28]. More recently, the $C^{1}$ triangular Argyris element was used by Yvonnet et. al. in [29] to model flexoelectricity in soft dielectrics at finite strains.

In this paper, we propose an immersed boundary hierarchical B-spline approach to numerically solve the governing equations of flexoelectricity in 2D and 3D. This method enables simulations on arbitrary geometries within a reasonable computational cost, unlike previous works in the literature. For the sake of simplicity we restrict ourselves to infinitesimal strains, although the same idea applies also to finite strains [34].

In this approach, the domain boundary is immersed into a fixed Cartesian mesh, and a hierarchical B-spline basis is built on top of it to discretize the primal unknowns (i.e. displacements and electric potential), fulfilling the smoothness requirement of the equations. The computational mesh does not fit to the embedded boundary, overcoming the rigidity of IGA approaches [26] that require Cartesian-like body-fitted meshes, difficult to generate for non-trivial geometries. In our case, mesh generation is straightforward regardless of the complexity of the domain shape. Moreover, a fixed mesh facilitates shape and topology optimization, avoiding re-meshing and the projection of the solution at each iteration. In this work, the domain boundary is represented explicitly by NURBS surfaces in 3D and NURBS curves in 2D, which can be exactly integrated by means of the NEFEM mapping [35], but any other geometrical description, such as e.g. level sets [36, 37] or subdivision surfaces [38], could also be considered.

Local mesh refinement to resolve local features can be implemented in a B-spline context with several approaches, such as T-Splines [39] and hierarchical B-splines (HB-splines) [40, 41, 42]. In this work we consider the latter, mainly due to its straightforward generalization to arbitrary dimensions and its relatively simple implementation.

For the first time to our knowledge, the complete set of boundary conditions is explicitly considered in a numerical solution of the flexoelectric boundary value problem. In the seminal Mindlin’s theory of strain gradient elasticity [43, 44, 45], which is the basis for deriving a stable flexoelectric theory [46, 23, 24, 25], additional non-local boundary conditions are required along non-smooth regions of the domain boundary (i.e. corners in 2D and edges in 3D). This is also the case for the flexoelectric theory. However, in practice, the numerical methods mentioned above neglect these non-local conditions or consider smooth enough domains so that they do not appear [30, 23, 24, 25, 26, 31, 29]. In this work, we show that non-local boundary conditions are mathematically required and we consider them in the formulation and implementation. We demonstrate that neglecting them can deteriorate the solution. In addition, we formulate the boundary conditions corresponding to sensing electrodes, which are common in electromechanical setups.

Within the unfitted framework, a Nitsche’s formulation is derived for the flexoelectric equations to weakly enforce essential boundary conditions, accounting also for the non-local condition from Mindlin’s theory. We show that not only the normal to the boundary, but also the curvatures play a role in the correct enforcement of boundary conditions. A Nitsche’s formulation is also proposed to enforce electrode boundary conditions.

The method converges optimally for high-order approximations of degree $p$ in the $L_{2}$ norm and $H_{s}$ semi-norms, for $s=1,\dots,p$. Namely, it achieves the optimal convergence rates $p+1-s$.

The paper is organized as follows. The variational formulation for flexoelectricity and the associated boundary value problem are presented in Section 2. The numerical approximation based on B-spline approximation and the immersed boundary method are presented in Section 3. Some illustrative numerical examples are given in Section 4. The importance of considering non-local boundary conditions is illustrated in the first example. In the second one, we perform a sensitivity analysis with respect to the Nitsche penalty parameters. In the third one, optimal convergence is tested with a synthetic problem. The remaining examples show 2D and 3D simulations, and compare with available analytical solutions and well-known benchmarks from the literature.

2 Variational formulation and associated boundary value problem

2.1 Notation and preliminary definitions

Let $\Omega$ be a physical domain in $\mathbb{R}^{3}$ . The domain boundary, $\partial\Omega$, can be conformed by several smooth portions as $\partial\Omega=\bigcup_{f}\partial\Omega_{f}$ (Fig. 8). At each point ${\boldsymbol{x}}\in\partial\Omega_{f}$ we define ${\boldsymbol{n}}^{f}$ as the outward unit normal vector. The boundary of the $f$-th portion of $\partial\Omega$ is denoted as $\partial\partial\Omega_{f}$, which is a closed curve. At each point ${\boldsymbol{x}}\in\partial\partial\Omega_{f}$ we define ${\boldsymbol{m}}^{f}$ as the unit co-normal vector pointing outwards of $\partial\Omega_{f}$, which is orthogonal to the normal vector ${\boldsymbol{n}}^{f}$ and to the tangent vector of the curve $\partial\partial\Omega_{f}$, ${\boldsymbol{s}}^{f}$ (see Fig. 8 and 8). The orientation of ${\boldsymbol{s}}^{f}$ is arbitrary and not relevant in the derivations next.

The operators that appear throughout this Section are defined next. The spatial derivative of a function $f$ with respect to the coordinate $x_{i}$ is denoted by $\partial_{i}$ or the subindex $i$ after a comma, that is, $\frac{\partial}{\partial x_{i}}(f)=\partial_{i}f=f_{,i}$. The gradient operator is denoted as $\nabla(~{}~{})$, and the divergence operator as $\nabla\text{\!}\cdot\text{\!}(~{}~{})$. For instance, for a second order tensor ${\mathbf{A}}$, they are defined as $[\nabla({\mathbf{A}})]_{ijk}\coloneqq\partial_{k}A_{ij}=A_{ij,k}$, and $[\nabla\text{\!}\cdot\text{\!}({\mathbf{A}})]_{i}\coloneqq\partial_{j}A_{ij}=A_{ij,j}$, respectively.

The symmetrized gradient $\nabla^{\text{sym}}{(~{}~{})}$ of a vector field ${\boldsymbol{a}}$ is defined as:

[TABLE]

On the domain boundary $\partial\Omega$ the derivative in the normal direction, namely the normal derivate, is denoted by $\partial^{n}$. For instance, the normal derivative of a vector field ${\boldsymbol{a}}$ is defined as $\partial^{n}({\boldsymbol{a}})\coloneqq\nabla({\boldsymbol{a}})\cdot{\boldsymbol{n}}$. The gradient and divergence operators can be decomposed on $\partial\Omega$ into their normal and tangential components as $\nabla(~{}~{})=\partial^{n}(~{}~{}){\boldsymbol{n}}+\nabla^{S}(~{}~{})$ and $\nabla\text{\!}\cdot\text{\!}(~{}~{})=\partial^{n}(~{}~{})\cdot{\boldsymbol{n}}+\nabla^{S}\text{\!}\text{\!}\cdot\text{\!}(~{}~{})$, respectively, where $\nabla^{S}(~{}~{})$ and $\nabla^{S}\text{\!}\text{\!}\cdot\text{\!}(~{}~{})$ denote the surface gradient and surface divergence operators, namely the projection of the gradient and divergence operators onto the tangent space of $\partial\Omega$. For a second order tensor ${\mathbf{A}}$ they are expressed as

[TABLE]

respectively, where ${\mathbf{P}}$ is the projection operator defined on $\partial\Omega$ as $\left[{\mathbf{P}}\right]_{ij}\coloneqq\delta_{ij}-n_{i}n_{j}$, being $\delta_{ij}$ the Kronecker delta.

The formulation involves second-order measures of the geometry, namely curvatures of $\partial\Omega$. The tensor that contains this information is known as the shape operator ${\mathbf{S}}$ (also known as curvature tensor) [47], defined on the surface $\partial\Omega$ as

[TABLE]

The mean curvature $H$ of a surface is an invariant, expressed in terms of ${\mathbf{S}}$ as

[TABLE]

With ${\mathbf{S}}$ and $H$ we define a tensor which arises in the formulation (see A), that we name second-order geometry tensor ${\mathbf{\tilde{N}}}$ and is defined as:

[TABLE]

On the curve $C^{ij}\coloneqq\partial\partial\Omega_{i}\cap\partial\partial\Omega_{j}$, $(i\neq j)$ we define the jump operator \left\llbracket~{}~{}\right\rrbracket acting on a given quantity ${\boldsymbol{a}}$ as the sum of that quantity evaluated at both sides of the curve, namely \left\llbracket{\boldsymbol{a}}\right\rrbracket\coloneqq{\boldsymbol{a}}^{i}+{\boldsymbol{a}}^{j}, where ${\boldsymbol{a}}^{k}$ is the value of ${\boldsymbol{a}}$ from $\partial\Omega_{k}$ (see Fig. 8). Note that a sum is considered in the definition of the jump, so that \left\llbracket{\boldsymbol{m}}\times{\boldsymbol{n}}\right\rrbracket vanishes in the case $C^{ij}$ is defined along a smooth region of $\partial\Omega$, where ${\boldsymbol{m}}^{i}=-{\boldsymbol{m}}^{j}$ and ${\boldsymbol{n}}^{i}={\boldsymbol{n}}^{j}$.

Finally, we denote the first variation of a certain functional $F[f_{1},\dots,f_{n}]$ with respect to the function $f_{i}$ as the functional

[TABLE]

where $\epsilon\in\mathbb{R}$, and the function $\psi_{i}$ is the variation of $f_{i}$ (hence denoted also as ${\delta\!f}_{i}$), defined on the same functional space as $f_{i}$. The variation of $F[f_{1},\dots,f_{n}]$ with respect to all the functions $f_{1},\dots,f_{n}$ is denoted as

[TABLE]

The second variation of $F[f_{1},\dots,f_{n}]$ with respect to the function $f_{i}$ in the direction $\psi_{i}$ is denoted by ${\delta^{2}_{f_{i}}F}[f_{1},\dots,f_{n};\psi_{i}]$, and is defined as the first variation of the functional ${\delta_{f_{i}}F}[f_{1},\dots,f_{n};\psi_{i}]$ with respect to $f_{i}$, i.e.

[TABLE]

where $\psi_{i}$ is the variation of the function $f_{i}$ for both variations of $F$.

2.2 Standard variational formulation

A continuum model for flexoelectric materials can be obtained by coupling a strain-gradient elasticity model with classical electrostatics, through the piezoelectric and flexoelectric effects. The state variables are the displacement field ${{\boldsymbol{u}}}$ and the electric potential $\phi$. The strain tensor and the electric field are given by

[TABLE]

[TABLE]

The bulk energy density $\mathcal{H}^{\Omega}$ in a flexoelectric material can be stated as [46, 23, 24, 25]

[TABLE]

The first two terms correspond to the mechanical energy density of a strain-gradient elastic material, in the Form II of the original paper of Mindlin [43] about strain gradient elasticity. The tensor $\varmathbb{C}_{ijkl}=\varmathbb{C}_{klij}=\varmathbb{C}_{jikl}=\varmathbb{C}_{ijlk}$ is the fourth-order elasticity tensor and $h_{ijklmn}=h_{lmnijk}=h_{jiklmn}=h_{ijkmln}$ is the sixth-order strain-gradient elasticity tensor. The third term is the electrostatic energy density, where $\kappa_{lm}=\kappa_{ml}$ is the second-order dielectricity tensor. The last two terms correspond to the piezoelectric and flexoelectric effects, where $e_{lij}=e_{lji}$ is the third-order piezoelectric tensor and $\mu_{lijk}=\mu_{ljik}$ the fourth-order flexoelectric tensor.

Alternative descriptions could also be considered, writing the bulk energy density in terms of other quantities instead of the electric field ${\boldsymbol{E}}$, such as the polarization ${\boldsymbol{p}}$ or the electric displacement ${\boldsymbol{\hat{D}}}=\varepsilon_{0}{\boldsymbol{E}}+{\boldsymbol{p}}$, where $\varepsilon_{0}$ is the vacuum permittivity constant. As argued in [15], all of them are valid, but the choice of ${\boldsymbol{E}}$ facilitates the derivation of equilibrium equations, which are simpler than those in terms of ${\boldsymbol{p}}$, allowing for simpler numerical methods for solving the associated boundary value problems. Some authors [16, 11, 10, 48] describe the energy density as $\mathcal{H}^{\Omega}[{\boldsymbol{\varepsilon}},\nabla{\boldsymbol{\varepsilon}},{\boldsymbol{p}},\nabla{\boldsymbol{p}}]$ accounting for the polarization gradient theory instead of strain gradient elasticity. We consider the simplified model in Eq. (14) since i) the quadratic term to the electric field gradient is assumed to be negligible for the problems we are interested in, and ii) the last term in Eq. (14) describes both the direct and converse flexoelectric effects as a Lifshitz invariant (see [49, 11, 23] for details), leading to the same equilibrium equations as if both terms were treated separately. We also note that either non-local mechanical or electrical effects (or both) must be introduced in the flexoelectric problem in order to get a meaningful energy density from a physical point of view, but also to get numerically stable formulations. It’s worth noting that incorporating the polarization gradient theory to the present formulation is straightforward by following the same derivations highlighted in the present paper.

The contribution from external loads is presented next. Being ${\boldsymbol{b}}$ the body force per unit volume, and $q$ the free charge per unit volume, their work per unit volume is

[TABLE]

Additional external loads are present on the domain boundary. In a strain gradient elasticity formulation [43], those are the traction ${\boldsymbol{t}}$ and the double traction ${\boldsymbol{r}}$, which are the conjugates of the displacement ${{\boldsymbol{u}}}$ and the normal derivative of the displacement $\partial^{n}({{\boldsymbol{u}}})$ on $\partial\Omega$, respectively. The electrical boundary load is the surface charge density $w$, which is conjugate of the electric potential $\phi$ on $\partial\Omega$. The work of the external loads per unit area is

[TABLE]

Moreover, as dictated by strain gradient elasticity theory [43], an additional force per unit length ${\boldsymbol{j}}$ arises at the edges $C$ of the boundary, i.e. at the union of the edges formed by the intersection of the portions of the boundary, in case $\partial\Omega$ is not smooth (see Fig. 8). That is, at $C=\bigcup_{f}\partial\partial\Omega_{f}$. The force ${\boldsymbol{j}}$ is the conjugate of the displacement ${{\boldsymbol{u}}}$ on $C$. Hence, its work per unit length is

[TABLE]

The total energy $\Pi_{0}[{{\boldsymbol{u}}},\phi]$ of a flexoelectric material is found by collecting all the internal and external energy densities as follows:

[TABLE]

The boundary of the domain $\partial\Omega$ can be split into several disjoint regions, corresponding to the different Dirichlet and Neumann boundaries. The external load is prescribed on the latter, whereas on the former its conjugate is prescribed. For the mechanical loads, we have $\partial\Omega=\partial\Omega_{u}\cup\partial\Omega_{t}$ and $\partial\Omega=\partial\Omega_{v}\cup\partial\Omega_{r}$. The electrical part is also split into $\partial\Omega=\partial\Omega_{\phi}\cup\partial\Omega_{w}$. The corresponding boundary conditions are

[TABLE]

where ${\boldsymbol{\bar{u}}}$, $\overline{{\boldsymbol{v}}}$ and $\bar{\phi}$ are the prescribed displacement, normal derivative of the displacement and electric potential at the Dirichlet boundaries, and ${\boldsymbol{\bar{t}}}$, ${\boldsymbol{\bar{r}}}$ and $\bar{w}$ the prescribed traction, double traction and surface charge density at the Neumann boundaries. The expressions ${\boldsymbol{t}}({{\boldsymbol{u}}},\phi)$, ${\boldsymbol{r}}({{\boldsymbol{u}}},\phi)$ and $w({{\boldsymbol{u}}},\phi)$ will be derived later as a result of the variational principle in Eq. (29).

The edges $C$ of $\partial\Omega$ are also split into $C=C_{u}\cup C_{j}$ corresponding to the Dirichlet and Neumann edge partitions, respectively. Here, $C_{u}$ is assumed to correspond to the curves within the classical Dirichlet boundary, namely $C_{u}=C\cap\overline{\partial\Omega_{u}}$, and $C_{v}=C\setminus C_{u}$. Edge boundary conditions are:

[TABLE]

where ${\boldsymbol{\bar{j}}}$ is the prescribed force per unit length at the Neumann edges, and the expression ${\boldsymbol{j}}({{\boldsymbol{u}}},\phi)$ will be derived later from Eq. (29). Many authors in the literature neglect the edge conditions in Eq. (22) [33, 46, 23, 24]. It is important to note that dismissing them is equivalent to considering homogeneous Neumann edge conditions, which may not be true on $C_{u}$ (Dirichlet edges). In this work, the edge conditions are kept in the formulation to ensure self-consistency and a well-defined boundary value problem.

The energy functional in Eq. (18) can be rewritten as follows, according to Eq. (19b)-(22b):

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The state variables $({{\boldsymbol{u}}},\phi)\in\mathcal{U}_{0}\otimes\mathcal{P}_{0}$, where

[TABLE]

fulfilling Dirichlet boundary conditions in Eq. (19a)-(22a).

The equilibrium states $({{\boldsymbol{u}}}^{*},\phi^{*})$ of the body correspond to the following variational principle:

[TABLE]

The Euler-Lagrange equations associated with the variational principle in Eq. (29) and the expressions ${\boldsymbol{t}}({{\boldsymbol{u}}},\phi)$, ${\boldsymbol{r}}({{\boldsymbol{u}}},\phi)$, $w({{\boldsymbol{u}}},\phi)$ and ${\boldsymbol{j}}({{\boldsymbol{u}}},\phi)$ from the Neumann boundary conditions in Eq. (19b)-(22b) are found by enforcing

[TABLE]

for all admissible variations $\delta{{\boldsymbol{u}}}\in\mathcal{U}_{0}$ and $\delta\phi\in\mathcal{P}_{0}$. The full derivation can be found in [15], and the resulting equations are given next [15, 46, 23, 24, 25]:

[TABLE]

and

[TABLE]

Note that the traction ${\boldsymbol{t}}({{\boldsymbol{u}}},\phi)$ in Eq. (32a) is an alternative expression to the one in [43], whose derivation can be found in A.

In Eq. (31) and (32d), the stress ${\boldsymbol{\hat{\sigma}}}({{\boldsymbol{u}}},\phi)$, the double stress ${\boldsymbol{\tilde{\sigma}}}({{\boldsymbol{u}}},\phi)$ and the electric displacement ${\boldsymbol{\hat{D}}}({{\boldsymbol{u}}},\phi)$ are defined as the conjugates to the strain ${\boldsymbol{\varepsilon}}({{\boldsymbol{u}}})$, the strain gradient $\nabla{\boldsymbol{\varepsilon}}({{\boldsymbol{u}}})$ and the electric field ${\boldsymbol{E}}(\phi)$, respectively, as follows:

[TABLE]

The positivity and negativity conditions on the second variations in Eq. (30b) lead to the following restrictions on the material tensors:

[TABLE]

being $n_{d}$ the number of spatial dimensions.

Remark*.*

This formulation corresponds to the 3D case. It also holds for 2D with the following considerations:

• •

The vectors ${\boldsymbol{n}}$ and ${\boldsymbol{m}}$ defined on the curve $\partial\Omega_{f}$ refer to the outward unit normal and tangent vectors, respectively (see Fig. 12).

• •

The edges $C$ in Eq. (18) correspond to corners of $\partial\Omega$, and therefore have dimension 0. As a consequence, their contribution to the external work $\int_{C}{\boldsymbol{j}}\delta{{\boldsymbol{u}}}\text{\,}\mathrm{d}\mathrm{s}$ is $\sum_{C}{\boldsymbol{j}}\delta{{\boldsymbol{u}}}$ in Eq. (18) and subsequent terms.

2.3 Variational formulation within an unfitted framework: The Nitsche’s method

The admissible space $\mathcal{U}_{0}\otimes\mathcal{P}_{0}$ of the state variables $({{\boldsymbol{u}}},\phi)$ is constrained on the Dirichlet boundaries, and therefore is not suitable for an unfitted formulation. In order to overcome this requirement, an alternative energy functional $\Pi[{{\boldsymbol{u}}},\phi]$ is proposed following Nitsche’s approach [50]:

[TABLE]

where $\Pi^{\text{Dirichlet}}[{{\boldsymbol{u}}},\phi]$ acts on the Dirichlet boundaries instead of $\Pi_{0}^{\text{Dirichlet}}[{{\boldsymbol{u}}},\phi]$, and incorporates Dirichlet boundary conditions in Eq. (19a)-(22a) weakly as follows:

[TABLE]

with the numerical parameters $\beta_{u},\beta_{v},\beta_{C_{u}},\beta_{\phi}\in\mathds{R}^{+}$.

Comparing Eq. (36) against Eq. (25), one can readily see that the expressions ${\boldsymbol{t}}({{\boldsymbol{u}}},\phi)$, ${\boldsymbol{r}}({{\boldsymbol{u}}},\phi)$, $w({{\boldsymbol{u}}},\phi)$ and ${\boldsymbol{j}}({{\boldsymbol{u}}},\phi)$ are now conjugate to the Dirichlet boundary conditions. The penalty terms inserted in each boundary integral are quadratic in the Dirichlet boundary conditions, and its purpose is to weakly enforce Dirichlet boundary conditions and to ensure equilibrium states $({{\boldsymbol{u}}}^{*},\phi^{*})$ being, respectively, actual minima and maxima of the energy functional with respect to ${{\boldsymbol{u}}}$ and $\phi$.

The variational principle associated to $\Pi[{{\boldsymbol{u}}},\phi]$ for the equilibrium states $({{\boldsymbol{u}}}^{*},\phi^{*})$ of the body is stated next:

[TABLE]

where $\mathcal{P}\coloneqq\mathcal{H}^{1}(\Omega)$, and $\mathcal{U}$ is the space of functions belonging to \big{[}\mathcal{H}^{2}(\Omega)\big{]}^{3} with $\mathcal{L}^{2}$-integrable third derivatives on the boundary $\partial\Omega_{u}$, to account for the integrals involving ${\boldsymbol{t}}({{\boldsymbol{u}}},\phi)$ in Eq. (32a). The variational principle in Eq. (37) leads to the same Euler-Lagrange equations in Eq. (31) and definitions of ${\boldsymbol{t}}({{\boldsymbol{u}}},\phi)$, ${\boldsymbol{r}}({{\boldsymbol{u}}},\phi)$, $w({{\boldsymbol{u}}},\phi)$ and ${\boldsymbol{j}}({{\boldsymbol{u}}},\phi)$ in Eq. (32d) as the constrained variational principle in Eq. (29).

The penalty parameters $\beta_{u},\beta_{v},\beta_{C_{u}},\beta_{\phi}$ in Eq. (36) have to be chosen large enough, but too large values would lead to ill-conditioning problems when finding the equilibrium states $({{\boldsymbol{u}}}^{*},\phi^{*})$ numerically. The derivation of stability lower bounds of the penalty parameters can be found in C. However, moderate values of the penalty parameters are enough to ensure convergence and enforce boundary conditions properly [51, 52, 53, 54].

2.4 Weak form of the boundary value problem

The weak form of the unfitted variational formulation is presented next. Vanishing of the first variation of Eq. (35) yields

[TABLE]

where

[TABLE]

being $\delta{{\boldsymbol{u}}}\in\mathcal{U}$ and $\delta\phi\in\mathcal{P}$ admissible variations of ${{\boldsymbol{u}}}$ and $\phi$, respectively.

The functionals ${\delta\Pi^{\Omega}}$ and ${\delta\Pi^{\text{Dirichlet}}}$ are conveniently rearranged as ${\delta\Pi^{\Omega}}[{{\boldsymbol{u}}},\phi;\delta{{\boldsymbol{u}}},\delta\phi]={\delta\Pi^{\Omega}_{\mathscr{B}}}[{{\boldsymbol{u}}},\phi;\delta{{\boldsymbol{u}}},\delta\phi]-{\delta\Pi^{\Omega}_{\mathscr{L}}}[\delta{{\boldsymbol{u}}},\delta\phi]$ and ${\delta\Pi^{\text{Dirichlet}}}[{{\boldsymbol{u}}},\phi;\delta{{\boldsymbol{u}}},\delta\phi]={\delta\Pi^{\text{Dirichlet}}_{\mathscr{B}}}[{{\boldsymbol{u}}},\phi;\delta{{\boldsymbol{u}}},\delta\phi]-{\delta\Pi^{\text{Dirichlet}}_{\mathscr{L}}}[\delta{{\boldsymbol{u}}},\delta\phi]$, with

[TABLE]

The weak form of the unfitted formulation for flexoelectricity reads:

[TABLE]

where

[TABLE]

2.5 Formulation including sensing electrodes

In electromechanics, conducting electrodes are frequently attached to the surface of the devices to enable either actuation or sensing. Actuators induce a deformation due to a prescribed electric potential, whereas sensors infer the deformation state by the measured change in the electric potential. In both cases, as the electrodes are made of conducting material, the electric potential in the electrode is uniform. The electrical Dirichlet boundary condition in Eq. (21) corresponds to actuating electrodes where the uniform electric potential is prescribed. In the case of sensing electrodes, the uniform electric potential is unknown and thus requires a special treatment as described next.

Let us consider a partition of the boundary distinguishing actuating and sensing electrodes, i.e.

[TABLE]

where $\partial\Omega_{\phi}$ and $\partial\Omega_{\Phi}$ correspond, respectively, to actuating and sensing electrodes on the boundary, respectively, and $\partial\Omega_{w}$ to the electrical Neumann boundary. The sensing boundary $\partial\Omega_{\Phi}$ is conformed by $N_{\text{sensing}}$ electrodes, namely $\partial\Omega_{\Phi}=\bigcup_{i=1}^{N_{\text{sensing}}}{\partial\Omega_{\Phi}}^{i}.$

The electric potential $\phi$ on sensing electrodes is constant, but unknown. Thus, at each electrode ${\partial\Omega_{\Phi}}^{i}$ a new state variable ${\Phi}^{i}\in\mathbb{R}$ is introduced, which is a scalar denoting the unknown constant value of the electric potential. In other words,

[TABLE]

Boundary conditions in Eq. (44) are weakly enforced by adding to the energy potential $\Pi[{{\boldsymbol{u}}},\phi]$ in Eq. (35) the work $\Pi^{\text{Sensing}}_{i}[{{\boldsymbol{u}}},\phi,\Phi^{i}]$ of each sensing electrode:

[TABLE]

and the associated variational principle for the equilibrium states $({{\boldsymbol{u}}}^{*},\phi^{*},\Phi^{1*},\dots,\Phi^{{N_{\text{sensing}}}*})$ of the body is

[TABLE]

Equation (45) has a similar form to the Nitsche terms $\Pi^{\text{Dirichlet}}[{{\boldsymbol{u}}},\phi]$ in Eq. (36), but the penalty term quadratic in Eq. (44) is omitted here, because if it is positive in sign then it could be made arbitrarily large with respect to $\phi$, and, conversely, if it is negative in sign then it could be made arbitrarily large with respect to $\Phi^{i}$. Vanishing of the first variation of the energy functional in Eq. (46) yields

[TABLE]

where

[TABLE]

being $\delta\Phi^{i}\in\mathbb{R}$ admissible variations of each $\Phi^{i}$. Finally, the weak form of the unfitted formulation for flexoelectricity accounting for sensing electrodes reads:

[TABLE]

3 Numerical approximation

3.1 B-spline basis

Fourth-order PDEs demand high-order continuity of the functional space for the numerical solution, i.e. the displacement and electric potential fields in the case of electromechanics. Usually, $C^{1}$-continuity of the solution is enough; however, in the unfitted approach presented in this work, many boundary integrals in the Nitsche’s weak forms (see Section 2) involve third-order derivatives, since the test functions do not vanish at the Dirichlet boundaries. It is clear that, in this case, $C^{1}$-continuous solutions are not smooth enough, and therefore we consider approximations of (at least) $C^{2}$-continuity.

A family of functions that provide high-order continuity is that of B-spline functions [55, 56, 57], which are smooth piecewise polynomials. Being $p$ the polynomial degree, they are by construction $C^{p-1}$-continuous throughout the domain. Therefore, cubic ($p=3$) or higher-order B-spline basis are suitable for the numerical approximation of the formulation in Section 2.

Let us consider a uniform B-spline basis; local mesh refinement will be further considered in Section 3.4. Without going into detail, the univariate uniform B-spline basis of degree $p$ consisting of $n$ basis functions is defined on the unidimensional parametric space $\xi\in[0,n+p]$ in terms of the uniform knot vector $\{\xi_{i}\}=[0,1,2,3,\dots,n+p]$. The $i$-th function of this basis is defined recursively as [55]:

[TABLE]

Due to the uniformity of the knot vector, the $i$-th B-spline function can be expressed as a translation of the first ([math]th) one as $B_{i}^{p}(\xi)=B_{0}^{p}(\xi-i)$. Figure 17 shows the function $B_{0}^{p}(\xi)$ of the basis for degrees $p=\{1,\dots,4\}$.

In the three-dimensional space, the $\boldsymbol{i}$-th B-spline function $B_{\boldsymbol{i}}^{p}(\boldsymbol{\xi})$ of a trivariate B-spline basis (where $\boldsymbol{i}$ is the trivariate index $[i_{\xi},i_{\eta},i_{\tau}]$) is defined as the tensor product of three univariate B-spline functions as

[TABLE]

which is defined on the three-dimensional parametric space $\boldsymbol{\xi}\in[0,n_{\xi}+p]\otimes[0,n_{\eta}+p]\otimes[0,n_{\tau}+p]$. Therefore, the parametric space is a cuboid which is defined globally on a Cartesian grid, in contrast with traditional Lagrangian basis present in standard FEM implementations, whose parametric space is defined elementwise.

In the physical space, the problem unknowns ${{\boldsymbol{u}}}$ and $\phi$ are approximated as

[TABLE]

where $\{\boldsymbol{{a^{u}}},\boldsymbol{{a^{\phi}}}\}$ are the degrees of freedom of the numerical solution (known as the control variables in B-spline nomenclature), $\boldsymbol{N}=[\boldsymbol{B}^{p}\circ\boldsymbol{\varphi}^{-1}]$ are the basis functions at the physical space, and $\boldsymbol{\varphi}(\boldsymbol{\xi})$ is the geometric mapping, which is a bijection that maps a given point $\boldsymbol{\xi}$ in the parametric space to a given point $\boldsymbol{x}$ in the physical space. Typically, the map $\boldsymbol{\varphi}(\boldsymbol{\xi})$ is expressed as the interpolation of a discretization of the physical space, namely:

[TABLE]

where $\boldsymbol{S}(\boldsymbol{\xi})$ are the basis functions for the interpolation of the geometry, and $\boldsymbol{\hat{x}}$ are points on the physical space defining the map (known as the control points in B-spline nomenclature).

Different choices of $\boldsymbol{S}(\boldsymbol{\xi})$ and $\boldsymbol{\hat{x}}$ are possible. However, since we want $\boldsymbol{N}(\boldsymbol{x})$ to be $C^{p-1}$-continuous, $\boldsymbol{S}(\boldsymbol{\xi})$ has to be $C^{p-1}$-continuous too, and the most natural choice is $\boldsymbol{S}(\boldsymbol{\xi})\coloneqq\boldsymbol{B}^{p}(\boldsymbol{\xi})$. Therefore, the map $\boldsymbol{x}=\boldsymbol{\varphi}(\boldsymbol{\xi})$ is defined globally. This fact hinders a conforming discretization of the physical space, since it requires an underlying rigid, Cartesian-like mesh in order to be mapped to the parametric space (as done in Isogeometric Analysis [58] and related works). In order to circumvent this requirement on the discretization of the physical space, we consider a different approach where the parametric space of the B-spline basis is not mapped to a conforming discretization of the physical space, but rather to a non-conforming one, naturally providing high-order continuity of the spanned functional space on arbitrary geometries. This concept is known as the immersed boundary approach and is introduced next.

3.2 Immersed boundary method

The main idea of the immersed boundary method, also known as the embedded domain method [59, 60], is to extend the physical domain $\Omega$ to a larger embedding domain $\Omega_{\square}=\Omega\cup\Omega_{\text{fict}}$ which is the one to be discretized, i.e. $\Omega_{\square}=\bigcup_{{\boldsymbol{c}}}\Omega^{{\boldsymbol{c}}}_{\square}$ (see Fig. 21). In this way, the discretization onto the physical space does not depend on the physical domain ($\Omega$) shape. In order to combine this approach with a B-spline basis, the embedding domain $\Omega_{\square}$ is defined as a cuboid, and discretized using a Cartesian-like mesh.

The physical boundary $\partial\Omega$ is allowed to intersect the cells $\Omega_{\square}^{{\boldsymbol{c}}}$ of the embedding mesh arbitrarily. Cells are classified into three different sets $\mathcal{I}$, $\mathcal{C}$ and $\mathcal{O}$, depending on their intersection with the physical domain $\Omega$:

1. i)

$\mathcal{I}\coloneqq\{\Omega_{\square}^{{\boldsymbol{c}}}:\Omega_{\square}^{{\boldsymbol{c}}}\subseteq\Omega\}$, the set of inner cells which remain uncut within the domain (Fig. 21), 2. ii)

$\mathcal{C}\coloneqq\{\Omega_{\square}^{{\boldsymbol{c}}}:\Omega_{\square}^{{\boldsymbol{c}}}\nsubseteq\Omega\text{ and }\Omega_{\square}^{{\boldsymbol{c}}}\cap\Omega\neq\emptyset\}$, the set of cells cut by the boundary (Fig. 21), 3. iii)

$\mathcal{O}\coloneqq\{\Omega_{\square}^{{\boldsymbol{c}}}:\Omega_{\square}^{{\boldsymbol{c}}}\cap\Omega=\emptyset\}$, the set of outer cells, which are neglected.

For the sake of convenience, and without loss of generality, in this work we consider structured cubic Cartesian meshes, which present several practical advantages. On the one hand, all cells are cubes, and a linear mapping can be considered in each cell, namely $\boldsymbol{\varphi}(\boldsymbol{\xi})=\boldsymbol{\breve{x}}_{{\boldsymbol{1}}}+h\boldsymbol{\xi}$, being $h$ the physical cell size and $\boldsymbol{\breve{x}}_{{\boldsymbol{1}}}$ the first corner of the cell. Thus, the Jacobian of the geometric mapping is constant, namely $\boldsymbol{J}(\boldsymbol{\xi})=\nabla_{\xi}\boldsymbol{\varphi}(\boldsymbol{\xi})=h{\mathbf{I}}_{3}$, where ${\mathbf{I}}_{3}$ is the identity matrix of rank 3. On the other hand, in a linear problem all the inner cells of the same size lead to the same elemental matrix, which is computed just once.

At this point, the geometry of $\Omega$ is only used for cell classification. Without going into details, this is usually accomplished by checking whether all vertices of each 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 vertices of each cell (see [61, 62, 63, 64]). For explicit boundary representation (e.g. CAD descriptions), this can be achieved by ray-tracing procedures (see [65, 66]). In this work we restrict ourselves to explicit boundary representation by means of NURBS surfaces in 3D and NURBS curves in 2D.

3.3 Integration on cut cells

Bulk integrals are numerically performed in each cell, i.e. inner ones $\Omega_{\square}^{{\boldsymbol{c}}}\in\mathcal{I}$ and also the physical part of cut ones $\Omega_{\square}^{{\boldsymbol{c}}}\cap\Omega$, for $\Omega_{\square}^{{\boldsymbol{c}}}\in\mathcal{C}$ (see Fig. 21 and 21). Standard cubature rules [67] apply for the former, but not for the latter which can have arbitrary shape. To this end, the physical part $\Omega_{\square}^{{\boldsymbol{c}}}\cap\Omega$ of every cut cell $\Omega_{\square}^{{\boldsymbol{c}}}\in\mathcal{C}$ is divided into several sub-domains (e.g. cuboids or tetrahedra) which are easily integrated. To sub-divide cut cells we rely on the marching cubes algorithm [68], which splits each cell into several conforming tetrahedra, although other conforming [63, 62] or non-conforming [69, 70] subdivision schemes are also possible. See [65] for details of our current implementation. Surface and line integrals are similarly performed on each corresponding sub-domain boundary.

Note that integration sub-domains in contact with the physical domain boundary $\partial\Omega$ might have curved faces or edges in the case $\partial\Omega$ is not flat. Hence, a linear cellwise approximation of the geometry leads to a geometric error of order $2$ which might spoil the optimal convergence of the method. Therefore, cell-wise polynomial approximations of the geometry of degree $p$ are required in general. Alternatively, we exploit the explicit NURBS representation of the geometry by resorting to the NEFEM approach [35, 71, 72, 73, 74] which captures the exact geometry without the need of any polynomial approximation [65].

Remark*.*

The discretization of the weak forms in Eq. (41) and (49), with B-spline basis functions and the mentioned numerical integration, leads to a linear system of equations for the coefficients of the approximation of the unknowns $\{{{\boldsymbol{u}}},\phi\}$, namely $\{\boldsymbol{{a^{u}}},\boldsymbol{{a^{\phi}}}\}$. This linear system typically suffers from ill-conditioning in the presence of cut cells with a small portion in the domain, i.e. when $|\Omega_{\square}^{c}\cap\Omega|\ll|\Omega_{\square}^{c}|$ for a given cell. Ill-conditioning arises basically due to: i) basis functions on the trimmed cell having very small contribution to the integral terms, and ii) basis functions being quasi-linearly dependent on the trimmed cell [75]. Moreover, ill-conditioning is more severe for high-order basis [75]. A detailed investigation on ill-conditioning of immersed boundary methods can be found in [75].

Several strategies have been proposed to alleviate ill-conditioning of trimmed cells, such as the ghost penalty method [76], the artificial stiffness approach [70, 69], the extended B-spline method [77, 78, 79, 54] or special preconditioning techniques specifically designed for immersed boundary methods [75], among others.

For uniform meshes, the extended B-spline approach by Höllig et al. [77, 78, 79, 54] is considered, due to its simple form and good performance. The main idea is to express the critical basis functions on the boundary as linear combinations of inner ones. The constrained basis has less degrees of freedom, but the conditioning and approximation properties are equivalent to those of body-fitted methods [77]. The extension to hierarchical meshes (see Section 3.4) follows the same idea but involves a more sophisticated implementation. In the numerical tests, for the sake of simplicity, hierarchical meshes are stabilized by means of a simple diagonal scaling preconditioning.

3.4 Local mesh refinement: Hierarchical B-spline basis

Hierarchical B-spline refinement was first introduced by Forsey and Bartels [40]. It can be understood as a technique for locally enriching the approximation space by replacing selected coarse B-splines (parents) with finer ones (children). It is based on a remarkable property of uniform B-splines: their natural refinement by subdivision. For a univariate B-spline basis of degree $p$, the subdivision property leads to the following two-scale relation [80]:

[TABLE]

where $\hat{\xi}\coloneqq 2\xi$.

In other words, a B-spline function $B_{i}^{p}(\xi)$ can be expressed as a linear combination of contracted, translated and scaled copies $\hat{B}_{j}^{p}(\hat{\xi})$ of itself [81], as illustrated in Fig. 26 for B-splines of different polynomial degree $p$. The extension to higher dimensions is trivial by means of the tensor product of univariate bases.

Without going into details, a hierarchical B-spline basis is defined from a uniform B-spline basis by replacing some basis functions with their corresponding children (see Fig. 29). This process can be performed recursively, leading to a parent-children hierarchy spanning several levels of refinement. Since each basis function spans several cells, basis refinement implies refinement of multiple cells. The change of focus from element refinement (as in conventional FE) to basis refinement is the key point, which allows maintaining the smoothness of the functional space. Further details can be found in [82, 41, 42, 81, 83] and references therein.

At the implementation level, the elemental matrices of inner cells can be computed just once per level of hierarchy by means of the subdivision projection technique developed in [82]. Therefore, hierarchical B-spline bases maintain the computational benefits of uniform meshes, as explained in Section 3.2, while allowing local mesh refinement.

4 Numerical results

Several numerical simulations are presented next to illustrate the performance of the method. The first example shows the effect of disregarding edge boundary conditions. A synthetic polynomial solution is considered, which can be exactly captured only if edge boundary conditions are enforced. In the second example we perform a sensitivity analysis of the solution and the condition number of the system matrix with respect to the Nitsche penalty parameters on the Dirichlet boundaries. The third example consists on an error convergence analysis considering a non-trivial 2D geometry with curved boundaries and corners, where optimal convergence rates are achieved for different approximation degrees. The fourth and fifth examples deal with two typical setups for flexoelectric characterization, namely a cantilever beam [10] and a truncated pyramid [12]. We compare our simulation results to previous solutions obtained in our group with the maximum entropy meshless method [23] and with approximate analytical solutions from the literature. In the sixth example we present a 3D simulation of a rod with varying semi-circular cross section under torsion, which could be used to measure the shear flexoelectric coefficient [84].

The material tensors in this section are defined next. The mechanics are described by an isotropic elasticity model, with a Young modulus $E$ and a poisson ratio $\nu$, and enriched with an isotropic strain-gradient elasticity model depending on a single length scale parameter $l$. In the 2D case, plane strain conditions are assumed. Electrostatics are described by an isotropic model with dielectricity constant $\kappa_{L}$. Piezoelectricity is described by a tetragonal symmetry model oriented in a certain principal direction ${\boldsymbol{d_{\text{piezo}}}}$. It depends on the longitudinal, transversal and shear piezoelectric coefficients $e_{L}$, $e_{T}$ and $e_{S}$, respectively. Flexoelectricity is described by a cubic symmetry model oriented in the Cartesian axes, with longitudinal, transversal and shear flexoelectric coefficients $\mu_{L}$, $\mu_{T}$ and $\mu_{S}$, respectively. The complete form of every material tensor can be found in B.

Finally, we briefly comment on the choice of penalty parameters in Eq. (36). C shows the derivation of theoretical stability lower bounds. However, for the sake of simplicity, in the following numerical examples, and analogously to other works in the literature [51, 52, 53, 54], we consider the penalty parameters $\beta_{u},\beta_{v},\beta_{C_{u}},\beta_{\phi}\in\mathds{R}^{+}$ in terms of a dimensionless parameter $\zeta\in\mathds{R}^{+}$ as follows:

[TABLE]

where $h$ denotes the physical cell size. We choose $\zeta=100$ for all numerical examples in this work, which is a suitable value as confirmed by the sensitivity analysis in Section 4.2.

4.1 Effect of non-local corner conditions

The effect of neglecting non-local corner conditions is illustrated in the following example. To this end, we simulate the flexoelectric effect on a $[-b,b]\times[-b,b]$ domain (Fig. 34). We consider the following synthetic solution $\{{{\boldsymbol{u}}}^{*},\phi^{*}\}$:

[TABLE]

which can be exactly represented by a cubic ($p=3$) B-spline basis, where $\bar{x}_{i}\coloneqq x_{i}/b$. The prescribed terms on the boundary and source terms on the bulk required for Eq. (56) to be the solution of the flexoelectric problem are computed by inserting Eq. (56) into Eq. (31) and (32d).

Dirichlet boundary conditions are considered at the boundary of the square for mechanics in Eq. (19),(20) and electrostatics in Eq. (21), i.e.

[TABLE]

According to the formulation in Section 2, additional mechanical Dirichlet conditions in Eq. (22) arise at the corners of $\partial\Omega$, namely

[TABLE]

Two cases are considered, depending whether corner conditions in Eq. (58) are enforced or not. In both cases the size $b$ of the domain is varied to assess the effect of corner conditions at different length scales. For the numerical approximation of the solution we consider a cubic ($p=3$) uniform B-spline basis on an unfitted mesh of $32\times 32$ cells (Fig. 34). The material parameters are chosen as follows:

$\quad E=$152\text{,}\mathrm{G}\mathrm{P}\mathrm{a}$;\quad\nu=0.33;\quad l=$1\text{,}\mathrm{nm}$;\quad\kappa_{L}=$141\text{,}\mathrm{nJ}\text{,}{\mathrm{V}}^{-2}\text{,}{\mathrm{m}}^{-1}$;\quad{\boldsymbol{d_{\text{piezo}}}}=x_{2};\quad e_{L}=$8.8\text{,}\mathrm{J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-2}$;$

$e_{T}=$-4.4\text{,}\mathrm{J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-2}$;\quad e_{S}=$4.4\text{,}\mathrm{J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-2}$;\quad\mu_{L}=$150\text{,}\mathrm{\SIUnitSymbolMicro J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-1}$;\quad\mu_{T}=$110\text{,}\mathrm{\SIUnitSymbolMicro J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-1}$;\quad\mu_{S}=$110\text{,}\mathrm{\SIUnitSymbolMicro J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-1}$.$

Figure 34 shows the $L_{2}$ error of the numerical solution at different length scales, computed as the $L_{2}$ norm of the error vector ${\boldsymbol{\upepsilon}}\coloneqq(u_{1}-{u_{1}}^{*},u_{2}-{u_{2}}^{*},\phi-\phi^{*})$ on $\Omega$. Results reveal that the synthetic solution is exactly captured (up to round-off errors) only in the case corner conditions are considered. Otherwise, a significant error is found, which grows as the domain size is reduced, i.e. at scales where strain-gradient elasticity and flexoelectricity couplings are relevant. The source of this error is not numerical, since the geometry is exactly represented and integrated, and the approximation space captures the analytical solution exactly.

The error introduced by neglecting non-local corner conditions is illustrated in Fig. 34 for the scale $b/b_{0}=100$, where $b_{0}\coloneqq l=$1\text{,}\mathrm{nm}$$ is a characteristic length of the problem used for normalization. The absolute value of the electric potential component of the error ${\boldsymbol{\upepsilon}}$, i.e. $|\phi-\phi^{*}|$, is depicted in $\log$-scale within $\Omega$. The error is concentrated around the corners of the domain, where the two boundary value problems defer. Similar behavior is observed for the mechanical components of the error ${\boldsymbol{\upepsilon}}$.

We conclude that non-local corner conditions are a fundamental part of the mathematical prescription of the physical problem of linear flexoelectricity in the presence of non-smooth domains. Ignoring them leads to solving a different boundary value problem, and the resulting discrepancies can be very significant below length scales where strain-gradient effects play a relevant role. Therefore, in the following numerical examples, non-local corner conditions are always properly considered.

4.2 Sensitivity analysis on the penalty parameter

We perform a sensitivity analysis on the penalty parameter $\zeta$ in Eq. (55). We consider the setup in the previous example, with two different meshes leading to minimum volume fractions $\chi_{\text{min}}\approx 93.7\%$ and $\chi_{\text{min}}\approx 0.9\%$ respectively, where $\chi_{\text{min}}\coloneqq\min_{c}\left(\frac{|\Omega_{\square}^{c}\cap\Omega|}{|\Omega_{\square}^{c}|}\right)$ denotes the volume fraction of the smallest cut cell in the mesh.

According to Eq. (30b), the second variations of the energy functional with respect to the mechanical and electrical unknowns are required to be positive and negative, respectively, consistent with the $\min$-$\max$ nature of the problem and to ensure stability of the formulation. Numerically, this is equivalent to checking that the minimum eigenvalue of ${\mathbf{K}}_{{{\boldsymbol{u}}}{{\boldsymbol{u}}}}$ and the maximum eigenvalue of ${\mathbf{K}}_{\phi\phi}$ are positive and negative respectively, where ${\mathbf{K}}_{{{\boldsymbol{u}}}{{\boldsymbol{u}}}}$ and ${\mathbf{K}}_{\phi\phi}$ are the submatrices of the system related to mechanical and electrical equations and unknowns, respectively.

Fig. 37 shows the minimum eigenvalue of ${\mathbf{K}}_{{{\boldsymbol{u}}}{{\boldsymbol{u}}}}$ and the maximum eigenvalue of ${\mathbf{K}}_{\phi\phi}$ as a function of $\zeta$. The $\min$-$\max$ condition stated above is fulfilled for $\zeta>10$. It is important noting that this threshold does not depend on the minimum volume fraction $\chi_{\text{min}}$ of the mesh. This result is in agreement with [53], where the threshold of the penalty parameter $\zeta$ is analytically proven to be mesh-independent after using cell-aggregation stabilization schemes on unfitted meshes, as the extended B-splines method we consider in this work (see remark in Section Remark).

The sensitivity of both the numerical error, ${\boldsymbol{\upepsilon}}$, and of the condition number of the system matrix, ${\kappa_{1}}$, on the penalty parameter $\zeta$ is shown in Fig. 40. For both meshes, machine-precision accuracy is achieved for moderate values of $\zeta$, i.e. $10<\zeta<1000$. However, for larger values of $\zeta$ the errors increase due to the increase in the condition number, which grows proportionally to $\zeta$. In view of the reported results, we consider $\zeta=100$ for the following numerical examples.

4.3 Accuracy and convergence properties of the method

In this example the convergence of the method is assessed for different high-order B-spline approximations, namely at $p=\{3,4\}$. We consider the following synthetic solution $\{{{\boldsymbol{u}}}^{*},\phi^{*}\}$:

[TABLE]

depicted in Fig. 43 on the domain $\Omega$, which is conformed by a circle of radius $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ with a square-shaped hole of size $0.4\text{\times}0.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ rotated 30 degrees with respect to the Cartesian coordinates; both geometries are centered at ${\boldsymbol{x}}=(0,0)$.

Dirichlet conditions are enforced at the boundaries and the corners of the domain. The boundary is exactly mapped by means of the NEFEM mapping, which makes the geometrical error vanish. The numerical integration is rich enough so that the integration error is negligible and does not pollute the convergence plots.

The material parameters are chosen as follows:

$\quad E=$100\text{,}\mathrm{G}\mathrm{P}\mathrm{a}$;\quad\nu=0.37;\quad l=$2\text{,}\mathrm{nm}$;\quad\kappa_{L}=$11\text{,}\mathrm{nJ}\text{,}{\mathrm{V}}^{-2}\text{,}{\mathrm{m}}^{-1}$;\quad{\boldsymbol{d_{\text{piezo}}}}=x_{2};\quad e_{L}=$8.8\text{,}\mathrm{J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-2}$;$

$e_{T}=$-4.4\text{,}\mathrm{J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-2}$;\quad e_{S}=$1.1\text{,}\mathrm{J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-2}$;\quad\mu_{L}=$0.5\text{,}\mathrm{\SIUnitSymbolMicro J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-1}$;\quad\mu_{T}=$1\text{,}\mathrm{\SIUnitSymbolMicro J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-1}$;\quad\mu_{S}=$0.5\text{,}\mathrm{\SIUnitSymbolMicro J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-1}$.$

Error convergence results for cubic $(p=3)$ and quartic $(p=4)$ B-spline bases are shown in Fig. 43. The error ${\boldsymbol{\upepsilon}}$ is measured in the $L_{2}$ norm and $H_{1}$, $H_{2}$ and $H_{3}$ semi-norms at six recursively-refined uniform meshes of cell size $h/h_{0}=\{$2\text{\times}{10}^{-3}$,\dots,$2\text{\times}{10}^{-8}$\}$, where $h_{0}\coloneqq$1\text{,}\mathrm{\SIUnitSymbolMicro m}$$ is a normalization factor. As expected, optimal convergence rates are (asymptotically) achieved with both B-spline bases. The asymptotic behavior of the error convergence rate is expected, and is due to the extended B-spline stabilization performed on trimmed basis functions near the boundary $\partial\Omega$ (see remark in Section Remark). For relatively coarse meshes, a small additional error is introduced since the approximation space is coarsened at the boundary. However, for fine enough meshes, this effect is negligible and error convergence rates tend to optimality [77, 78].

4.4 Transversal transduction: cantilever beam under bending

Bending a cantilever beam is a natural way to mobilize transversal strain-gradients and, consequently, to trigger the transverse flexoelectric effect. For this reason, electroactive beam bending has been extensively used in experiments for the characterization of the transversal flexoelectric coefficient [85, 86, 87, 88, 89]. This setup has also been modeled with approximate analytical [48] and numerical [23] models.

Figure 47 depicts the geometrical model of a cantilever beam. The aspect ratio of the beam is fixed to $L/a=20$, where $L$ represents the length of the beam and $a$ its width. We consider different lengths $L$ in order to capture the size-dependent nature of the flexoelectric coupling. Mechanically, the beam is clamped on its left-end and undergoes a point force $F$ on its top-right-corner. Electrically, it is grounded on its right-end, and the other edges are considered charge-free. The corresponding boundary conditions are:

[TABLE]

Strain-gradient elasticity is neglected to isolate the effect of piezoelectricity and flexoelectricity couplings. Three different cases are considered for the electromechanics: i) piezoelectricity only, ii) flexoelectricity only and iii) combined piezoelectricity and flexoelectricity.

The beam accommodates the mechanical load by bending, which produces a linear distribution of the axial strain $\varepsilon_{11}$ along the transversal ($x_{2}-$) direction which is well known in classical elasticity. Namely, both i) strain and ii) strain gradients are generated, which are the triggers for direct piezoelectric and flexoelectric effects, respectively. Therefore, a non-zero electric field is generated on the sample as a consequence of the mechanical loading. For a piezoelectric beam, we expect the electromechanical response to be the same regardless of the size of the sample. However, for a flexoelectric or flexo-piezoelectric beam it should grow inversely to the scale due to the size-dependent nature of the strain gradient field and thus of flexoelectric coupling [23].

In order to quantify the energy conversion, we define the electromechanical coupling factor $k_{\text{eff}}$ as

[TABLE]

which is a positive, dimensionless scalar that indicates the relationship between dielectric and mechanic energies required to accommodate the external mechanical load.

An analytical estimation of $k_{\text{eff}}$ for the flexo-piezoelectric beam can be found in the literature [48] as

[TABLE]

where $\varepsilon_{0}\approx$8.854\text{\times}{10}^{-12}\text{,}\mathrm{C}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-1}$$ is the vacuum permittivity constant, and it has been assumed that the material parameters $\nu=l=\mu_{L}=\mu_{S}=e_{L}=e_{S}=0$ (see B).

One can also define the normalized effective piezoelectric constant $e^{\prime}$ [48] as

[TABLE]

which indicates the ratio between the current $k_{\text{eff}}$ and the $k_{\text{eff}}$ that would be obtained if the beam was purely piezoelectric. By combining Eq. (62) and (63), analytical estimations of $e^{\prime}$ are obtained for flexoelectric and flexo-piezoelectric beams as

[TABLE]

where $a^{\prime}\coloneqq-ae_{T}\mu_{T}^{-1}$ is the normalized beam thickness.

For the numerical approximation of the solution let us consider a cubic ($p=3$) hierarchical B-spline basis on a mesh with two levels of refinement around the left-end and top-right-corner of the beam, as depicted in Fig. 47 and 47. The material parameters are the following:

$\nu=l=\mu_{L}=\mu_{S}=e_{L}=e_{S}=0;\quad E=$100\text{,}\mathrm{G}\mathrm{P}\mathrm{a}$;\quad\kappa_{L}=$11\text{,}\mathrm{nJ}\text{,}{\mathrm{V}}^{-2}\text{,}{\mathrm{m}}^{-1}$;\quad e_{T}=$-4.4\text{,}\mathrm{J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-2}$;\quad\mu_{T}=$1\text{,}\mathrm{\SIUnitSymbolMicro J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-1}$.$

The normalized effective piezoelectric constant $e^{\prime}(a^{\prime})$ is depicted in Fig. 52, which shows very good agreement between our numerical results and the analytical estimations in Eq. (64). The flexoelectric effect is negligible at large scales; as a consequence, for large $a^{\prime}$ the electromechanical response tends to vanish for the flexoelectric beam, whereas the flexo-piezoelectric one tends to behave as purely piezoelectric. On the contrary, at smaller scales the flexoelectric effect is much more relevant and leads to an enhanced electromechanical transduction in both cases.

Figures 52-52 depict the spatial distribution of the mechanically-induced electric field ${\boldsymbol{E}}$ at scale $a^{\prime}=1.76$ for the cases of piezoelectric, flexoelectric and flexo-piezoelectric beams. In the piezoelectric case, the distribution is skew-symmetric (divergent) with respect to the neutral axis of the beam, in accordance with the axial strain $\varepsilon_{11}$ distribution. Highest values appear close the left-end and away from the neutral axis. The other two cases have a similar $k_{\text{eff}}$ but very different distributions of the electric field ${\boldsymbol{E}}$. The electric field on the flexoelectric beam points downwards, and remains almost constant along the transversal direction while increases close to the clamped end. The flexo-piezoelectric beam, however, presents an inhomogeneous distribution which can be thought as a combination of the two previous ones. Depending on the scale, the piezoelectric effect dominates the flexoelectric one or viceversa. Here, $a^{\prime}=1.76$ is intentionally chosen since it leads to comparable electromechanical effects.

4.5 Longitudinal transduction: truncated pyramid under compression

Another frequent experimental setup is the truncated pyramid compression, widely used by experimentalists to characterize the longitudinal flexoelectric coefficient [12, 90, 91, 92]. Although analytical expressions are not available for this more complex setup, numerical solutions have been developed by our group [23].

Figure 57 depicts the geometrical model of a flexoelectric truncated pyramid of height $a$ and bases $a$ (top) and $3a$ (bottom). The angle between bases and lateral boundaries is $\pi/4$. Mechanically, the bottom basis is fixed and a compressive force $F$ is uniformly distributed on the top one. Electrically, it is grounded on the top basis, and a sensing electrode is placed at the bottom. The corresponding boundary conditions are:

[TABLE]

where $V$ is a priori unknown but constant.

Since the top and bottom bases have different sizes, they undergo different compressive tractions and a longitudinal strain gradient in the $x_{2}$-direction arises as a result, which triggers the longitudinal flexoelectric effect. Another source of strain gradient is the bottom layer being fixed, which results in an inhomogeneous distribution of the traction along the bottom boundary. Therefore, a non-zero electric field is generated on the sample as a consequence of the mechanical loading. A direct measure of the electromechanical transduction is the value $V$ of the sensing electrode at the bottom of the truncated pyramid. More interestingly, one can compute the effective electric field $E_{\text{Eff}}$ measured as the voltage difference between electrodes over the height of the pyramid, i.e.

[TABLE]

Numerical results are obtained with a cubic ($p=3$) B-spline basis on a uniform unfitted mesh of cell size $h\approx 0.01195a$, and shown in Fig. 57 for $F=$4.5\text{,}\mathrm{N}\text{,}{\mathrm{mm}}^{-1}$$, with the following material parameters:

$l=\mu_{S}=e_{L}=e_{T}=e_{S}=0;\quad E=$100\text{,}\mathrm{GPa}$;\quad\nu=0.37;\quad\kappa_{L}=$11\text{,}\mathrm{nJ}\text{,}{\mathrm{V}}^{-2}\text{,}{\mathrm{m}}^{-1}$;\quad\mu_{L}=\mu_{T}=$1\text{,}\mathrm{\SIUnitSymbolMicro J}\text{,}{\mathrm{V}}^{-1}\text{,}{\mathrm{m}}^{-1}$.$

Very good agreement with previous works in the literature [23] is reported. The size-dependent nature of the flexoelectric coupling is evidenced in Fig. 57, which shows the effective electric field $E_{\text{Eff}}$ as a function of the size of the truncated pyramid. Electromechanical transduction of the device takes place mainly at the micro- and nanoscale, whereas it is not relevant at larger scales.

In order to illustrate the complexity of the physics, the vertical strain $\varepsilon_{22}$ and electric potential $\phi$ distributions at scale $a=7.5~{}$\mathrm{\SIUnitSymbolMicro m}$$ are depicted in Fig. 57 and 57, respectively. A highly inhomogeneous distribution of the strain takes place, specially near the corners of the grounded electrode on top of the device, which causes large strain-gradients triggering the flexoelectric effect. As a consequence, the electric potential distribution is also inhomogeneous within the domain. For this reason, simplified 1D models are not reliable to simulate the flexoelectric truncated pyramid, and numerical simulations are required [23].

4.6 Shear transduction: conical semicircular rod under torsion

Unlike the longitudinal and transversal flexoelectric coefficients, the shear coefficient has been scarcely characterized experimentally. One reason is that, in many setups such as in the cylindrical rod torsion, shear strain gradients are effectively mobilized but the overall net polarization vanishes, and therefore no flexoelectric measurement can be effectively done. An alternative setup proposed by Mocci et. al. [84] consists on a conical rod with semicircular cross section under torsion, where a net angular polarization arises thanks to the longitudinal variation of the cross section.

Figure 60 shows the geometrical model of the conical semicircular rod, with a length of $100~{}$\mathrm{\SIUnitSymbolMicro m}. The radii of the semicircular bases are $26.3~{}$\mathrm{\SIUnitSymbolMicro m} and $7.5~{}$\mathrm{\SIUnitSymbolMicro m}, and their centers are located at ${\boldsymbol{x}}_{O}=(0,0,0)~{}$\mathrm{\SIUnitSymbolMicro m} and ${\boldsymbol{x}}_{o}=(100,0,0)~{}$\mathrm{\SIUnitSymbolMicro m}$$.

The larger semicircular basis is clamped and grounded, and torsion is enforced at the opposite basis by prescribing the displacement field. The corresponding boundary conditions are:

[TABLE]

where $\alpha$ is the tangent of the prescribed torsional angle.

The mechanical response of the rod is composed by several effects, including non-constant twisting (in-plane rotation) and warping (out-of-plane displacement). Without going into the details, one can think of the rod undergoing $\varepsilon_{12}$ and $\varepsilon_{13}$ shear strains varying along the $x_{1}$ direction, hence triggering the shear flexoelectric effect along the $x_{2}-x_{3}$ planes.

Numerical simulations are performed with a cubic $(p=3)$ trivariate B-spline basis on an unfitted uniform mesh of cell size $1.778~{}$\mathrm{\SIUnitSymbolMicro m}$$ (see Fig. 60). The prescribed torsion is set to $\alpha=0.1$, which corresponds to a counterclockwise torsion of about $5.7\text{\,}\mathrm{\SIUnitSymbolDegree}$. The material constants are set to match those of barium strontium titanate, a strongly flexoelectric ceramic, in its paraelectric phase:

[TABLE]

In order to isolate the shear component of the flexoelectric effect, two simulations are performed. In the first one, only shear flexoelectricity is taken into account, namely $\mu_{L}=\mu_{T}=0$. In the second one, the complete flexoelectricity tensor ${\boldsymbol{\mu}}$ is considered.

Numerical results are shown in Fig. 63. The electric potential takes positive values for $x_{2}>0$ and negative values otherwise, being more prominent near the free end. An effective electric field arises in the polar direction contained in the $x_{2}-x_{3}$ plane [84], which can be readily seen by plotting the electric potential in a cross section of the rod (see Fig. 66). This distribution allows us to measure the electric potential difference between both sides of the rod, and therefore can be used to measure the shear flexoelectric coefficient [84].

Results do not vary much by considering or disregarding the longitudinal and transversal coefficients of the flexoelectric tensor (see Subfig. (a) and (b) in Fig. 63 and 66). In order to quantify it, the voltage difference at the corners of the cross sections in Fig. 66, namely at ${\boldsymbol{x}}_{+}=(90.5,9.265,0)~{}$\mathrm{\SIUnitSymbolMicro m} and ${\boldsymbol{x}}_{-}=(90.5,-9.265,0)~{}$\mathrm{\SIUnitSymbolMicro m}, is evaluated for the two cases, yielding:

[TABLE]

which shows that considering the longitudinal and transversal coefficients of the flexoelectric tensor affects the voltage difference only by 2.87%. Therefore, it is apparent that the flexoelectric behavior of this setup is mainly controlled by the shear flexoelectric coefficient $\mu_{S}$.

5 Concluding remarks

We have developed a computational approach with unfitted meshes to simulate the electromechanical response of small scale dielectrics at infinitesimal strains, including the piezoelectric and flexoelectric electromechanical couplings. The high-order nature of the equations is addressed by smooth B-spline basis functions on a background Cartesian mesh, which can be hierarchically refined to resolve local features. The unfitted nature of the method allows B-spline-based simulations on arbitrary domain shapes, which can be explicitly represented by means of NURBS surfaces (3D) or curves (2D), and exactly integrated by means of the NEFEM mapping. Therefore, the method is suitable for the rational design and optimization of nanoscale electromechanical devices with no geometrical limitations. The formulation has been detailed for infinitesimal deformations. The ideas described here are currently being extended to a finite deformation setting suitable for the study of flexoelectricity in soft materials as well [34].

Our work highlights two features of the flexoelectric formulation that have been scarcely commented in the literature, namely a) the correct way of considering non-smooth boundaries (i.e. corners in 2D) and their corresponding non-local boundary conditions, and b) the role that the curvature of the boundary plays in physical quantities such as the mechanical tractions.

The Nitsche’s method has been particularized to the flexoelectric problem to weakly enforce essential boundary conditions and sensing electrode conditions. Optimal high-order error convergence rates are reported on non-trivial geometries featuring both curved boundaries and corners.

We have simulated several electromechanical setups that are traditionally used to quantify the longitudinal and transversal components of the cubic flexoelectric tensor and are standard benchmarks in the literature of linear flexoelectricity, such as the bending of a cantilever beam and the compression of a truncated pyramid. In all cases, our results match the ones in the literature. Additionally, we perform a 3D simulation of the torsion of a conical semicircular rod, to illustrate the ability of the proposed method to deal accurately with complex geometries. The results are in excellent agreement with [84].

The proposed computational framework can assist the design and optimization of a new generation of nanoscale electromechanical devices such as actuators, sensors and energy harvesters, allowing any complexity on the domain shape.

Acknowledgments

This work was supported by the Generalitat de Catalunya (“ICREA Academia” award for excellence in research to I.A., and Grant No. 2017-SGR-1278), and the European Research Council (StG-679451 to I.A.).

Appendix A Particular expression for the traction

Equation (32a) in Section 2 presents a particular expression for the traction ${\boldsymbol{t}}({{\boldsymbol{u}}},\phi)$, which is not the usual one in strain-gradient elasticity or flexoelectricity formulations in the literature [43, 44, 45, 23, 93, 46, 30, 11, 16], which is:

[TABLE]

We obtain Eq. (32a) by expanding and rearranging terms in Eq. (68) in the following way:

[TABLE]

where ${\mathbf{S}}$ is the shape operator of a surface and $H=\frac{1}{2}Tr(\,{\mathbf{S}}\,)$ its mean curvature [47], as defined in Section 2.

Equation (A) reveals that the traction ${\boldsymbol{t}}({{\boldsymbol{u}}},\phi)$ has two contributions: the former involves stress measures, i.e. ${\boldsymbol{\hat{\sigma}}}({{\boldsymbol{u}}},\phi)$ and $\nabla\text{\!}\cdot\text{\!}{\boldsymbol{\tilde{\sigma}}}({{\boldsymbol{u}}},\phi)$, dotted with the normal vector ${\boldsymbol{n}}$ (a first-order measure of the geometry), whereas the latter involves the double stress measure ${\boldsymbol{\tilde{\sigma}}}({{\boldsymbol{u}}},\phi)$ dotted with second-order geometry measures, namely the second-order geometry tensor defined as ${\mathbf{\tilde{N}}}:={\mathbf{S}}-2H{\boldsymbol{n}}\otimes{\boldsymbol{n}}$. Thus, it is clear that high-order physics are intrinsically linked to high-order geometrical measures of the domain, such as the curvature of its boundary.

Appendix B Material tensors

In the following, $n_{d}$ refers to the number of dimensions of the physical space $\Omega$ (either 2 or 3). Material tensors are described component-wise, and only the non-zero components are specified.

Isotropic elasticity is represented by the fourth-order tensor $\varmathbb{C}$, which depends on the Young modulus $E$ and the Poisson ratio $\nu$ as

[TABLE]

where the parameters $C_{L}$, $C_{S}$ and $C_{T}$ are

[TABLE]

in the 3D and plane strain 2D cases, and

[TABLE]

for the plane stress 2D case.

We consider an isotropic simplified strain gradient elasticity model [94], which is a particular case of the general model of strain gradient elasticity in [44]. The strain gradient elasticity tensor is represented by the sixth-order tensor ${\mathbf{h}}$, which depends on the Young modulus $E$, the Poisson ratio $\nu$ and a single internal length scale $l$ in the following form:

[TABLE]

with parameters $C_{L}$, $C_{S}$ and $C_{T}$ defined in Eq. (71) and (72).

Isotropic dielectricity is represented by the second-order tensor ${\boldsymbol{\kappa}}$, which depends on the parameter $\kappa_{L}$ as

[TABLE]

Piezoelectricity is represented by the third-order tensor ${\boldsymbol{e}}$. Tetragonal symmetry is considered, which has a principal direction and involves longitudinal, transversal and shear couplings represented by the parameters $e_{L}$, $e_{T}$ and $e_{S}$, respectively. For a material with principal direction ${\boldsymbol{x}}_{1}$, the piezoelectric tensor ${\boldsymbol{e}}^{<{\boldsymbol{x}}_{1}>}$ reads

[TABLE]

The piezoelectric tensor ${\boldsymbol{e}}$ oriented in an arbitrary direction ${\boldsymbol{d}}$ is obtained by rotating ${\boldsymbol{e}}^{<{\boldsymbol{x}}_{1}>}$ as

[TABLE]

where ${\mathbf{R}}$ is a rotation matrix that rotates ${\boldsymbol{x}}_{1}$ onto ${\boldsymbol{d}}$.

Flexoelectricity is represented by the fourth-order tensor ${\boldsymbol{\mu}}$. Cubic symmetry is considered, which leads to a flexoelectric tensor involving longitudinal, transversal and shear couplings represented by the parameters $\mu_{L}$, $\mu_{T}$ and $\mu_{S}$, respectively. We refer to [22] for an extensive analysis of other possible symmetries for the flexoelectric tensor. The components of the flexoelectric tensor ${\boldsymbol{\mu}}^{<{{\boldsymbol{x}}}>}$ of a material oriented in the Cartesian axes are the following:

[TABLE]

The flexoelectric tensor ${\boldsymbol{\mu}}$ oriented in an arbitrary orthonormal basis is obtained by rotating ${\boldsymbol{\mu}}^{<{\boldsymbol{x}}>}$ as

[TABLE]

where ${\mathbf{R}}$ is the rotation matrix that rotates the Cartesian basis to the desired orthonormal basis.

Assuming the material models presented in this Appendix, the restrictions on material tensors in Eq. (34) simplify to restrictions on material coefficients as follows:

[TABLE]

Appendix C Computation of penalty parameters

Nitsche’s method involves numerical penalty-like parameters to enforce Dirichlet boundary conditions. The stability of the formulation is guaranteed for large values of the parameters. Too large values would lead to ill-conditioning of the system matrix but, differently to penalty methods, usually moderate values provide good results. Lower bounds of the penalty parameters can be assessed globally (constant for the whole mesh) [95] or locally (cell-wise) [96], being the latter more appealing due to i) lower condition number of the resulting algebraic system and ii) lower computational cost [75].

Typically, for elliptic PDE, lower bounds are found by studying the coercivity of the formulation (see for instance [95]). In the case of flexoelectricity, it is a saddle point problem corresponding to the coupling between mechanical (positive definite) and electrical (negative definite) problems. Therefore, coercivity is met by checking the positivity and negativity of the second variations of the energy functional with respect to the mechanical and electrical unknowns, respectively, as stated in Eq. (30b):

[TABLE]

From Eq. (80), one can readily see that the second variations of the energy functional depend solely on $\delta{{\boldsymbol{u}}}$ or $\delta\phi$. Therefore, the coupling physics (i.e. piezoelectricity and flexoelectricity) do not play any role in determining the coercivity of the formulation. As a consequence, the coercivity analysis can be performed for the uncoupled mechanical and electrical problems independently, namely

[TABLE]

being $\Pi^{\text{Mechanical}}\coloneqq\Pi|_{\delta\phi=0}$ and $\Pi^{\text{Electrical}}\coloneqq\Pi|_{\delta{{\boldsymbol{u}}}={\boldsymbol{0}}}$.

The mechanical part of the problem corresponds to a strain-gradient elasticity formulation, for which we derive next the lower bounds arising from the coercivity analysis. The electrical part of the problem corresponds to a (negative) Poisson equation, which has been largely studied in the literature [95] and whose corresponding lower bounds are well known and not presented here.

C.1 Deriving conditions for coercivity of the strain-gradient elasticity formulation

The bilinear form of the strain-gradient elasticity formulation corresponds to the mechanical part of the flexoelectric bilinear form in Eq. (42a):

[TABLE]

It is easy to verify that the second variation ${\delta^{2}_{{{\boldsymbol{u}}}}\Pi^{\text{Mechanical}}}[\delta{{\boldsymbol{u}}}]$ corresponds to $\mathcal{B}[\delta{{\boldsymbol{u}}},\delta{{\boldsymbol{u}}}]$. Since $\mathcal{B}[{{\boldsymbol{u}}},\delta{{\boldsymbol{u}}}]$ can be expressed as $\mathcal{B}[{{\boldsymbol{u}}},\delta{{\boldsymbol{u}}}]=\sum_{c}\mathcal{B}^{\Omega_{\square}^{c}\cap\Omega}[{{\boldsymbol{u}}},\delta{{\boldsymbol{u}}}]$, coercivity of $\mathcal{B}[{{\boldsymbol{u}}},\delta{{\boldsymbol{u}}}]$ is cell-wise met by proving coercivity of $\mathcal{B}^{\Omega_{\square}^{c}\cap\Omega}[{{\boldsymbol{u}}},\delta{{\boldsymbol{u}}}]$, $\forall\Omega_{\square}^{c}\in\mathcal{I}\cup\mathcal{C}$. Therefore, a sufficient cell-wise coercivity condition is:

[TABLE]

Since Eq. (83) is trivially fulfilled on inner cells, the coercivity analysis is carried out only on cut cells:

[TABLE]

where the penalty parameters $\beta^{c}_{u},\beta^{c}_{v},\beta^{c}_{C_{u}}\in\mathds{R}^{+}$ correspond to each cell $\Omega_{\square}^{c}$.

Applying the Cauchy-Schwarz inequality and the Young’s inequality leads to:

[TABLE]

which holds $\forall\epsilon_{u},\epsilon_{v},\epsilon_{C_{u}}\in\mathds{R}^{+}$. Let us consider now the mesh-dependent constants $K_{u},K_{v},K_{C_{u}}\in\mathds{R}^{+}$ such that

[TABLE]

for all admissible $\delta{{\boldsymbol{u}}}$, which can be computed as detailed in C.3. Then, Eq. (85) leads to:

[TABLE]

Thus, the conditions

[TABLE]

are sufficient conditions for the coercivity of $\mathcal{B}^{\Omega_{\square}^{c}\cap\Omega}[\delta{{\boldsymbol{u}}},\delta{{\boldsymbol{u}}}]$, for any $\epsilon_{u},\epsilon_{v},\epsilon_{C_{u}}\in\mathds{R}^{+}$.

C.2 Conditions on the penalty parameters

In order to get explicit bounds of the penalty parameters, let us define $\alpha_{1}\coloneqq\epsilon_{u}/\epsilon_{v}\in\mathds{R}^{+}$ and $\alpha_{2}\coloneqq\epsilon_{u}/\epsilon_{C_{u}}\in\mathds{R}^{+}$. The first coercivity condition in Eq. (88) is rewritten as

[TABLE]

By means of Eq. (89), we can rewrite the remaining conditions in Eq. (88) as

[TABLE]

Since Eq. (90) ensure coercivity for any $\epsilon_{u}\in\mathds{R}^{+}$, we can express the bounds of the penalty parameters as a function of $\alpha_{1}$ and $\alpha_{2}$ only.

The bilinear form $\mathcal{B}^{\Omega_{\square}^{c}\cap\Omega}$ is coercive if $\beta_{u}^{c},\beta_{v}^{c},\beta_{C_{u}}^{c}$ satisfy

[TABLE]

for any $\alpha_{1},\alpha_{2}\in\mathds{R}^{+}$.

At this point, it is worth mentioning that Eq. (91) correspond to a particular cell $\Omega_{\square}^{c}\in\mathcal{C}$ such that i) $\Omega_{\square}^{c}\cap\partial\Omega_{u}\neq\emptyset$ (classical Dirichlet boundary), ii) $\Omega_{\square}^{c}\cap\partial\Omega_{v}\neq\emptyset$ (non-local Dirichlet boundary) and iii) $\Omega_{\square}^{c}\cap{C_{u}}\neq\emptyset$ (non-local Dirichlet edges), i.e. for cut cells with the complete set of Dirichlet boundary conditions. In the case of a cell with a Neumann boundary condition, it is easy to verify that Eq. (91) still holds if we just consider the bounds and constants corresponding to the Dirichlet boundary conditions. For instance, in a cell where double tractions ${\boldsymbol{r}}$ are prescribed on the boundary, $\Omega_{\square}^{c}\cap\partial\Omega_{v}=\emptyset$ and we obtain the following bounds:

[TABLE]

Any choice of $\alpha_{1},\alpha_{2}\in\mathds{R}^{+}$ leads to a coercive formulation; however, the condition number of the resulting linear system might be affected. Further investigation is required to assess suitable values for $\{\alpha_{1},\alpha_{2}\}$ for providing condition numbers of the system as low as possible.

C.3 Computation of the mesh-dependent constants $K_{u}$, $K_{v}$ and $K_{C_{u}}$

Mesh-dependent constants $K_{u}$, $K_{v}$ and $K_{C_{u}}$ satisfying Eq. (86) can be taken as the largest eigenvalues of the following generalized eigenvalue problems that arise from the spatial discretization of Eq. (86):

[TABLE]

The matrices $\mathbf{B_{u}},\mathbf{B_{v}},\mathbf{B_{C_{u}}}$ and $\mathbf{V}$ are defined at each cell $\Omega_{\square}^{c}$ as the discretization of the weak forms $\mathcal{B}_{u}^{\Omega_{\square}^{c}\cap\Omega}$, $\mathcal{B}_{v}^{\Omega_{\square}^{c}\cap\Omega}$, $\mathcal{B}_{C_{u}}^{\Omega_{\square}^{c}\cap\Omega}$ and $\mathcal{B}_{\Omega}^{\Omega_{\square}^{c}\cap\Omega}$ in Eq. (86), respectively, with

[TABLE]

The numerical solution of the generalized eigenvalue problems requires careful consideration [75], since a) the matrix $\mathbf{V}$ is always singular (it is based solely on derivative quantities of the test functions) and b) the matrix $\mathbf{V}$ can be bad-conditioned if the corresponding volume fraction $\chi_{c}$ of the cell is very small. We refer to [75] for further details.