Second-Order Collocation-Based Mixed FEM for Flexoelectric Solids

: Flexoelectricity is an electromechanical coupling between the electric ﬁeld and the mechanical strain gradient, as well as between the mechanical strains and the electric ﬁeld gradient, observed in all dielectric materials, including those with centrosymmetry. Flexoelectricity demands C 1 -continuity for straightforward numerical implementation as the governing equations in the gradient theory are fourth-order partial differential equations. In this work, an alternative collocation-based mixed ﬁnite element method for direct ﬂexoelectricity is used, for which a newly developed quadratic element with a high capability of capturing gradients is introduced. In the collocation method, mechanical strains and electric ﬁeld through independently assumed polynomials are collocated with the mechanical strains and electric ﬁeld derived from the mechanical displacements and electric potential at collocation points inside a ﬁnite element. The mechanical strain gradient and electric ﬁeld are obtained by taking the directional derivative of the independent mechanical strain and electric ﬁeld gradients. However, an earlier proposed linear element is unable to capture all mechanical strain gradient components and, thus, simulate ﬂexoelectricity correctly. This problem is solved in the present work by using quadratic shape functions for the mechanical displacements and electric potential with fewer degrees of freedom than the traditional mixed ﬁnite element method. A Fortran user-element code is developed by the authors: ﬁrst, for the linear and, after that, for the quadratic element. After verifying the linear element with numerical results from the literature, both linear and quadratic elements’ behaviors are tested for different problems. It is shown that the proposed second-order collocation-based mixed FEM can capture the ﬂexoelectric behavior better compared to the existing linear formulations.


Introduction
With the theoretical identification of the flexoelectric effect by Mashkevich and Tolpygo in 1957 [1] and the quantification of flexoelectric coefficients by Kogan in 1964 [2], a new chapter in the research on microelectromechanical systems (MEMS) began. Like other electromechanical effects, such as piezoelectricity, this opened the possibility of constructing new types of sensors and actuators or energy-harvesting devices. Unlike piezoelectricity, which only occurs in noncentrosymmetric crystals, flexoelectricity breaks the inversion symmetry of the crystal structure by inducing mechanical strain gradients. Therefore, flexoelectricity is observed in all dielectrics and is a universal electromechanical effect [3]. This property makes flexoelectricity particularly relevant for materials with high dielectric parameters, such as ferroelectrics [4,5].
Flexoelectricity can be subdivided into direct and converse flexoelectric effects. The first describes the generation of an electric field due to mechanical strain gradients, whereas the converse effect defines the coupling between the electric field gradients and mechanical strains. Typically, the flexoelectric coefficients have small values. High mechanical strain gradients must be induced for flexoelectricity to become a dominant effect. Because flexoelectricity is inversely related to the length scale, this is more easily accomplished in micro-and nanostructures. Due to the ongoing miniaturization in microelectronics, size-dependent effects play an increasingly important role, making flexoelectricity an essential topic in recent years. Various experiments including cantilever beam and truncated pyramid setups [6][7][8][9] have been conducted to quantify flexoelectricity. In order to utilize the flexoelectric effect in electronic devices, analytical formulations have been developed, e.g., for flexoelectricity in one-dimensional nanosized cantilever beams. Based on the theory developed in [10], the formulations for a Bernoulli-Euler beam are derived, in which piezoelectricity and direct flexoelectricity are considered [11]. However, with more complex structures involved, the problems become very difficult to solve analytically. Therefore, different numerical methods were developed to solve this issue. The most popular one is the conventional finite element method (FEM); however, it is not directly suitable for simulating flexoelectricity because of the requirement of C 0 -continuous elements. Consequently, mesh-free formulation [12], isogeometric analysis [13,14], moving least square [15], the hierarchical B-spline method [16], and mixed FEM have been developed to simulate flexoelectricity. For the latter method, numerical results for the "plate with a hole" problem [17,18] and the "infinite length tube" problem [19] were obtained. Furthermore, recently in [20], numerically robust mixed finite elements were proposed for modeling size-dependent flexoelectric behavior in piezoelectric solids to highlight mutual interactions of piezoelectricity and flexoelectricity. It is also possible to use C 1 -continuous elements with higher-order shape functions and additional degrees of freedom (DOFs) as given in [18].
If choosing mixed FEM for flexoelectric behavior simulation, the constraints between the mechanical displacement field and its gradient are enforced by Lagrange multipliers [17]. These Lagrange multipliers are set as extra DOFs. For a 2D quadrilateral element, 27 DOFs are used for mechanical displacements and electric potential on all nine nodes, 16 DOFs are needed for displacement gradients at corner nodes and four Lagrange multipliers are involved at the center node [19]. With this high number of DOFs involved, the stiffness matrix, which has to be inverted to solve the numerical problem, becomes weighty, leading to a low computational efficiency for this element. Furthermore, the standard mixed/hybrid FEM has some drawbacks, like the requirement to satisfy the Ladyzhenskaya-Babuška-Brezzi condition [21] which is, however, not possible a priori.
In the present work, the collocation-based mixed FEM is used. In this method, the mechanical strains derived from the nodal displacements are collocated with independently estimated mechanical strains assumed to be a polynomial. The collocation is done on specific points inside the finite element. These points should be located at the Gaussian points, as pointed out in [22]. The collocation-based mixed FEM was extended for the collocation of an electric and magnetic contribution reducing the sensitivity to mesh distortion and aspect ratio compared to the displacement-based elements [23]. Ref. [24] extended the collocation method for higher-order and 3D elements. A collocation mixed FEM to simulate flexoelectric behavior was developed in [25]. It uses C 0 -continuity and is capable of accounting for the size effect of the mechanical strain gradients. The mechanical strain gradients are thereby obtained by taking the directional derivatives of the independently assumed mechanical strains. This linear collocation-based element only has 12 DOFs corresponding to two mechanical displacements and one electric potential DOF at each node of a four-noded element. Ref. [26] was recently extended for the simulation of semiconductors.
In the proposed manuscript, the collocation-based mixed FEM [25] is further extended, introducing biquadratic elements. For this, independent mechanical strains and the electric field are assumed as a quadratic polynomial. This increases the number of DOFs to 27. However, there are fewer DOFs involved than in traditional mixed FEM [17,19,27]. In the current work, a Fortran user element (UEL) for ABAQUS is developed from scratch to simulate flexoelectricity by using linear and quadratic elements. First, the linear element is verified by comparing the results given in [25] with the results of the linear element developed by the authors, for which few specific mechanical strain gradient components are neglected. Later, the performance of the newly developed quadratic element was compared to the linear element and studied in detail. The numerical examples presented in our study are both the cantilever beam and the truncated pyramid problems.
The paper is organized as follows: Section 2 introduces the governing equations for flexoelectric material behavior. Sections 3.1 and 3.2 present the collocation mixed FEM formulations for the linear element and introduces the quadratic element by outlining the new developments in the collocation method. Section 4.1 presents the result in which the cantilever beam is used as a numerical example. There, in Section 4.1.1 the validation of the results of the developed linear UEL with the numerical results obtained in [25] is shown. The newly developed quadratic element is validated in Section 4.1.2 by neglecting specific mechanical strain gradient components and comparing the obtained results with the results of the linear element for the cantilever beam problem. Additionally, a convergence of linear and quadratic elements are studied in Section 4.1.3. By comparing the linear and the quadratic element in Section 4.1.4, considerable differences in the results are observed and analyzed. Furthermore, in Section 4.1.5, the material behavior for both element types is investigated for a wide range of flexoelectric coefficients. Also in Section 4.1.6, possibilities for further improvements in the physical correctness of the obtained result for the quadratic element are discussed. In Section 4.2, the numerical results for the truncated pyramid are presented, in which the difference in the electric response is outlined. The final conclusions and discussion are presented in Section 5.

Flexoelectricity: Constitutive and Governing Equations
The internal energy density U, which accounts for flexoelectricity [10], can be summarized as where c ijkl and a ij are the components of the elastic and dielectric permittivity tensors, respectively, e kij are the components of the piezoelectric tensor, f ijkl , and b klij are the component of the direct and converse flexoelectricity tensors, whereas g ijklmn and h ijkl are the components of the strain gradient elasticity (SGE) tensor and the higher-order electric field gradient tensor, respectively. The mechanical strain tensor ε, mechanical strain gradient tensor η and the electric field E are defined as where u i are the components of the mechanical displacement and φ is the electric potential. From the internal energy density Equation (1) the constitutive equations for stress σ ij , higher-order stress τ kji , electric displacement D i and higher-order electric displacement Q ij are obtained as partial derivatives with respect to ε ij , η kji , E i and E i,j , respectively: In order to account for the size effect of the mechanical strain gradients and electric field gradients, internal length scale parameters l and q are introduced. These parameters are multiplied with the elastic stiffness coefficients c ijkl and the dielectric coefficients a ij to get the higher-order elastic parameters g jklmni and the higher-order electric parameters h ijkl of Equations (5b) and (5d), respectively [25,28]: Here, δ ij denotes the Kronecker delta and c jkmn and a ik are given for orthotropic material as For a 4 mm tetragonal crystal structure, the piezoelectric coefficients e kij are given by [25,29] where the indices 1 and 3 describe two dimensions of a 2D Cartesian coordinate system in which 3 is the poling direction. So the stress due to piezoelectricity can be written as The direct and converse flexoelectric tensors have low symmetry. For a cubic crystal, they can be expressed by three independent coefficients [30]. For converse flexoelectricity, the coefficients b klij are represented as b 1 , b 2 and b 3 [30]: The stress due to converse flexoelectricity can be written as For direct flexoelectricity, the coefficients f ijkl can be reduced further to two independent components f 1 and f 2 when isotropic material is considered [19,30]: By applying the Voigt notation, the constitutive Equations (5a)-(5d) can be presented as [25] with material tensors from Equations (6)-(10), (12) and (14)   The governing equations for a piezoelectric solid with direct and converse flexoelectric effects in the absence of the body forces and free charges can be written as [18,25] σ ij,j − τ ijk,jk = 0, (17a) The Dirichlet boundary conditions are prescribed as on the boundary Γ with normal derivative of mechanical displacement w i and normal derivative of electric potential p. Γ u , Γ w , Γ φ and Γ p describe mechanical displacement boundary, boundary of normal derivative of mechanical displacement, electrical potential bound-ary, and boundary of normal derivative of electric potential, respectively. The bar over the variables indicates that it is a prescribed value. The Neumann boundary conditions are where R i denotes the higher-order traction and Z is the higher-order surface charge. Traction boundary, higher-order traction boundary, electrical surface charge boundary and higher-order surface charge boundary are denoted as Γ t , Γ R , Γ S and Γ Z , respectively: The traction vector t i and the electrical surface charge S are defined as [25] where ρ i = n k π j τ ijk and α = n i π j Q ij . π i describes the Cartesian component of the unit tangent vector on boundary Γ, and δ(x) is the Dirac delta function. ||ρ i (x k )|| and ||α(x k )|| mean the jump at a corner x k on the oriented boundary contour Γ.

Collocation Mixed Finite Element Method (CMFEM)
For the numerical simulation of the direct and converse flexoelectricity, one needs to calculate mechanical strains, mechanical strain gradients, electric field, and electric field gradients. As mechanical strain gradients and electric field gradients are the second derivatives of mechanical displacement and electric potential, respectively, C 1 -continuity is required in traditional FEM. To avoid C 1 -continuous elements due to their complexity, mixed FEM is traditionally used for simulating flexoelectricity [17,19,27]. To reduce the number of DOFs and increase computational efficiency, CMFEM, used in this study, was proposed for flexoelectricity and extensionally developed by [25]. In CMFEM, independently assumed mechanical strains and independently assumed electric field are introduced and collocated with the mechanical strains and electric field computed from the mechanical displacement u i and electric potential φ. u i and φ are interpolated from the nodal values u a i and φ a inside the finite element, respectively: Here, n is the total number of nodes and r = r s T is the position vector in the local curvilinear coordinate system. In order to calculate the gradient in the global coordinate system, the following coordinate transformation equation is needed ( Figure 1) Therefore, the Jacobian matrix J in Equation (24) is calculated with the derivatives of the shape functions N a and the coordinates of each node x a i . We have The global derivative of shape function N a with respect to x 1 and x 3 are introduced as b a 1 and b a 3 , respectively, by multiplying Equation (24) with shape functions N a : The mechanical strains and electric field are derived from the mechanical displacements and the electric potential, respectively, Here, B a,T ε11 (r), B a,T ε33 (r) and B a,T ε13 (r) are the strain-displacement conectivity matrix components

CMFEM for Linear Element
The idea of the collocation method is that a polynomial is assumed for each component of the independent mechanical strain and electric field. These polynomials consist of a P-vector and a coefficient vector β or α Here, the superscript "In" describes an independent quantity. The shape of the P-vector and the coefficient vector depends on whether the collocated, independent component is assumed to be constant, linear, or quadratic. For the linear element, all components of mechanical strain and electric field are collocated linearly. Consequently, the independent P-vector is defined as follows: The corresponding coefficient vectors have to be assumed accordingly, so that the number of components is the same as for the assumed P vector The left top indices distinguish different coefficients. Because the coefficients i α and i β for the assumed polynomials are not known initially, these have to be computed by using the collocation method. For this, the quantities which are derived from nodal values are set to be equal to the independently assumed quantities at collocation points x c . The collocation points must be at specific positions inside the finite element, which have to be the Gaussian points to pass the patch test [22,23]. Consequently, ε ij and E i are collocated at the 2×2 Gaussian points 5, 6, 7, and 8 in Figure 2. We have By assembling P-vectors for all collocation points, A matrix is obtained. A contains the coordinates of all considered collocation points, and hence can be precalculated for each Gaussian point as where The subscript of r and s denotes the collocation points. The B a -matrices in Equations (35) and (36) can be written as follows: From Equations (35) and (36) the coefficient vectors are obtained and inserted into Equations (30) and (31) as where u and φ contain the mechanical displacement and electric potential for all nodes, respectively, and the newly introduced B a u -and B a φ -matrices are being calculated as By taking the directional derivative of P, mechanical strain gradients and electric field gradients are obtained [25] as where (44) Applying the Voigt notation for the mechanical strain gradients and the electric field gradients, the following equations are obtained: The variational condition describes the equivalence of the variational internal work and the variational external work. In gradient theory for FEM, it can be derived as [10] In order to implement described formulations, the stiffness matrix for the flexoelectric material has to be derived. Therefore, the expressions for mechanical strains, mechanical strain gradients, electric field, and electric field gradients obtained by the collocation method and the constitutive Equations (5a) and (5d) are inserted into the variational condition Equation (47): Here, e denotes an element and V e the volume of an element.t a andR a are given as where the superscript is the node number. Because the variation δu and δφ can be arbitrary, Equation (48) splits into two equations which can be written for each element as with Algorithm 1 presents schematically how the collocation method is implemented in UEL.
The inverse Jacobian matrix for all integration points (ip) 2: 3: compute B u,N according to Equation (27 compute P ,1 and P ,3 Equation (44) The directional derivatives 14: of the P-vector for current ip 15: 16:

CMFEM for Quadratic Element
In this subsection, CMFEM is derived for a quadratic element. In this quadratic element, quadratic shape functions [31] are used for interpolating the mechanical displacement and electric potential. This creates the possibility to assume the independent mechanical strains ε ij , as well as the independent electric field E i to be quadratic. So in this study, they are collocated over nine collocation points which are the 3 × 3 Gaussian points 10,11,12,13,14,15,16,17,18 in Figure 3 [24]. For such an element, P-vector could be defined as P(r) = 1 r s rs r 2 s 2 r 2 s rs 2 r 2 s 2 .
The corresponding coefficient vectors take the following form where the left top index indicates that these are not the same coefficients as in Equations (33a) and (33e) : The A-matrix changes to A = A quad with A quad as The B a -matrices given in Equations (35) and (36) are required to be The advantage of this quadratic element is that it can capture all mechanical strain gradient components in contrast to the linear element, which is only capable of calculating strain gradient components η 133 , η 311 , η 131 and η 313 . In the linear element, the components η 111 and η 333 cannot be captured. This means that the generated electric potential due to flexoelectricity and the element stiffness is essentially more accurate with this quadratic element.

Numerical Examples
This section investigates the performance of the linear and the proposed quadratic element. To begin with, the developed linear element is verified in the literature. A first numerical example serves the cantilever beam problem. Thereafter, the correctness of the quadratic element is verified by comparing it with the linear element results, setting particular strain gradient components to zero in the quadratic element. In a convergence study, the performance of both element types is evaluated, and the differences in the results are investigated and analyzed. Thereafter, a parameter study in which the influence of flexoelectric coefficients f 1 and f 2 on the deflection and the generated electric field measures both linear and quadratic elements. Finally, the performance of the newly developed quadratic element is illustrated on the truncated pyramid under the compression problem and verified with the literature.

Verification of Linear Element
In this subsection, the results obtained from the developed UEL, based on Section 3.1, are compared with the numerical results of the linear collocation mixed FEM element given in [25]. For this purpose, a nanoscaled cantilever beam problem with a transverse force on its end as in [25] is used as an example. The problem is assumed to be 2D plane strain.
As depicted in Figure 4, the beam has a length of L = 500 nm, a width of B = 10 nm a height of H = 20 nm. The whole left side of the beam is horizontally fixed and one node is fixed in the vertical direction. At the free end of the beam over the right side, a force F is applied. Furthermore, the free end of the beam is electrically grounded. The boundary conditions are summarized in Table 1.   [25].

Side/Corner
Mechanical BC Electrical BC Prescribed Force Piezoelectric material PZT-5H is used. As in [25], the material is chosen to be isotropic with parameters presented in Table 2. Table 2. Material properties for the cantilever beam [25].
A Poisson's ratio of ν = 0 and not the stated value of ν = 0.2 is used for the cantilever beam simulation in [25]. This can be shown by analytical calculation of the deflection for the case of the linear elastic beam and comparing the result with the analytical solution given in the deflection plot. So ν is set accordingly in the current work and only direct flexoelectricity is considered. In [25], 1600 elements for the cantilever beam are used. The same number of elements is used in the current simulations to compare both results better. Figure 5 shows the distribution of the deflection components u 3 and u 1 , mechanical strain ε 11 , and mechanical strain gradient η 311 inside the cantilever beam. One can see that not only u 3 but also u 1 is built up and so there exists mechanical strain ε 11 . The gradient of mechanical strain ε 11 along x 3 is much higher than any other mechanical strain gradient component, leading to the generation of the electric field. Deflection u 3 of the cantilever beam along its length is presented in Figure 6 for different flexoelectric coefficients It can be seen that by accounting for the mechanical strain gradient effect and increasing the flexoelectric coefficient f 1 , the beam becomes stiffer and the deflection decreases. The electrical response is shown in Figure 7. Electric potential is mainly built up on the fixed side of the beam where mechanical strain gradient η 311 attains the highest value. In Figure 8, the electric field along the length of the cantilever beam is illustrated. One can see that the electric field has its maximum value at the fixed end and reduces linearly to the free end. From Figures 6 and 8 it is apparent that the results obtained with the developed code for the linear elements are in good agreement with the results presented in the literature.
However, with the linear element, only some components of the mechanical strain gradient tensor can be calculated correctly. This means some components are nearly zero, even if they should have a nonzero value. The reason for this lies in the coefficients of the independently assumed mechanical strains. In the case of linear elements, the mechanical strain gradient components η iij are always zero because by taking the directional derivative of ε ij along x i only components of the P-vector remain, the coefficients of which are zero. The idea is to get more accurate results by extending CMFEM by using quadratic shape functions for mechanical displacement and electric potential so that the element can capture all mechanical strain gradient components.

Comparison of Linear and Quadratic Elements
At this point, it is essential to compare the performance of linear and quadratic elements with relevant boundary value problems (BVP). The numerical results are analyzed and compared with those in the corresponding literature. In this subsection, only those strain gradient components are considered which are nonzero in the linear elements. Other mechanical strain gradient components are hardcoded as zeros in the UEL for the quadratic element. From this, we can make the closest comparison with the known linear element behavior (cf. Section 4.1.1). The full calculation without such restrictions are presented in Section 4.1.4.
In Figures 9 and 10, the displacement u 3 and electric field E 3 are presented for different flexoelectric coefficients. Both quantities are measured at x 3 = H /2 along the length of the beam. It can be seen that with an increasing value of f 1 , the mechanical displacement reduces and the electric field increases. Although there is a negligibly small difference in the results of the linear element and the analytical solution, the results of the quadratic element appear to be in very good agreement with the linear element.

Convergence Study
A convergence study is made to compare the linear and quadratic elements. The relative error between numerical results for developed linear and quadratic elements and the analytical solution is analyzed for different mesh sizes. For this convergence study the same parameters for the cantilever beam problem is used as described in Tables 1 and 2 with f 1 = 0.1 µC /m and f 2 = 0.
The relative error e(X) of the numerical simulations is calculated by using the l 2 -norm X s represents the numerical results and X is the analytical solution. As a reference solution, the analytical results for the deformation of a Bernoulli-Euler beam are chosen where only the result for direct flexoelectricity is available. The deflection of the flexoelectric cantilever beam is given by [11] The exponent of L is corrected in [25] and the correct sign is taken from [11]: C 4 = − FL 2 6λ 2 g 11 I . The analytical solution for the generated electric field is given by [11] However, this analytical solution depends only on x 1 , whereas the numerical results are a function of both x 1 and x 3 . Consequently, the quadratic element is considered more accurate than the analytical solution [25]. For the convergence analysis, eight different meshes are analyzed, including 25, 100, 225, 400, 625, 900, 1225, and 1600 elements, respectively. Figure 11 illustrates the finest and roughest mesh used. Figure 11. Coarse (top) and fine (bottom) mesh used for the convergence analysis of the cantilever beam problem with total element number n. Not all elements are shown for better visualization. Figure 12 shows the relative error for the displacement u 3 and electric field E 3 for both element types with different average element length h, which is calculated from the total number of elements in the current mesh n as h = √ HL /n. It can be seen that the results for the linear elements are strongly mesh-dependent, and hence the relative error is higher for a mesh model with fewer elements. The relative error of the quadratic element is lower than the linear element, regardless of the coarseness of the meshed model. As presented, the accuracy of the results of the quadratic element is way higher than that of the linear element. This is expressed by the fact that the low accuracy of the linear elements when using a coarse mesh increases by increasing the fineness of the mesh, but even for a very fine mesh, it does not reach the performance of the quadratic elements for a very coarse mesh. The accuracy of quadratic elements, on the other hand, is the same for the coarsest and finest mesh.

Investigation of Differences between Linear and Quadratic Elements
For the analysis performed in this and further sections, the mechanical strain gradient components η iij are no longer set to zero as in Section 4.1.2. Furthermore, apart from a Poisson's ratio of ν = 0.3, the boundary conditions of Table 1 and material properties of  Table 2 are used.
The deflection and the electric field generated in the cantilever beam are presented for SGE with flexoelectricity f 1 = 0.1 µC /m in Figures 13 and 14, respectively. It can be observed that the compliance of the beam for the quadratic element is increased compared to the linear element. Due to oversimplification of the derivations for the linear elements, the electromechanical coupling is overestimated by the linear element and therefore the electric field for the linear element is higher compared to that of the quadratic element. Displacement of u 3 = 40 nm is prescribed for the cantilever beam to find the cause for the differences in the current investigation results. The resulting mechanical strains are presented in Figure 15 and is measured at x 1 = 0.1L along the height of the beam. It can be seen that mechanical strain ε 11 rises, ε 33 reduces and ε 13 is nearly the same. Strain ε 33 insignificantly differs at the edges for the linear elements compared to the quadratic element. The reason for this is explained in Section 4.1.6 (cf. Figure 27). The magnitude of all mechanical strain components is the same for the linear and quadratic elements. For the stress components, both types of elements gave similar results ( Figure 16). The reason for the small difference in stress components σ 11 , σ 22 and σ 33 at the edges is the same as mentioned above for strain component ε 33 . However, differences in calculated mechanical strain gradients using linear and quadratic elements are illustrated in Figure 17. For comparison reasons, in these cantilever beam simulations the material parameters of [25] are used where only the influence of f 1 is considered which means that f 2 is set to be zero. Consequently, according to Equation (68), the mechanical shear strains have in this example no direct influence on the magnitude of the generated electric field E 3 . Therefore, only the gradients of ε 11 , ε 33 are considered for both element types in the following analysis. In this numerical example, among all the components of mechanical strain gradients, only the two components η 311 and η 333 are of relevant magnitude. Strain gradient η 311 is a consequence of the distribution of mechanical strain ε 11 that results from the bending of the beam caused by the prescribed displacement u 3 . When u 3 is prescribed, the bending of the beam is the same for the linear and quadratic element, so η 311 is identical in both cases. Because of the usage of the nonzero Poisson's ratio, mechanical strain ε 33 is built up that varies along x 3 inside the cantilever beam, as shown in Figure 18. This results in a strong mechanical strain gradient η 333 at the beam's left side, which reduces toward the free end. In contrast to the quadratic element, the linear element is unable to capture this mechanical strain gradient component. A comparison of the mechanical strain gradient η 333 distribution is presented in Figure 19. Linear, σ 11 Linear, σ 33 Linear, σ 22 Linear, σ 13 Quadratic, σ 11 Quadratic, σ 33 Quadratic, σ 22 Quadratic, σ 13  The inability of the linear element to capture this mechanical strain gradient yields comparatively different higher-order stress in contrast to quadratic elements visible in Figure 20. The higher-order stress τ is measured along the length at the midheight of the beam. It can be observed that both higher-order stress components τ 311 and τ 333 reduce for the quadratic element, in particular τ 333 nearly vanishes. The components τ 111 and τ 133 remain zero. This reduces the beam's stiffness, which explains the beam's higher deflection for the quadratic element ( Figure 13). Based on Equation (15c), the higher-order stress exclusively depends on the contributions from flexoelectrical coupling and pure SGE as where the flexoelectric contribution is denoted with "f " and that from SGE with "g". For τ 311 and τ 333 , flexoelectric contributions are and the SGE contributions can be written as The flexoelectric contribution to higher-order stress is the same for τ 311 and τ 333 when f 2 = 0, cf. Equation (66). τ f 311 and τ f 333 reduce proportionally to the electric field E 3 as seen in Figure 22. Furthermore, the SGE contribution to higher-order stress differs for linear and quadratic elements because of the nonzero component η 333 in the quadratic elements (cf. Equation (67)). As a result, the contribution of SGE and flexoelectricity effectively cancel each other for higher-order stress τ 333 as illustrated in Figure 20. The resulting change in electric displacement that is measured along the height at x 1 = 0.1 L can be seen in Figure 23. For the quadratic element, electric displacement D 1 reduces while D 3 is zero for both element types along the height of the beam. The reason, for D 3 to come out to be zero, is that at the top and bottom side of the cantilever beam no electric potential is prescribed. If D 3 is zero and there is no piezoelectricity involved, Equation (15b) simplifies to the following formula for the electric field E 3 : The resulting electric field for the quadratic element is drastically reduced compared to the linear element. The reason for this is that the mechanical strain gradient component η 333 can be captured with quadratic elements and has a nonzero value. Additionally, η 333 has a sign opposite to η 311 , as can be seen in Figure 17. This explains the difference in intensity of the electric field in Figure 14. In this numerical example f 2 is omitted as in [25].

Parameter Study of Flexoelectric Coefficients
In this subsection, the influence of flexoelectric coefficients in case of the cantilever beam problem described in Section 4.1.1 is studied. The boundary conditions of Table 1 and material properties of Table 2 with Poisson's ratio ν = 0.3 are used. For both element types, two simulation series are performed in which the flexoelectric coefficient f 1 varies between 0.01 and 50.0 µC /m. For one case scenario, f 2 = 0 as commonly supposed in the literature [25,32,33], whereas for the other, f 1 and f 2 are equal. Figure 24 shows the maximum displacement u 3 for different flexoelectric coefficients. The flexoelectric coefficient axis has a logarithmic scale to cover a wider span of flexoelectric coefficients. One can see that for the linear element the flexoelectric influence is negligible when f 1 and f 2 are smaller than 0.02 µC /m. From f 1 = 0.05 µC /m, the influence of flexoelectricity on the beam's stiffness increases until at f 1 = 2 µC /m the beam becomes so stiff that the maximum deflection saturates nearing zero. When f 1 and f 2 are equal, this saturation value is reached for even smaller flexoelectric coefficients. On the other hand, the quadratic element's maximum deflection reaches a saturation value at 39 nm when f 2 = 0. When f 2 = f 1 , the value of the flexoelectric coefficients has negligible influence on the beam's stiffness.
The electric field for different values of the flexoelectric coefficient f 1 , measured at the fixed end, is presented in Figure 25. It can be seen that in the area of the small flexoelectric coefficients, the magnitude of the electric field increases alongside the stronger flexoelectric coupling. However, the electric field magnitude reaches a maximum value and with further increased flexoelectric coefficients, the generated electric field inside the beam gets closer to zero again.
Because from a physical standpoint the flexoelectric coefficients describe the ability of the center atom in a unit cell to break the symmetry and to generate an electric field as a result of the application of strain gradients, it is expected that for higher flexoelectric coefficients this center atom can move easier. From a numerical standpoint, strain gradient components η 311 and η 333 arise from the bending of the cantilever beam. However, if Poisson's ratio is zero [12,16,25], η 333 is not generated. When SGE is involved, η 311 and η 333 build up higher-order stress τ 311 and τ 333 , as stated in Equation (67). Due to flexoelectricity, both strain gradient components contribute to the generation of an electric field E 3 (cf. Equation (68)). However, because of the two-way coupling of flexoelectricity, this generated electric field E 3 again influences τ 311 and τ 333 , which affects the strain gradient components η 311 and η 333 . As can be seen in Equation (66), the pure flexoelectric contribution to higherorder stress components τ 311 and τ 333 is identical when f 2 = 0. When considering that for high flexoelectric parameters the pure flexoelectric contribution to higher-order stresses becomes more dominant compared to the SGE contribution, this leads to the corresponding strain gradients in a way that the magnitude of both strain gradients η 311 and η 333 become closer in magnitude. As can be seen in Figure 17, η 311 and η 333 have the opposite sign, which leads to their sum nearing zero for higher flexoelectric coefficients. With C 1 -continuity, η 333 not only influences the strains inside each element separately, but also influences the other results of the cantilever beam. This would lead to an electric field which saturates at a specific value. When the strains are not continuous from element to element, the magnitude of η 333 due to the presence of E 3 is overestimated, as discussed in Section 4.1.6. This leads to the sum of η 311 and η 333 nearing zero at a faster rate than the flexoelectric coefficients increase. According to Equation (68), this leads to a reduction of electric field which explains why E 3 is nearing zero for high flexoelectric coefficients. In case of C 1 -continuity, this behavior would be more balanced out in a way that the electric field reaches a saturation value. Because by using the linear element, η 333 is always zero, for high flexoelectric parameters, the sum of both strain gradient components only reaches zero when η 311 becomes zero, too. Mechanical strain gradient η 311 describes the change of ε 11 along x 3 resulting from the bending of the cantilever beam as visualized in Figure 5. One can see in Figure 26 that in the linear case in which there is nearly no deflection at all (cf. Figure 24), the stiffness of the beam becomes so high that there is nearly no bending of the beam. In contrast, there is no difference in the normalized bending profile between low-and high-flexoelectric coefficients for the quadratic element. This means that η 311 is proportionally reduced compared to the quadratic element. To sum up this subsection, the linear element is not suitable for flexoelectric simulations because of its inability to capture all strain gradients, though for correct flexoelectric simulations, it is necessary to calculate all strain gradient components correctly. By using quadratic shape function for the DOFs and assuming the independent mechanical strain and electric field of a quadratic shape, the quadratic element can compute all strain gradient components. Thereby, it can simulate a physically more realistic flexoelectric behavior. However, the quadratic element can be further improved by ensuring C 1 -continuity as discussed in Section 4.1.6.

Discussion about the Advantages of Quadratic Element
As mentioned in Section 4.1.4, the linear and quadratic elements' behavior partially differs. This is because the quadratic element can capture all mechanical strain gradient components, which is a priori not possible with the linear element. A higher resolution of ε 33 is presented in Figure 27. The mechanical strain gradient η 333 behavior for linear and quadratic cases is represented by the slope of ε 33 for each element. For the linear element, the mechanical strain in each element is constant; hence, the mechanical strain gradient is zero, whereas for the quadratic element, the mechanical strain inside each element is no longer a constant. The gradient in each element is slightly higher than the overall mechanical strain gradient along the height of the beam. Consequently, the linear element is not suitable for the description of some relevant strain gradient components. The quadratic element is able to describe all of the components and only overestimates slightly those mechanical strain gradient components that the linear element does not capture. However, calculating the correct mechanical strain gradients is important because in coupled systems like flexoelectricity, incorrectly calculated strain gradients lead to drastically different results of mechanical displacement, electric potential, and their derived quantities. The quadratic element is capable of simulating a realistic flexoelectric material behavior. Further improvement can be achieved by ensuring C 1 -continuity.

Truncated Pyramid Compression Problem
Another typical example of modeling flexoelectricity is a truncated pyramid [8,34]. Because of the different sizes of the top and bottom surfaces, applied stress on the top surface distributes toward the bottom surface, generating mechanical strain gradients.
Similar to the cantilever beam example, this is a plane strain problem for direct flexoelectricity (q = 0, b 1 = b 2 = b 3 = 0). The truncated pyramid has a height h, a top surface of length h and a bottom surface of length 3h. The angle between the bases and the oblique surfaces is π /4. The bottom surface is fixed in all directions. For compressing the pyramid, a constant linear force F is prescribed on the top surface, as illustrated in Figure 28. Furthermore, the pyramid is electrically grounded on the top surface. Due to the different sizes of the top and bottom surfaces, strain gradients build up which generate an inhomogeneous electric potential φ distribution. Typically, opposite the grounded surface, a sensing electrode is attached. This electrode equals the electric potential on this surface to a constant value which is not known a priori. The boundary conditions as in [16] are summarized in Table 3.
The material properties used are shown in Table 4. All simulation results are obtained for a mesh model containing 2500 elements. The results in Figures 29 and 30 are obtained for h = 7.5 µm. The contour plots of the mechanical strain ε 33 and electric potential φ are shown in Figure 29 for the linear and quadratic element. It can be observed that the mechanical strain is identical for both elements. In contrast, the electric potential using linear element is reduced by approximately 230% compared to the quadratic element and the reference solution [35]. Not only are higher maximum values reached by using the linear element, but the nature of the electric field distribution is physically invalid (cf. Figure 29 bottom right), even leading to a negative electric potential at the bottom part of the pyramid.

Quadratic Element
Linear Element . Contour plots of the electric potential φ in the truncated pyramid by [16,35] vs. the use of the developed quadratic element. The color scale for the contour of the quadratic element is adjusted for a higher contrast image. Figure 30 shows the contour of the electric potential φ for the same simulation using the quadratic element as in Figure 29 but with an adjusted scale for comparison with the results obtained using B-spline method [16,35]. One can see that the nature of the distribution of the electric field for the quadratic element is in good agreement with the results of the B-spline formulation. However, small differences in the electric potential distribution can be observed. For example at midheight, the electric potential gradient in x 1 -direction from the center to the oblige sides is higher using quadratic elements. This phenomena is physically more realistic, as can be seen in the computational results of [36].
Furthermore, an effective electric field E eff = V /h can be calculated, which is the voltage difference between the top and bottom bases. Figure 31 shows the effective electric field for different sizes of the truncated pyramid. For this simulation series, the line force F mentioned in Table 3 is used while h is different for each simulation. For different pyramid sizes the mechanical strain distribution is similar, whereas for high values of h the mechanical strain gradients are smaller compared to small pyramids. Moreover, the electric field is smaller for bigger pyramid sizes and becomes negligible for h greater than 200 µm. This phenomena is in accordance with the size effect of flexoelectricity. It is evident that the results of the quadratic element and [16,35] are in good agreement. . Effective electric field E eff of the truncated pyramid for linear element, quadratic element and [16] for different geometric sizes.

Conclusions
In the manuscript, a second-order collocation-based mixed FEM (CMFEM) for flexoelectricty is established. First, CMFEM formulation for linear elements is presented and tested. Because first-order CMFEM has deficiencies and cannot realistically capture flexoelectric material behavior, a quadratic CMFEM element is proposed and implemented. Here, quadratic shape functions are used for the DOFs and, consequently, independent mechanical strains and electric field are assumed to be quadratic polynomials for the collo-cation. The new elements are implemented as Fortran user elements and tested on several boundary value problems.
The achieved results and drawn conclusions can be summarized as follows.
• A linear CMFEM element is developed and verified with numerical results for a 2D cantilever beam problem from the literature. • The correctness of the newly proposed quadratic element is verified through comparison with the linear element and examples from the literature. To make the closest comparison, the strain gradient components, which cannot be computed by the linear element, are hardcoded to zero in the quadratic UEL. • Based on a convergence study, it is shown that the quadratic element is more accurate than the linear element. In particular, for meshes with very few elements, the accuracy is drastically improved with the usage of the newly developed quadratic element. • By no longer hardcoding specific strain gradient components to zero in the quadratic UEL, CMFEM elements' performance is analyzed. There, it is presented that mechanical strain gradient η 333 cannot be captured by the linear element, but has a relevant magnitude, leading to a strong overestimation of flexoelectric coupling and mechanical stiffness. • The abovementioned findings are further investigated by a parameter study in which cantilever beam's deflection and generated electric field are studied for a wide range of flexoelectric coefficients. There, it is found that using the linear element, the cantilever beams' deflection nears zero for high flexoelectric coefficients and the remaining deflection is build-up only from shear stress. For the quadratic element, on the other hand, the deflection saturates at a nonzero value mimicking a physically realistic behavior. • In the flexoelectric truncated pyramid example, although the linear element produced physically invalid results, the results of the quadratic element are in a very good agreement with the literature. The newly proposed quadratic element yields realistic results compared to those obtained in experiments [8,34] and a distribution of the electric field similar to [12,16]. This is additionally outlined in a simulation series wherein the flexoelectric size effect is visualized.
By using the newly developed second-order formulation, it is possible to model flexoelectric behavior accurately with less computational cost than traditional mixed FEM. This opens the possibility to apply the proposed quadratic element to the computationally demanding nonlinear electromechanical problems. It should be mentioned that there is scope for further improvement by ensuring C 1 -continuity.

Data Availability Statement:
The source code together with the data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.