Three-Dimensional Isogeometric Analysis of Flexoelectricity with MATLAB Implementation

Flexoelectricity is a general electromechanical phenomenon where the electric polarization exhibits a linear dependency to the gradient of mechanical strain and vice versa. The truncated pyramid compression test is among the most common setups to estimate the flexoelectric effect. We present a three-dimensional isogeometric formulation of flexoelectricity with its MATLAB implementation for a truncated pyramid setup. Besides educational purposes, this paper presents a precise computational model to illustrate how the localization of strain gradients around pyramidal boundary shapes contributes in generation of electrical energy. The MATLAB code is supposed to help learners in the Isogeometric Analysis and Finite Elements Methods community to learn how to solve a fully coupled problem, which requires higher order approximations, numerically. The complete MATLAB code which is available as source code distributed under a BSD-style license, is provided in the part of Supplementary Materials of the paper.


Introduction
In dielectric crystals with non-centrosymmetric crystal structure such as quartz and ZnO, electrical polarization is generated upon the application of uniform mechanical strain. This property of certain materials, which is known as piezoelectricity, is caused by relative displacements between the centers of oppositely charged ions. Flexoelectricity which has recently attracted significant attention; however, differs from the piezoelectricity in two aspects: 1) it is a more general effect which is available in all dielectrics including those with centrosymmetric crystal structures and 2) the induced electrical polarization is related to the gradient of mechanical strain and is thus a size dependent effect [Sharma, Maranganti and Sharma (2007); Yudin and Tagantsev (2013)]. The compression of a truncated pyramid is a common setup to measure flexoelectric properties of dielectrics [Zhu, Fu, Li et al. (2006); Abdollahi, Millan, Peco et al. (2015); Huang, Shu, Kwon et al. (2014); Kwon (2017)]. The pyramidal boundary shapes introduce strain gradients and thus induce electrical voltage. By recording these two quantities, one can quantify flexoelectricity [Cross (2006)]. Sometimes the experimental and theoretical measurements of flexoelectric coefficients differ, noticeably [Sharma, Landis and Sharma (2010)]. One potential reason is lack of a precise mathematical model. Available analytical models [Sharma, Landis and Sharma (2010); Maranganti, Sharma and Sharma (2006); Eliseev, Morozovska, Glinchuk et al. (2009); Mao and Purohit (2014)] are not capable to address the complexity of the deformation fields; particularly around the edges of the pyramid. There are also computational models [Qi, Huang, Fu et al. (2018); Shen and Hu (2010); Mao, Ai, Xiang et al. (2016); Yvonnet and Liu (2017); Mao, Purohit and Aravas (2016); Nguyen, Zhuang and Rabczuk (2018)]. Abdollahi et al. [Abdollahi, Peco, Millan et al. (2014); Abdollahi, Peco, Millian et al. (2015)] presented a meshfree model in 2D. Thai et al. [Thai, Rabczuk and Zhuang (2018)] presented a large deformation isogeometric approach for flexoelectricity. He et al. [He, Javvaji and Zhuang (2019)] implemented element-free Galerkin method to characterize flexoelectricity in a composite material. Sidhardh et al. [Sidhardh and Ray (2018)], presented a numerical model to obtain the effective properties of flexoelectric fiber-reinforced nanocomposite. The authors of this paper already presented an IGA model for flexoelectricity in two-dimensional space [Ghasemi, Park andRabczuk (2017, 2018)]. It is noteworthy to mention that, only a few of available computational models are in 3D [Abdollahi, Millan, Peco et al. (2015); Deng, Deng and Shen (2018) 2019)]. This is our motivation to extend our previous work into three-dimensional space for precisely obtain the mechanical strains near the edges of the pyramid. Proposed by Hughes and his co-workers [Hughes, Cottrell and Bazilevs (2005)], the basic idea behind IGA was to unify Computer Aided Design (CAD) and Computer Aided Engineering (CAE). IGA shows some advantages in comparison to the classical FEM. Among them one can point to the exact representation of the geometry, ease in adaptivity and mesh refinement, accuracy in imposing the essential boundary conditions and the higher and controllable continuity at the inter-element boundary. Here, IGA benefits us with compact support high order B-spline basis functions to discretize the fourth order partial differential equations of flexoelectricity; which demand at least 1 continuous basis functions in Galerkin method [Abdollahi, Millan, Arroyo et al. (2014)]. We present a MATLAB implementation for a 3D isogeometric formulation of flexoelectricity. We provide necessary tutorials for each section of the code, making it straightforward to follow. Hopefully, it will contribute to the popularity of the topic. The remainder of this paper is organized thus: Section 2 summarizes the discretized governing equations of flexoelectricity, Section 3 explains MATLAB implementation, Section 4 discusses about results and Section 5 offers concluding remarks.

A summary of the governing equations and discretization 2.1 Governing equations
A summary of the governing equations of the flexoelectricity is presented in this section. More details are available in [Abdollahi, Peco, Millan et al. (2014); Ghasemi, Park and Rabczuk (2017)] and references therein. Accounting for the flexoelectricity, the enthalpy density, ℋ, can be written as where is the fourth-order elasticity tensor, is the mechanical strain, is the electric field, is the fourth-order flexoelectric tensor, ℎ is the sixth-order straingradient elasticity tensor and is the second-order dielectric tensor.
The different stresses / electric displacements including the usual ( � / � ), higher-order ( � / � ) and physical ( / ) ones are then defined through the following relations which are the governing equations of the flexoelectricity. By imposing boundary conditions and integration over the domain, Ω, the total electrical enthalpy is Using Hamilton's principle, we finally have which is the weak form of the governing equations of the flexoelectricity. In Eq. (8), is the mechanical displacements, is the electric potential, ̅ is the prescribed mechanical tractions and is the surface charge density. Γ and Γ are boundaries of Ω corresponding to mechanical tractions and electric displacements, respectively.

IGA discretization
There are two different spaces in IGA namely the physical space ( Fig. 1(c)) and parameter space ( Fig. 1(b)). Knot vectors discretize the parameter space. A knot vector in one dimension is a non-decreasing set of coordinates in the parameter space. We call knot vectors corresponding to , and directions as = { 1 , 2 , … . , + +1 } , = { 1 , 2 , … . , + +1 } and = { 1 , 2 , … . , + +1 } respectively. is called the knot index and , , ∈ ℝ are the ℎ knot in each direction. , and are the number of basis functions; while , and are their polynomial orders. The parent element ( Fig. 1(a)) is used for numerical integration. If knots are equally-spaced in the parametric space, they are said to be uniform. A knot vector is said to be open if its first and last knots appear + 1 times. ,0 ( ) = � and for = 1,2,3, … assuming that 0 0 = 0. The derivatives of basis functions, , ′ ( ) are given by ,  Basis functions can be enriched by increasing their order. They are − continuous (i.e., there are − continuous derivatives) across element boundaries, where is polynomial order and is multiplicity of each knot value. This higher continuity is our main motivation to use B-spline basis functions to discretize the fourth order PDE of flexoelectricity. Fig. 3 shows that, the support of a B-spline function of order is + 1 knot spans; it means , is non-zero over { , + +1 }.There is a notion of -refinement in IGA, which is an order elevation to followed by a knot insertion, to obtain −1 continuity of the basis functions at inserted knot. Assuming , ( ), , ( ) and , ( ) to be univariate B-spline basis functions of order , and corresponding to knot vector , and , B-spline basis functions in three dimensional space, , , , , ( , , ) are presented as  Fig. 1(c)) in IGA are used to discretize the geometry and define the degrees of freedom. They do not necessarily lie on the surface itself, but define its envelope. As shown in Fig. 1, each element in the physical space is the image of a corresponding element in the parameter space (and vice versa) via mapping and −1 , respectively. Solution fields (e.g., ℎ ) can be similarly approximated via tensor product of nodal values (either coordinates or solution fields) with their corresponding basis functions where � , , � is the corresponding control points of the element.
where the superscript denotes nodal parameters at the control points; � contain the spatial derivatives of the B-spline basis functions. and are Hessian matrices containing the second derivatives of the basis functions: The discrete system of Eq. (8) is expressed as In Eqs. (14a)-(14f), the subscript, , in Ω , Γ and Γ denotes the ℎ finite element where Ω = ⋃ Ω . and are mechanical and electrical load vectors. We consider isotropic elasticity and permittivity tensors and adopt cubic symmetry for flexoelectric tensor. Thus, , , and read as follows:

MATLAB implementation
Due to difference in surface areas, the truncated pyramid experiences different tractions on the top and bottom surfaces. This results in strain gradient and consequently, electrical polarization. We consider a pyramid subject to a uniform pressure on its top surface with dimensions and material properties listed in Tab. 1. The boundary conditions are demonstrated in Fig. 4. We use cubic B-spline elements, unless otherwise specified.

Element structure and connectivity
The function CUBIC_Mesh generates cubic elements with eight vertices in the parametric space. The inputs of the function are number of elements in , and directions as well as the number of nodes per element which is eight. The function returns the nodes coordinates and connectivity as configured in Fig. 5. The function cppolygonPrint_3D_truncated returns coordinates of the control points. For a typical quadratic element, the corresponding control points are identified according to the pattern shown in Fig. 6. The function pageomapping_3D_truncated (Listing 1) maps all elements nodes from the parameter space onto the physical space.

Assembling stiffness and force matrices
The function formStiffness_3D_Truncated_flexo (Listing 2) returns the global stiffness matrix of the pyramid. The variable, depicts the total number of degrees of freedom, which is four times the total number of control points. In fact, on each control point there are three components of the displacement field and an electric potential, [ ] . The solution field is also obtained as )×abs(wt×detJacob×J2); end % loop over Gauss points end %loop over elements

Boundary conditions
To apply the equipotential boundary condition on the top surface of the pyramid, penalty stiffness method is used as described in Listing 5. The electric potential degrees of freedom corresponding to the top nodes are tied to each other one by one in both and directions. The script applyBC_3D (see Listing 6) is used to impose mechanical boundary conditions as well as the zero-electric potential boundary condition on the bottom edge.

Post processing
The function VTKPostProcess3d, is adopted to visualize the numerical results using ParaView which is an open-source, multi-platform data analysis and visualization platform. The electric potential and strain through thickness direction are plotted for each case. The results show acceptable conformity with results presented in Abdollahi et al. [Abdollahi, Millan, Peco et al. (2015)]. It is observable that an increase in causes larger strain gradients and consequently more electric polarization. The same trend is also found for = 3 as shown in Fig. 8.

Extensions
A suite of extensions can be thought for the presented code. For instance, change in supporting and boundary conditions can be mentioned. Flexible supports instead of rigid type can cause bending of the pyramid, besides its compression. The inverse flexoelectric effect can be investigated by applying electric potential on the electrodes and measure the deformation of the pyramid. A flexoelectric multi-pyramid composite can be designed using multi patches.

Concluding remarks
We present a three-dimensional isogeometric formulations of flexoelectricity with its MATLAB implementation for a truncated pyramid setup. We take advantages of the Bspline basis functions to discretize the fourth order partial differential equations of flexoelectricity; which demand at least 1 continuous basis functions in a Galerkin method. Our numerical simulation clarifies that the strain gradients are highly localized around the edges of the pyramid and their magnitudes increase by an increase in the pyramid area ratio, . We provide MATLAB code as Supplementary Materials of the paper with explanations to facilitate the popularity of the topic.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.
The 3 × 3 matrix is called which is the Jacobian matrix of the transformation. For a quadratic B-spline basis functions corresponding to an element with 9 control points, we have 9 basis functions 1 , 2 , … , 9 . Thus, To calculate the second derivatives we write: Putting Eqs. (A5a)-(A5f) in the matrix form we will have: