Characterizing Flexoelectricity in Composite Material Using the Element-Free Galerkin Method

This study employs the Element-Free Galerkin method (EFG) to characterize flexoelectricity in a composite material. The presence of the strain gradient term in the Partial Differential Equations (PDEs) requires C1 continuity to describe the electromechanical coupling. The use of quartic weight functions in the developed model fulfills this prerequisite. We report the generation of electric polarization in a non-piezoelectric composite material through the inclusion-induced strain gradient field. The level set technique associated with the model supervises the weak discontinuity between the inclusion and matrix. The increased area ratio between the inclusion and matrix is found to improve the conversion of mechanical energy to electrical energy. The electromechanical coupling is enhanced when using softer materials for the embedding inclusions.


Introduction
An energy harvester utilizing the electromechanical coupling effect has been applied in various applications, ranging from sensors [1][2][3] to biomedical devices [4,5] at both micro-and nano-scales [6][7][8]. The well-known electromechanical coupling effect, piezoelectricity, generates the electrical polarization under mechanical deformation only in the non-centrosymmetric material. Flexoelectricity is another type of electromechanical coupling, which describes the coupling between the electrical potential and strain gradient. Unlike piezoelectricity, flexoelectricity can exist in any dielectric material regardless of the material's inner structure. In the meantime, experimental studies [9,10] observed an unexpected giant flexoelectric effect in barium strontium titanate, which suggested that at the micro-/nano-scale, the flexoelectric effect outperforms the piezoelectric effect. These superior properties of flexoelectricity at a small scale have drawn extensive research attention towards the fundamental investigation of this electromechanical coupling from the molecular scale to the continuum scale. Atomic-level density functional theories [11][12][13] and molecular dynamics simulations [14,15] explored the mechanism of inducing flexoelectric polarization under deformation. However, continuum scale studies can lead to better understanding of the practical applicability of the flexoelectric property in designing sensors and actuators.
The characterization of flexoelectricity with existing numerical methods, like the Finite Element Method (FEM), requires remarkable modifications to handle the strain gradient term. Zhang et al. [16] used conventional FEM to evaluate the piezo-and flexo-electric effects by solving the partial differential equations for the displacement and electric potential in a decoupled manner. The conventional FEM works on the basis of C 0 continuity and is uncertain when higher order strain terms are present (challenging to hold C 1 continuity). Nanthakumar et al. [17,18] used the Mixed Finite Element Formulation (MFEF) to maximize the flexoelectric effect through topology optimization of barium titanate material. However, the MFEF method suffers from high computational cost [19,20]. Isogeometric Analysis (IGA) is another numerical method, which ensures the required C 1 continuity by employing the higher-order shape functions [21][22][23]; whereas, IGA strongly depends on the geometrical symmetry conditions to reduce the computational cost [24,25]. Abdollahi et al. [26,27] developed a meshless method to study the flexoelectric response of dielectric material both in cantilever beam and truncated pyramid shapes. Their study further suggested that the simplified analytical model is unable to capture the flexoelectricity in multi-dimensional geometries. These studies have used structures with a ceramic material such as barium titanate or strontium titanate. Applications of these materials have a limited scope due to the small operational strain and high-stress conditions [28].
Experimental studies revealed that composites made up of two or many materials are promising to support high stress operating conditions. For example, nano-indentation tests of a bilayer cantilever beam structure have shown good mechanical properties [29,30]. A study on the polymer-based composite with graphene oxide has shown enhanced piezoelectric performance with greater flexibility [31]. Synthesized polymeric composite with piezoelectric zirconia titanate material promised to exhibit good energy-harvesting solutions [32]. However, these findings critically depend on the piezoelectric properties of the material. On the other hand, it is possible to create the piezoelectric composites without using a material with piezoelectric properties. The flexoelectric fiber-reinforced composite [33,34] and multi-material-based flexoelectric composites [35] are some example studies in this direction. The flexoelectric composites develop large strain gradients near the material discontinuity when uniformly deformed, which generates the electrical polarization through flexoelectricity. However, the numerical modeling other than conventional FEM requires extra care when dealing with the material and geometry discontinuity between the composite constituents. Lagrange multipliers [36,37] and the global enrichment approach [38,39] are the few techniques to treat the material discontinuity. The applicability of these methods is cumbersome when using complex geometries in 2D structures [40,41]. The Extended Finite Element Method (XFEM) has been used for the level set technique to describe the weak discontinuity between the inclusion and matrix composite [42][43][44]. In addition, the level set technique also allows using the multiple randomly-distributed inclusions inside the matrix.
The current work employs the Element-Free Galerkin (EFG) method [45][46][47][48][49] along with the level set technique to characterize flexoelectricity in a composite. It is assumed that the composite is a combination of embedding matrix material and a non-piezoelectric material (inclusion). The possible generation of electrical voltage due to the compressive loading of the composite is investigated. This study is extended to find the dependence of the electrical voltage on the concentration of inclusions in the embedding matrix. The numerical model is validated with standard examples and then utilized to explore the flexoelectricity in a composite. Section 2 composes the details of the EFG. Numerical validation and example case studies are presented in Section 3. Section 4 concludes the work.

Simulation Method
The enthalpy density H for a dielectric solid with the piezoelectric and flexoelectric effect is [50][51][52]: where E i = −φ ,i is the electric field; φ i being the electric potential in the i direction; ij is the mechanical strain; C ijkl is the fourth-order elastic moduli tensor; e ikl is the third-order tensor of piezoelectricity; µ ijkl is the fourth-order flexoelectric tensor; and κ ij is the second-order dielectric tensor. Appendix A covers further details about the weak form of Equation (1) and possible boundary conditions. The resulting governing equations are solved for unknown variables displacement u and electric potential φ using Moving Least Squares approximation (MLS). Belytschko et al. [53][54][55] introduced the MLS approximation within the EFG. According to the MLS, the unknowns (u and φ) for the center node are approximated by the neighboring nodes in its support domain (Figure 1a). If the support domain encounters a material discontinuity, corresponding support nodes are enriched. The weak discontinuity across the material interface in the composite is modeled with the level set function with absolute sign distance function Ψ(x). Figure 1b shows the schematic illustration of the level set function and its first-order-derivative across the interface. The local approximation of displacement u and electric potential φ at location x are given as: where Φ(x) is the shape function at position x, u and φ are the nodal displacements and electric potential, and N and M are the number of support and enriched nodes, respectively. More details about Φ(x) are given in Appendix B. Hereafter, we write the approximation in a simplified form The derivatives for unknowns u and φ (Equation (2)) contribute by both support nodes (u std , φ std ) and enriched nodes (u enr , φ enr ), which are given as:

Support nodes
Center node The discrete equilibrium equations involving Equation (3) are given as: where: where e refers to element number. Further details about the individual B and H matrices are included in Appendix C. The Lagrange multiplier method is used for imposing mechanical and electrical boundary conditions. In this study, the plane strain condition is assumed. Material property matrices C, κ, e, µ are given as:

Numerical Examples
In this section, we validate the numerical model with benchmark problems: electromechanical characteristics of a cantilever beam under mechanical and electrical loading conditions and coupling effects in a mechanically-compressed truncated pyramid. Note that these problems are free from the material discontinuity. The obtained results are validated with the reported literature. The complete numerical model with material discontinuity and local enrichment is employed to estimate the characteristics of the composite material.  Table 1 reports the various material parameters [21].   (Figure 2a), and the electric potential is constrained to zero. Electromechanical coupling induces the electrical energy under point load deformation. The conversion from mechanical to electrical energy (k 2 e f f ) is defined as:

Cantilever Beam
The present model is simplified by assuming that the transversal piezoelectric (e = e 31 ) and flexoelectric (µ = µ 12 ) components are the only non-zero in Equations (12) and (13). The Poisson effect is also ignored. The results of this simplified model are compared with the analytical model derived by Majdoub et al. [50]. The analytical solution for k e f f is: where the normalized piezoelectric constant is: where k piezo is obtained by neglecting flexoelectricity coefficient µ in Equation (15). Figure 3a plots the comparison between the present model and the analytical solution, where h = − eh/µ is the normalized beam thickness for the open circuit mechanical loading condition. The variation between e and h from the present model agrees with the analytical solution from Equation (16). This proves that the present model correctly estimates the electromechanical coupling in a non-piezoelectric beam under bending deformation. The coupling between mechanical and electric energy is strongly mediated by the flexoelectricity due to the low piezoelectric constant of the material (refer to Table 1). Figure 3b shows the distribution of the generated electric potential under the flexoelectric effect.

Electrical Loading
In this section, we study a cantilever beam (Figure 2b) under pure electric loading, which served as an actuator. The bottom edge of the cantilever beam (Figure 2b) enforced an electrical loading of −20V and grounded the top edge. There is no external mechanical loading in the cantilever beam. Figure 4 displays the displacement and electric potential profiles. The electric field across the beam in the y-direction is plotted in Figure 5. The electric field gradients near the top (negative sign) and bottom (positive sign) surface of the cantilever beam are due to the effect of converse flexoelectricity. This phenomenon generates mechanical stress at the top and bottom surface and thus deforms the beam. The displacement and electric potential distribution on the cantilever beam are in good agreement with earlier findings based on IGA analysis [21].
x direction displacement (m) y direction displacement (m)

Compress Truncated Pyramid
In this section, we have extended the validation for the compression of a truncated pyramid ( Figure 6). Due to its different top (a 1 ) and bottom (a 2 ) edge lengths, the applied uniform force F generates different tractions at the top and bottom edges, which results in a longitudinal strain gradient. The top edge of the pyramid is grounded, and a uniform force of −6e6 N is applied to it. The aspect ratio h of the pyramid is set as 75 µm, and the bottom edge length is 225 µm. The remaining material parameters are listed in Table 1. Figure 7 shows the developed strain in the y-direction and the electromechanical coupling-induced electric potential profiles. The numerical values and the field distributions are in good agreement with earlier reports [21]. This represents that the present model accurately estimated the electromechanical coupling.

Flexoelectricity in Composite
In this section, we demonstrate the possibility of inducing electric polarization in a composite system. The composite system is a combination of an embedding matrix in a square shape and a non-piezoelectric material (inclusion) in a circular shape, as shown in Figure 8. The composite under a uniform mechanical loading induces electrical polarization due to the local strain gradients near the material discontinuity. The dielectric constant of the inclusion material is about 10% of the embedding matrix. It is assumed that the Young's modulus of the matrix material (E mat ) is lower than the Young's modulus of the inclusion material (E inc ). Three different Young's modulus ratios ( E inc E mat = 10, 100, 1000) are used to understand its influence on the energy transfer ratio between mechanical energy and electrical energy. The edge length of the square domain is L = 10 µm with center inclusion radius r 1 = 1.5 µm. The remaining parameters (Table 1) Figure 9 presents the strain and electric potential profile for the domain with a center inclusion. The non-uniform strain field near the inclusion boundary generates the polarization gradient and electric potential. Figure 10 plots the strain gradient profile: yy,y and xx,x along the horizontal and vertical center line of the square domain, respectively. In both directions, a high strain gradient is seen near the boundary of the inclusion, which is due to the different material toughness of the matrix and the inclusion. Since the inclusion material is a non-piezoelectric material, the induced electrical potential is a consequence of flexoelectricity. Therefore, the linear relationship between electrical potential and strain gradient generates the strong potential near the inclusion, as seen from Figure 9b.
In general, the volume percentage of the inclusions has a vast impact on the overall composite properties [56]. In order to investigate the volume percentage effect on the electromechanical coupling, we have repeated the simulations with many inclusions under the same loading condition. The number of inclusions is varied according to the area ratio from 0.5-2.5%. The area ratio is defined as the total area of inclusions divided by the area of a square domain. For each area ratio, several inclusion configurations (center location for inclusions) are examined, and later, the results are averaged. Note that, for inserting many inclusions, the radius is decreased to 0.4 µm. Figure 11a shows the strong electric potential near the inclusion boundary for a domain with 10 randomly-distributed inclusions. This is due to the non-uniform strain distributions near the inclusion boundary. A strong electric potential is noted in the region between the nearby inclusions. This corresponds to the interaction of the non-uniform strain fields around the neighboring inclusions. Figure 11b shows the mesh configuration. Figure 12 plots the variation of the total electrical energy transfer rate with the area ratio. The energy transfer increases with the increase of the inclusion area ratio. This behavior is expected since higher irregularities of strain are induced when many inclusions are present. Figure 12 also indicates that domains with an identical inclusion area ratio with higher E inc E mat induce higher energy conversion from mechanical to electrical energy. This suggests that softer matrix materials enhance the electromechanical coupling characteristics.

Conclusions
In this study, the proposed model demonstrates the possibilities of inducing electromechanical coupling in a nano-composite material without the presence of the piezoelectric effect. The electromechanical coupling is modeled using the EFG method. Strain gradient involving partial differential equations is numerically solved using the MLS approximation. The C 1 continuity due to the strain gradient term is fulfilled by choosing the special weight function in the EFG. The level set techniques are employed to handle the weak discontinuity between the inclusion and matrix. The present numerical model captures the finite size flexoelectric effect for a mechanically-loaded cantilever beam. The converse flexoelectric effect-induced mechanical deformation of an electrically loaded cantilever beam is in good agreement with earlier reports. The results of the truncated pyramid further validate the current model. The non-uniform strain fields near the inclusion boundary induce the electrical polarization and thus electric potential due to flexoelectricity. We also found that the magnitude of the electromechanical coupling is largely dependent on the area ratio between the inclusion and the matrix. The higher the inclusion area ratio, the stronger the electromechanical coupling. Furthermore, a softer matrix material can also enhance the electromechanical coupling when compared to a stiff material. These findings help in designing a nano-composite utilizing the flexoelectric effect.

Appendix A. Theory of Flexoelectricity
The enthalpy density H for a dielectric solid with piezoelectric and flexoelectric effects is [50,51]: where E i = −φ ,i is the electric field; φ i is the electric potential; ij is the mechanical strain; C ijkl is the fourth-order elastic moduli tensor; e ikl is the third-order tensor of piezoelectricity; f iljk and d ijkl are the fourth-order direct and converse flexoelectric tensors; and κ is the second-order dielectric tensor. Sharma et al. [52] defined the flexoelectric tensor µ ijkl as the difference between d iljk and f ijkl by the application of integration by parts and Gauss divergence theorem to Equation (A1), which is: Rewriting Equation (A1) leads to Equation (1). The stress (σ ij ) and electric displacement (D i ) from Equation (1) when considering the piezoelectricity effect is given as:σ Due to the presence of flexoelectricity, the higher-order stress (σ ijk ) and electric displacement (D ij ) read:σ The total stress and electric displacement from piezoelectric and flexoelectric effects are summarized as: The essential and natural electric boundary condition are: whereφ and w are the applied electric potential and surface charge density, ∂Ω represents the boundary of the domain, and n i is the unit normal to the boundary ∂Ω. The mechanical boundary conditions are given as: whereū andt k are prescribed mechanical displacement and traction. The remaining boundary conditions (normal derivation of displacement and higher-order tractions) resulting from the strain gradient have been set to zero under the assumption of a homogeneous condition. Rewrite Equations (A3) and (A4) as: and integration over the domain Ω gives: where H is the total electrical enthalpy. The external work by surface mechanical and electrical forces is given by: Finally, the weak form of the mechanical and electrical equilibrium derived from the Hamilton principle for the static problem yields: Substituting Equations (A3)-(A5) into Equation (A11) yields: The unknowns (displacement and electric potential) in Equation (A12) are approximated using MLS.

Appendix B. Details of the Shape Function
The shape function Φ I (x) associates with Node I and a point x under MLS is: where p(x) is the second-order polynomial, which is: The quadratic spline weight function w ensures C 3 continuity inside an element and C 2 continuity between elements [40]. The mathematical expression for w is: where: d is the predefined search radius of the support domain and d equals three-times the nodal spacing. The moment matrix A(x) has the form: A sufficient number of support nodes ensures the non-singularity of matrix A(x). The enrichment function Ψ(x) has the form [42]: Ψ(x) = abs(ψ(x)); where ψ(x) = min where n c is the total number of inclusions inside the domain, x i c is the center coordinate of the i th circular inclusion, and r i c is the radius of the i th circular inclusion.