Research on Output Characteristics of Microscale BST Laminate Structure Based on Mixed Finite Element Method

The flexoelectric effect, which is sensitive to size, refers to the phenomenon of coupling between the strain gradient and electrical polarization and involves higher-order derivatives of physical quantities such as displacement, and the analytical process is complicated and difficult. Therefore, in this paper, a mixed finite element method is developed considering the effects of size effect and flexoelectric effect on the electromechanical coupling behavior of microscale flexoelectric materials. Based on the theoretical model of enthalpy density and the modified couple stress theory, the theoretical model and finite element model of microscale flexoelectric effect are established, and the Lagrange multiplier is used to coordinate the higher-order derivative relationship between the displacement field and its gradient, and the C1 continuous quadrilateral 8-node (displacement and potential) and 4-node (displacement gradient and Lagrange multipliers) flexoelectric mixed element. By comparing the numerical calculation results and analytical solutions of the electrical output characteristics of the microscale BST/PDMS laminated cantilever structure, it is proved that the mixed finite element method designed in this paper is an effective tool for studying the electromechanical coupling behavior of flexoelectric materials.


Introduction
In the 1960s, people discovered the non-traditional "piezoelectric behavior" in nonuniformly deformed materials-the flexoelectric effect, that is, the electrical polarization of the dielectric material induced by the strain gradient change (positive flexoelectric effect) and caused by the electric field gradient. It produces a phenomenon of strain response (inverse flexoelectric effect) [1]. Different from piezoelectric materials, flexoelectric materials do not require pre-polarization, have no temperature limit, and have scale sensitivity. They have wide application prospects in the fields of micro-and nano-scale energy harvesting, sensing and driving, and have attracted much attention from researchers [2]. Some achievements have been made in the preparation of flexoelectric materials [3], performance testing and characterization [4], theoretical modeling of electromechanical coupling effects [5], numerical simulation methods [6], and application exploration [7,8].
The non-locality of the strain gradient or polarization (electric field) gradient makes the boundary value problem with respect to the flexoelectric effect of a solid dielectric material dominated by a fourth-order partial differential equation, which is only applicable when models such as one-dimensional or axisymmetric models (Simple geometric problems such as cylinders, disks or plates) can be solved analytically, while the flexoelectric response problems of those complex structures need to be solved by corresponding numerical methods. It is only in the process of numerical solution that the higher-order partial differential equations require that each element meet the higher-order continuity condition; that is, the shape function not only needs to satisfy the function continuity on the edge of 2 of 11 the element but also the first-order derivative function of the function needs to satisfy the continuity function.
Shu et al. [9] established a mixed 2D finite element with displacement gradients and Lagrangian multipliers as extra-nodal degrees of freedom to analyze the electrical polarization generated in a single crystal flexoelectric cantilever beam and used the C 0 continuous shape. The function achieves the same convergence as the C 1 continuous function. Amanatidou and Aravas [10] developed some mixed elements, which can obtain relatively accurate solutions in some classical problems (such as single-crystal cantilever beams), but the theory requires additional material intrinsic parameters that are difficult to determine experimentally. Couple stress theory [11] is a method for studying microscale effects, which has the advantage of having few additional material intrinsic constants and is relatively easy to determine experimentally. Regarding the general strain gradient elasticity [12] problem, Hutchinson [13] and Shu [14] took the couple-stress rotation angle as an additional nodal degree of freedom and solved the flexoelectric higher-order partial differential equation by the mixed finite element method, which verifies that the couplestress Feasibility of theoretically dealing with flexoelectric problems.
Aiming at the high-order continuity problem encountered in the numerical analysis of flexoelectric response, Abdollahi et al. [15] first used the Local Maximum Entropy (LME) meshless method to try to solve the modeling of flexoelectric single crystal cantilever structure. In order to simplify the model and reduce the difficulty of modeling, it deliberately ignores the influence of part of the flexoelectric effect on the boundary value problem. Based on the work of Amanatidou et al., Mao et al. [16] constructed a mixed formulation of flexoelectricity with displacement and displacement gradient as independent variables and developed a 2D cell to model BVPs (boundary value problems). However, the model they developed has an additional degree of freedom of the polarization node, which still has many problems that the constants are difficult to determine. Deng et al. [17] proposed a mixed finite element method with displacement, electric potential, displacement gradient and Lagrange multipliers as degrees of freedom to study the flexoelectric coupling response of dielectric materials and found that after considering the strain gradient, it can make the strain field smoother near the boundary. Nguyen et al. [18] used the IGA (Isogeometric analysis) method to solve the high-order continuum problem of the flexoelectric effect and focused on the influence of the flexoelectric effect on the potential of the nanowire, but the research on the scale effect at the micro-nano scale was not comprehensive enough. The above work has well promoted the numerical modeling research of flexoelectric effect, but the influence of size effect is ignored in the above research, and the size effect of flexoelectric materials is only studied based on the strain gradient theory, resulting in many parameters that are difficult to determine. In fact, with the continuous development of modern processing technology, the structure of micron-scale MEMS functional devices in practical applications is constantly emerging, and the preparation process is becoming more and more mature so that researchers must refocus on the research of micro-scale flexoelectric effect. BST ceramics are the preferred dielectric materials for the preparation of flexoelectric functional structures (the flexoelectric coefficient is much higher than other dielectric materials), and combined with flexible materials (such as PDMS) to form a flexoelectric composite/laminate structure, it can give full play to its power-to-electricity conversion performance can effectively avoid the disadvantage that the micro-scale structure is difficult to form and damaged.
In this paper, the modified couple stress theory is introduced into the flexoelectric enthalpy density theoretical model, and a mixed finite element method is developed to explore the flexoelectric effect by taking the cantilever structure laminated by microscale ceramic material BST and flexible material PDMS as an example. And the influence of scale effect on the output voltage of microscale flexoelectric laminate structure and the energy harvesting efficiency of cross-scale flexoelectric structure can be improved through the parametric design of the structure.

Electrical Enthalpy Density Theoretical Model of Microscale Flexoelectric Effect
Shen and Hu established the enthalpy density theory to describe the mechanical energy density, dielectric energy density, piezoelectric energy density and flexural electrical energy density of dielectric elastomers according to the energy method [19,20], and the specific form is as follows: where C ijkl is the elastic coefficient tensor, k ij is the dielectric constant, ε jk is the strain tensor, E k is the electric field vector, e ijk is the piezoelectric constant tensor, g ijkl , f klij are the positive and inverse flexoelectric coefficient tensors, respectively, ε ij,l are the strain gradient tensors, are the b ijkl is the non-local electric field coupling coefficient tensor. In order to accurately characterize the scale characteristics of dielectric solid materials at the microscale, a microscale factor term C 0 ijkl χ ij χ kl /2 is added on the premise of mutual verification with the experimental results, and the positive and inverse flexoelectric coefficients of the dielectric materials are considered objectively g ijkl = −f ijkl . Then a new enthalpy density function is obtained: Among them, C 0 ijkl = Gl 2 is the microscopic elastic stiffness matrix, l is the intrinsic constant of the material, which only depends on the microstructural characteristics of the material and is related to the microscale properties of the material, but has nothing to do with the specific structural size of the material, G is the Lame constant; the couple stress term χ ij represents that as the relevant scale gradually decreases, the ion polarization effect increases, which in turn affects the enthalpy density of the material; at the same time, the contribution of the orientation polarization effect, which accounts for a large proportion in the macroscopic scale, to the structural, electric enthalpy density decreases with the decrease of the scale, which is embodied in the change of the macroscopic strain term ε ij .
Studies have shown that there is no piezoelectric effect in the centrosymmetric crystal point group dielectric solid materials, which can be made according to the stress tensor σ ij , high-order stress tensor τ jkl , electric displacement vector D i , couple stress tensor m ij , electric quadrupole moment Q ij and electric enthalpy density function [21], the micro-scale flexural-electrical coupling constitutive equation with only flexoelectric effect is obtained:

Constrained Variation Principle
For a two-dimensional plane problem, the curvature tensor and couple stress tensor are: where rotation vector θ 1 = θ 2 = 0, θ 3 = 1 2 ∂v ∂x − ∂u ∂y . On the boundary it should satisfy: Integrate the inner term and boundary term of the constitutive relation of microscale flexoelectric materials, and integrate E = −∇ 1 ϕ and χ = ∇ 1 θ, the enthalpy function of the two-dimensional enthalpy model can be obtained: where p is the surface force vector, d is the surface force couple, and b is the charge surface density. For a two-dimensional plane problem, the displacement gradient vector can be expressed as: Then the corresponding strain, strain gradient, rotation angle and rotation gradient can be expressed as a functional form with displacement gradient: Substitute Equation (8) into Equation (6) and introduce the constraint relationship between displacement and displacement gradient Equation (7), set the Lagrange multiplier to be α, then the induction enthalpy function becomes: The formula obtained by dividing the enthalpy density function after the variation is valid for any δu, δϕ, δψ and δα, then we can get: The boundary terms must satisfy: (11) Equation (11) is the Navier equation of the flexoelectric effect in the microscale flexoelectric structure, and its boundary conditions are: u = u, on u; ϕ = ϕ, on ϕ; ψ = ψ, on ψ (12) where u, ϕ and ψ are the known displacement, potential and displacement gradient on the boundary. For the BST material of the m3m point group, the forward flexoelectric coefficient matrix is: For the isotropic m3m point group crystalline material, the elastic coefficient matrix C 0 is: The dielectric coefficient matrix κ, and the specific form of the operator tensor that appears in the calculation process is as follows: Figure 1 shows the mixed element of quadrilateral four-nodes and quadrilateral eightnodes selected in this paper. The dielectric coefficient matrix κ , and the specific form of the operator tensor that appears in the calculation process is as follows: Figure 1 shows the mixed element of quadrilateral four-nodes and quadrilateral eight-nodes selected in this paper. The displacement is 8-node interpolation, and each element node has two degrees of freedom: (16) where (8) i N represents the 8-node shape function.

Microscale Flexoelectric Mixed Cell Construction
The potential also selects an 8-node interpolation function, and each element node has 1 degree of freedom: The displacement gradient selects 4-node interpolation, and each element node has 4 degrees of freedom: The displacement is 8-node interpolation, and each element node has two degrees of freedom: represents the 8-node shape function. The potential also selects an 8-node interpolation function, and each element node has 1 degree of freedom: The displacement gradient selects 4-node interpolation, and each element node has 4 degrees of freedom: (18) where N (4) i denotes a 4-node shape function. The Lagrange multiplier selects 4-node interpolation, and each element node has 4 degrees of freedom: The constraints between the corresponding basic variables are: (20) According to Equations (13)- (20), it can be deduced that the equivalent strain matrices B ε , B E , B η , B u and B χ respectively are:

Output Voltage Analysis of Single Functional Layer BST/PDMS Laminated Structure
On the basis of the previous section of this article, the mechanical and electrical coupling behavior of the micro-scale flexoelectric laminated cantilever beam with a flexible layer attached to the surface is explored by following the constitutive relationship of the micro-scale flexoelectric material. The PDMS flexible layer attached to the beam surface can effectively reduce the fabrication difficulty of micro-scale flexoelectric composite beams and play a protective role in static and dynamic tests. Figure 2 shows a finite element schematic diagram of a microscale flexoelectric laminated cantilever beam structure with a fixed left end and a free right end. The cross-section of the laminated beam is rectangular, with a BST flexoelectric layer in the middle, and its upper and lower surfaces are coated with a high-temperature electrode layer with a negligible thickness, and a PDMS flexible layer is attached to the outside of the electrode layer. The beam length is 2 cm, the BST layer thickness is h B , and the flexible layer thickness is both h P /2 = 5 µm, intrinsic size characteristic parameters of the material is l = 5 µm. The material parameters of the laminated beams are shown in Table 1.  , and the displacement between the structural layer elements is used as the constraint boundary condition. The upper and lower surfaces of the BST sensitive layer are equipotential surfaces. Figure 3 shows the variation of the output voltage of the microscale flexoelectric laminated beam with the thickness of the BST layer under different static loads calculated by the mixed finite element method and the analytical method, where the voltage represents the output voltage in the open-circuit environment. As can be seen from the figure: ① The agreement between the analytical solution and the numerical solution is relatively high. The force couple, uniform load and concentrated force reach their peaks at the thickness of the BST layer of 4.3 µm. At this time, the errors are 0.7%, 0.5% and 1.1%, respectively, which verifies the finite element method. accuracy and universality; ② The variation trend of the voltage output of the laminate structure under the three loads remains the same. The thickness of PDMS remains unchanged, and the output voltage of the BST/PDMS laminate structure first increases and then decreases with the increase in the thickness of the BST layer. This is due to the combined effect of the increase in the thickness of the BST, the increase in the specific gravity of the mechanoelectrical conversion functional layer of the structure and the decrease in the size effect (strain gradient).      1 The agreement between the analytical solution and the numerical solution is relatively high. The force couple, uniform load and concentrated force reach their peaks at the thickness of the BST layer of 4.3 µm. At this time, the errors are 0.7%, 0.5% and 1.1%, respectively, which verifies the finite element method. accuracy and universality; 2 The variation trend of the voltage output of the laminate structure under the three loads remains the same. The thickness of PDMS remains unchanged, and the output voltage of the BST/PDMS laminate structure first increases and then decreases with the increase in the thickness of the BST layer. This is due to the combined effect of the increase in the thickness of the BST, the increase in the specific gravity of the mechanoelectrical conversion functional layer of the structure and the decrease in the size effect (strain gradient).    From the variation law of the output voltage of the structure with BST thickness shown in Figure 5, it can be seen that the numerical solution is in good agreement with the analytical solution; as the thickness of the BST layer increases, the output voltage of the overall structure increases first and then decreases; The output voltage of the composite beam structure reaches the maximum value when the thickness of the BST layer is 3 From the variation law of the output voltage of the structure with BST thickness shown in Figure 5, it can be seen that the numerical solution is in good agreement with the analytical solution; as the thickness of the BST layer increases, the output voltage of the overall structure increases first and then decreases; The output voltage of the composite beam structure reaches the maximum value when the thickness of the BST layer is 3 µm, and the error between the numerical solution and the analytical solution is about 0.7%, which is within the acceptable range.  Setting the thickness of the BST layer (3 µm) and the uniform load 100 μN/m q = − applied to each BST layer, the variation law of the output voltage of the BST/PDMS flexoelectric laminated cantilever structure with the number of BST layers under different PDMS layer thicknesses is shown in Figure 6. Setting the thickness of the BST layer (3 µm) and the uniform load q = −100 µN/m applied to each BST layer, the variation law of the output voltage of the BST/PDMS flexoelectric laminated cantilever structure with the number of BST layers under different PDMS layer thicknesses is shown in Figure 6. Setting the thickness of the BST layer (3 µm) and the uniform load 100 μN/m q = − applied to each BST layer, the variation law of the output voltage of the BST/PDMS flexoelectric laminated cantilever structure with the number of BST layers under different PDMS layer thicknesses is shown in Figure 6. It can be seen from Figure 6 that when the thickness of the PDMS layer is constant, the output voltage of each BST layer transducer element decreases with the increase in the number of layers of the microscale flexoelectric laminate structure; when the thickness of the PDMS layer is 20 µm (When the thickness of the BST layer is about 5 times), the output voltage of the structure decreases slightly with the increase in the number of BST layers. This is because the increase in the number of BST flexoelectric functional layers leads to an increase in the overall stiffness of the BST/PDMS laminate structure, which in turn leads to a decrease in the output voltage of each BST transducer element. It can be seen from Figure 6 that when the thickness of the PDMS layer is constant, the output voltage of each BST layer transducer element decreases with the increase in the number of layers of the microscale flexoelectric laminate structure; when the thickness of the PDMS layer is 20 µm (When the thickness of the BST layer is about 5 times), the output voltage of the structure decreases slightly with the increase in the number of BST layers. This is because the increase in the number of BST flexoelectric functional layers leads to an increase in the overall stiffness of the BST/PDMS laminate structure, which in turn leads to a decrease in the output voltage of each BST transducer element.

Analysis of Output Voltage of Multilayer BST Laminated Structure
Selected h P = 5h B , the output voltage of the BST/PDMS laminate structure with different layers of BST sensitive layers under the action of uniform load is shown in Figure 7. It can be seen that the output voltage of the BST/PDMS laminated structure with different BST layers increases first and then decreases with the increase in the thickness of the BST layer, and the variation law is consistent; with the increase in the thickness of the BST layer, the corresponding single-layer BST thickness at the peak of the output voltage of the structure is smaller. As the number of layers increases gradually, the weakening amplitude of the influence of the size effect on the output voltage of the structure after the thickness of the overall structure increases is smaller than that of the flexoelectric effect, and a sufficient number of superimposed layers is sufficient to compensate for the scale effect and the flexural effect. The energy loss caused by the electrical effect at this time, the thickness of the BST transducer layer can be appropriately increased when the laminated structure reaches the optimal output state, which can reduce the difficulty and loss of preparation in practical engineering applications to a certain extent.
Considering the difficulty of preparation of BST sensitive layer and the difficulty of attachment of too thin PDMS flexible layer in practical engineering applications, a BST sensitive layer with a thickness of 10 µm was selected as the basic element, and the PDMS layer was still 5 times the thickness of BST layer. The output voltage of the microscale flexural electric laminated beam energy harvesting device with different layers and thicknesses is calculated by the finite element method, as shown in Figure 8. It can be seen that with the increasing number of layers of microscale flexoelectric lamination beam under a constant thickness, the output voltage of the overall structure shows a linear growth trend, which proves that it is feasible to construct an energy harvesting device with both microscale effect and energy output characteristics with the thickness of BST and PDMS.
plitude of the influence of the size effect on the output voltage of the structure after the thickness of the overall structure increases is smaller than that of the flexoelectric effect, and a sufficient number of superimposed layers is sufficient to compensate for the scale effect and the flexural effect. The energy loss caused by the electrical effect at this time, the thickness of the BST transducer layer can be appropriately increased when the laminated structure reaches the optimal output state, which can reduce the difficulty and loss of preparation in practical engineering applications to a certain extent. Considering the difficulty of preparation of BST sensitive layer and the difficulty of attachment of too thin PDMS flexible layer in practical engineering applications, a BST sensitive layer with a thickness of 10 µm was selected as the basic element, and the PDMS layer was still 5 times the thickness of BST layer. The output voltage of the micro-scale flexural electric laminated beam energy harvesting device with different layers and thicknesses is calculated by the finite element method, as shown in Figure 8. It can be seen that with the increasing number of layers of microscale flexoelectric lamination beam under a constant thickness, the output voltage of the overall structure shows a linear growth trend, which proves that it is feasible to construct an energy harvesting device with both microscale effect and energy output characteristics with the thickness of BST and PDMS.

Conclusions
Based on the Lagrange multiplier method, the displacement and displacement gradient are constrained, and the C 1 continuity of the element boundary is taken into account, a new mixed element is constructed. According to the parameter transformation method, a two-dimensional mixed finite element method related to the size effect and flexoelectric effect is developed.
By adjusting the proportion of BST layer and PDMS layer in the micro-scale flexoelectric laminated beam, the output voltage of the laminated beam structure with different structural sizes under static load was studied, and a three-layer micro-scale flexoelectric laminated beam structure was simulated to calculate the output voltage of the whole structure under uniform load, which was compared with the analytical solution and found that the error was very small. The effectiveness and universality of the finite element method are verified.
Finally, through the simulation of the multilayer microscale flexoelectric energy harvesting device, the influence trend of the scale effect and flexoelectric effect on the output voltage of the laminated beam in the multilayer structure state is compared, and the convenience of the finite element method can be seen by comparing with the theoretical analysis method.

Conclusions
Based on the Lagrange multiplier method, the displacement and displacement gradient are constrained, and the C 1 continuity of the element boundary is taken into account, a new mixed element is constructed. According to the parameter transformation method, a twodimensional mixed finite element method related to the size effect and flexoelectric effect is developed.
By adjusting the proportion of BST layer and PDMS layer in the micro-scale flexoelectric laminated beam, the output voltage of the laminated beam structure with different structural sizes under static load was studied, and a three-layer micro-scale flexoelectric laminated beam structure was simulated to calculate the output voltage of the whole structure under uniform load, which was compared with the analytical solution and found that the error was very small. The effectiveness and universality of the finite element method are verified.
Finally, through the simulation of the multilayer microscale flexoelectric energy harvesting device, the influence trend of the scale effect and flexoelectric effect on the output voltage of the laminated beam in the multilayer structure state is compared, and the convenience of the finite element method can be seen by comparing with the theoretical analysis method.
Author Contributions: Conceptualization, Y.L. and H.L.; methodology and modelling simulation, Y.L. and T.P.; data processing, T.P. and H.L.; original draft preparation, T.P.; writing-review and editing, Y.L., H.L. and T.P.; supervision, Y.L. All authors have read and agreed to the published version of the manuscript.
Funding: This research is funded by the National Science Foundation of China, grant number 12072133.

Data Availability Statement:
The data presented in this study are contained within the article. No additional data were created or analyzed in this study.

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