Asymptotic Homogenization Applied to Flexoelectric Rods

In this manuscript, the equilibrium problem for a flexoelectric one-dimensional composite material is studied. The two-scales asymptotic homogenization method is used to derive the homogenized formulation of this problem. The manuscript offers a step-by-step methodology to derive effective coefficients and to solve local problems. As an illustrative example, results reported in the literature for piezoelectric composites are obtained as a particular case of the formulation derived here. Finally, three flexoelectric/piezoelectric composites are studied to illustrate the influence of the flexoelectric property on the effective coefficients and the global behavior of the structure.


Introduction
In the last 35 years, the study of piezoelectric composites has increased because of their use in several areas of engineering [1], contributing to the development of new mechanical structures that are used in high-technology devices [2,3]. The wide range of applications of piezoelectric composites has advanced the development of mathematical, experimental, and computational models related to the study of the properties of these materials, see for instance References [4][5][6]. Piezoelectric composites have been extensively studied, but only in very recent years, flexoelectric materials have taken an important place in the new scientific works, see Reference [7].
Piezoelectricity is usually expressed as an interaction between mechanical strain and one of the electrical variables: the electric field, the electric displacement, or the electric polarization. In Reference [8], the author examines the consequences of considering an additional, linear, electromechanical effect: an interaction between the strain and the polarization gradient. Experiments have shown that when large amplitude mechanical disturbances propagate through a dielectric medium a voltage is developed across the ends of the sample, see Reference [9]. Some of these experiments with nonpiezoelectric elements have been shown to produce a polarization charge upon shock loading. Strain-gradient-induced polarization is known as the flexoelectric effect. According to Reference [10], flexoelectricity is an electromechanical effect of dielectric materials whereby they exhibit a spontaneous electrical polarization induced by a strain gradient (inhomogeneous deformation).
The flexoelectric effect can be considered as a high-order electromechanical phenomenon with respect to the piezoelectric effect [7].
Flexoelectricity was first theoretically predicted in the 1950's and described from a phenomenological standpoint in the 1960's [11]. Recently, flexoelectric materials have gained prominence due to the development of new methods for manufactured materials with coupled mechanical and electrical behavior [11] and their applications in membrane structures [12][13][14].
In References [15,16], important contributions to the study of flexoelectric materials have been presented highlighting the differences between piezoelectric and flexoelectric structures. The authors emphasize the significance of the flexoelectric effect discussing several applications which include materials with large flexoelectric coefficients as well as cases in which electromechanical coupling of piezoelectricity is not present and in applications related to soft materials. Moreover, some important applications of the flexoelectric effect are focused on the study of biological membranes and the development of piezoelectric structures without using piezoelectric materials, sensing, actuating, or energy harvesting.
Many authors have studied the flexoelectric mechanical equilibrium equations (for instance References [11,17]), bringing together different approaches to solve the problem at hand. One of the most common methods in solid mechanics to approximate the solution of the equilibrium problem of a solid (elastic, piezoelectric, flexoelectric, etc.) is the two-scales asymptotic homogenization method (AHM) [18,19]. The multiscale asymptotic homogenization method has been widely used to derive the effective properties of composite materials with different mechanical properties. In References [18,[20][21][22], the AHM is used to derive the effective properties of elastic composites (laminates, fibrous, wavy, etc.). These ideas are extended to the study of piezoelectric structures [23]; viscoelastic composites [24]; and thermo-piezoelectric materials [25]. To the best of the authors knowledge, a methodology to find the effective properties of flexoelectric materials has not been presented before. For that reason, this work acts to extend those methodologies (that find effective properties) to include the case of flexoelectric materials using the two-scales asymptotic homogenization method, which provides the first step to subsequently study three-dimensional curvilinear flexoelectric structures.
In general, the two-scale asymptotic homogenization method is used to find the effective properties of periodic composites, i.e., for any mechanical property p(x), x ∈ R n , there is a T ∈ R n s.t. p(x + T) = p(x) for all x ∈ R n . For non-periodic media, statistical and numerical methods have been developed to find representative volume elements and the associated effective properties [26][27][28]. In Reference [29], the self-consistent method in explicit form (effective field) and implicit form (effective medium) as well as the AHM are used to derive the effective piezoelectric properties of fibrous composites with randomly positioned fibers and periodically distributed fibers, showing good agreement among the different methods.
The work is organized as follows. In Section 2, a theoretical framework for the equilibrium of a one-dimensional periodic flexoelectric composite material is presented. The asymptotic expansions of the displacement and the electric potential are introduced in Section 3, using the two-scales asymptotic homogenization method, and a methodology to derive the expression of the effective coefficients depending on the so-called local function expressions is illustrated. The process to obtain the analytic solutions of the local problems is shown. A variational formulation of the local problem is proposed and the system is solved using the classical finite element method (FEM). As an important benchmark, the model reproduces the solution of the local problems reported in Reference [30] for a one-dimensional piezoelectric structure. Finally, in Section 4, the theoretical results derived in the previous sections are used to study some numerical cases. The effective coefficients reported in References [19,30] for a two-material piezoelectric structure made of lead zirconate titanate (PZT) C91 and P-82 are derived using the model developed here, which makes use of the solution of the local problem variational formulation (FEM). In order to extend the study for the case of randomly distribute one-dimensional flexoelectric composites, a two-constituents structure Ω N is considered for which the materials properties are randomly distributed along the rod with a binomial probability. Using the concept of correlation length for the particular case of one-dimensional binomial distributed structure, as reported in Reference [31], the methodology is then extended to the case of non-periodic structures. As it is shown in Reference [29], the effective properties of random composites converge to the effective properties of periodic structures when the length of Ω N increases. In addition, the combination of flexoelectric/piezoelectric/non-piezoelectric composites are studied in order to illustrate the influence of the flexoelectricity on the effective properties. One of the composites is considered as the combination of two flexoelectric materials, barium titanate (BaTiO 3 ) and strontium titanate (SrTiO 3 ) [32][33][34]. Other composites are combinations of these two flexoelectric materials with a exclusively piezoelectric material, polyvinylidene fluoride (PVDF) [35,36], a non-piezoelectric polymer described in Reference [37] and PZT-7A reported in Reference [38]. From these material combinations important conclusions are derived with respect to the influence of the flexoelectric property in the global behavior of the composites. Finally, a comparison between the solutions of a heterogeneous and homogeneous problem for a one-dimensional flexoelectric structure is presented. The finite element method is used to solve the problems of two different examples considering the same prescribed boundary conditions but different external forces.

Flexoelectric One-Dimensional Problem
A one-dimensional flexoelectric composite rod Ω is considered. The constitutive relations between stress σ and electrical displacement D with strain and electric field E are written in the following form (see Reference [17]), where C(x), e(x), κ(x), and µ(x) denote the stiffness, piezoelectric, permittivity, and flexoelectric tensors, respectively. Linear deformations of the composite, the strain, and the electric field are studied in terms of the displacement u and the electric potential φ as Assume the material parameters C, e, κ, µ are rapidly oscillating periodic functions along the rod with respect to the variable y, i.e., where y = x/ε ∈ Y is the local variable, since 0 < ε << 1 is a very small parameter that characterizes the periodicity of the structure, and Y is the periodic cell. The composite has the structure presented in Figure 1.
To simplify the equations, material parameters and functions are considered to depend on the position variable x or local variable y, i.e., u ≡ u(x), φ ≡ φ(x), C ≡ C(y), e ≡ e(y), κ ≡ κ(y), µ ≡ µ(y). The equilibrium equations presented in Reference [17] for one-dimensional solids take the following expression: with boundary conditions where the functions f and ρ are the external forces and body charge density, respectively. On the other hand, the values u 0 , S 1 , w 0 , r 1 , ϕ 0 , and τ 1 are the prescribed displacement, traction, displacement derivative, high-order traction, electric potential and charge, respectively. Perfect contact conditions at the interface are considered, i.e., function displacement u is differentiable and classical stress C + eE, higher order stress µ dE dx (see Reference [17]), electric potential φ and electric displacement D are considered continuous functions at the interfaces of the materials. The boundary value problem (5)-(8) is cast as a linear system of third order ordinary differential equations for which, under the perfect contact conditions as mentioned above, existence and uniqueness of the solution on the interval [0, 1] is guaranteed [39,40].

Asymptotic Homogenization Method
The boundary value problem (5)-(8) has rapidly oscillating coefficients. In order to approach the solution of the problem, an equivalent homogenized problem must be obtained. To derive the expression of the homogenized system, the two-scales asymptotic homogenization method is used. In Reference [19], the asymptotic expansion to the solution of the piezoelectric problem is reported as where v 0 (x) and φ 0 (x) only depend on the global variable x and the functions N k , Π k , Ψ k , and Θ k are ε-periodic continuous functions that only depend on the local variable y, i.e., N k (y + ε) = N k (y), . Now, the expansions (9)-(10) are introduced into the equilibrium Equations (5)-(6), for which one needs to consider the derivative operator where (·) x and (·) y denote the partial derivative with respect to x and y, respectively. After some simple manipulations, (5) and (6) take the following expressions: respectively, where L km , P km , Q km , and R km denote the coefficient of the m derivative of v 0 and φ 0 that multiplies ε k . For instance, L −11 is the coefficient of the first derivative of v 0 that multiplies ε −1 . As a result of the continuity requirements, when ε → 0, the coefficients of k = −1 are equated to zero. The non-identical zero coefficients are the so-called local problems. In what follows, we characterize these problems.

Local Problems
The following local problems have to be considered. The system of local problems (LQ): with interface conditions denotes the jump at the interface. Furthermore, the system of local problems (PR): with interface conditions The LQ system is derived from the coefficients multiplying the derivatives of v 0 in each equation and it relates the local functions N 1 , Ψ 1 , N 2 , and Ψ 2 . Similarly, the PR system is obtained from the coefficients of the derivatives of φ 0 and it relates the local functions Θ 1 , Π 1 , Θ 2 , and Π 2 .

Methodology to Solve Local Problems Using a System of Linear Equations
To obey conditions (13)-(20), the LQ and PR systems need to be solved. More specifically, dy 2 , and d 2 Θ 2 dy 2 need to be obtained. Later, we show that the local functions N 1 , N 2 , Ψ 1 , Ψ 2 , Π 1 , Π 2 , Θ 1 , and Θ 1 are not necessary to find the homogenized problem.
The system of local problems is solved using ideas described in Reference [20]. Here, we summarize the methodology. Let us focus on the LQ system. Integrating Equations (13) and (15) the following system is derived: where λ α ≡ λ α (x), α = 1, 3. Thus, the solutions of the system are To find expressions for λ α in terms of the materials parameters C, e, κ and µ, one needs to use the average operator with respect to the local variable. This average operator is given by where V Y = |Y|. From the periodicity of N 1 and Ψ 1 and taking the average operator in both sides of the Equations (21) and (22), the following system for λ α is obtained: It follows that λ 1 and λ 3 are solutions of the system (24) and (25). Considering C, e, κ, and µ to be constant functions along each material, one gets that (14) and (16) respectively. Finally, solutions of the LQ system are given by (21), (22) and where ∆ = e 2 + Cκ and λ i , i = 1, 2, 3, 4, satisfy the system of equations Similarly, solutions of the PR system (17)- (20) are obtained. Expressions for these solutions are given in Appendix A.

Variational Formulation of the LQ System. Finite Element Method (FEM)
The finite element method is a strong and convenient numerical tool to solve systems of ordinary or partial differential equations. It is an alternative method to approximate the solution for the LQ systems (13)- (16). This is the procedure. Let us take Y = [0, γ) ∪ [γ, 1], where γ ∈ (0, 1). Consider the test functions v i , i = 1, 2, 3, 4, multiplying the LQ problem. Integrating along Y, using integration by parts and taking into account the continuity of the function at the interface, the following variational formulation for the LQ system is derived: Periodic conditions for the local functions and perfect contact at the interface are considered. Now, using Equations (31)-(34) and the standard finite element method [41] on the domain Y, the local functions are approximated.

Effective Coefficients and Homogenized Problem
To find the effective coefficients, the average operator is used on the coefficients of ε 0 in Equations (11) and (12), i.e., one deals with L 0m , P 0m , Q 0m , and R 0m . Taking into account the periodicity of the local functions as well as C, µ, e, and κ, the non-identical zeros are Substituting the expressions of the solutions of the LQ and PR systems, which are given in (21), (22), (28), (29), and the Appendix Equations (A1)-(A4), into the Equations (35)-(42), the following relations are obtained Therefore, the effective coefficients of Equations (11) and (12) are The original problem (5) and (6), by means of the homogenization technique, is transformed into the following homogenized problem (in terms of the effective coefficients): where σ 0 , E 0 , and D 0 are the effective stress, effective electric field, and effective electrical displacement, respectively.

Model Validation with the Particular Case of Piezoelectric Materials
To validate the model, our goal is to reproduce the results of the one-dimensional piezoelectric composite analyzed in Reference [30]. The constituents of this one-dimensional piezoelectric structure are lead zirconate titanate (PZT) C-91 and P-82. The materials properties are given in Table 1 of Reference [30]. For this particular case, the flexoelectric coefficient µ associated with the components of the structure is considered zero. Therefore, the equilibrium problems (5)-(8) take the expression (1) and (2) of Reference [30]. Using the asymptotic expansion (9) and (10) and following the methodology described in the previous section, the non-vanishing local problems are (13), (15), (18), and (20), which corresponds with Equations (23) and (24) in Reference [30]. Finally, the expressions of the effective coefficients for this piezoelectric composite are given by Equations (43), (44), and (45).
In Table 1, the values reported for the effective coefficients C, e, κ of this piezoelectric composite considered in References [19,30] are compared with the results obtained using the present model (Section 3) and the results obtained using FEM (solving the variational formulation of the LQ and PR systems). The volume fraction of C-91, for the numerical calculation, is considered to be V 1 ∈(0, 0.4, 0.8, 1). To avoid singularities of solutions of LQ and PR systems, the value of flexoelectric property for the piezoelectric composites is taken as µ = 10 −20 . Suffice to say that the results reported in Reference [19,30] are reproduced using both the AHM (present model) and FEM (variational formulation) methods. Therefore, the piezoelectric composite results are obtained as a particular case of the flexoelectric composite formulation developed here. Table 1. Comparison of the effective coefficients C, e, κ for a piezoelectric composite computed using the methodology described in Section 3, the finite element method (FEM), and the results reported in References [19,30].

Flexoelectric Composite Rod without Ideal Periodicity
It is interesting to study the the above-described methodology for the case of random composites, i.e., composites without periodicity. A two-elements flexoelectric rod where the constituents are L 1 (BaTiO 3 ) and L 2 (SrTiO 3 ) with a fixed length of 1 unit is considered. Assume Ω N = [0, N], for every 0 ≤ n < N, n ∈ N, the region [n, n + 1] is occupied by the layer L 1 with probability p = 0.7146. The numbers N 1 and N 2 represent the amount of components L 1 and L 2 in Ω N , respectively, and N = N 1 + N 2 . The average operator is considered This average operator is a particular case of the Equation (2.2.1) of Reference [31] for one-dimensional structures. From the results described in Section 3 and considering the average operator (51), the effective properties (43)- (46) are easily computed for a structure without periodicity. To compare the behavior of the effective properties for random composites with the ones for periodic structures, a periodic composite Ω is considered. The constituents of Ω are BaTiO 3 with volume fraction p = 0.7146 and SrTiO 3 with volume fraction 1 − p. In Table 2, the material properties of BaTiO 3 and SrTiO 3 are given. Figure 2 illustrates the behavior of three experiments with random composites considering a probability p = 0.7146 for L 1 and a composite with periodicity and a volume fraction of p = 0.7146 for BaTiO 3 for N ∈ [1, 2 25 ]. The three random composites were randomly generated with MATLAB using the so-called rand function. Figure 2 shows that the values of effective properties for non-periodic composites coincide with the values of the effective coefficients for a periodic structure when the number of layers increase. It follows that the methodology is valid for structures with no periodicity, making it possible to approximate the effective properties for random composites considering the periodic structures. Similar results are obtained in Reference [42], for the particular case of an elastic Fibonacci laminate composite.

Influence of the Flexoelectric Property on the Effective Coefficients
In this section, the influence of the flexoelectric property on the effective coefficients is studied. A bi-material one-dimensional composite is considered. In Table 2   In Figure 3, the influence of the flexoelectricity in combination with an active polymer PVDF which has piezoelectric properties on the effective characteristics is studied using three composites. The first composite is made of the flexoelectric material BaTiO 3 and the piezoelectric active polymer PVDF; the second one is the composition of the flexoelectric composite SrTiO 3 and PVDF; the materials of the third composite are the piezoelectric constituents PZT-7A and PVDF (see Reference [38]). With the values given in Table 2 and the methodology described in Section 3, the effective coefficients for each composite are obtained. Figure 3 illustrates the behavior of the effective coefficients for different volume fractions of the composites. The effective coefficient C continues increasing for the three combinations of materials. In addition, as is shown in the figure, the influence of the flexoelectric materials on the piezoelectric and the dielectric effective coefficients is remarkable. For both BaTiO 3 -PVDF and SrTiO 3 -PVDF, the piezoelectric effective coefficient e increases until V 1 = 0.9 and V 1 = 0.8, respectively, from which the property decays. On the other hand, for PZT-PVDF, the piezoelectric effective tensor decreases for all values of V 1 . Similarly to C, the dielectric effective coefficient κ increases with the volume fraction for the three combinations of materials, although for the composites formed by a flexoelectric material this property is greater than that of the PZT-PVDF composite. The effective flexoelectric coefficient µ varies only for the composites with a flexoelectric component (either BaTiO 3 or SrTiO 3 ).  To study the influence of non-flexoelectric and non-piezoelectric constituents in the effective properties for composites with a flexoelectric component, three structures are studied. One of the structures is the combination of the flexoelectric materials BaTiO 3 and SrTiO 3 ; the constituents for the second structure are SrTiO 3 and a non-piezoelectric polymer; the last structure is made of strontium titanate (SrTiO 3 ) and a piezoelectric material PZT-7A. In Figure 4, the effective properties of these three composites for different volume fractions are shown. An important detail to highlight is that for all the effective properties in Figure 4, the values of the effective coefficients of the SrTiO 3 -polymer are lower than the values of the effective coefficients of the other two combinations. The non-piezoelectric polymer properties diminish the values of the composite effective coefficients, compared with the other composites. On the other hand, the combination of the two flexoelectric materials reinforce the values of C, e, κ, and µ, compared with the composite made of the non-piezoelectric polymer. It is shown in Figure 4 that the effective piezoelectric property e is almost zero for the combination of SrTiO 3 with a non-piezoelectric polymer. A similar situation is shown for the flexoelectric parameter, i.e., e remains close to zero for the composites that involve a non-flexoelectric component (either SrTiO 3 -polymer or SrTiO 3 -PZT). Therefore, it can be concluded that for composites made of flexoelectric material with non-flexoelectric material, the global behavior of the effective coefficients of the structure is similar to a non-flexoelectric material. A parallel situation occurs when considering the composition of a piezoelectric with a non-piezoelectric structure as it is illustrated in Figure 4. Hence, for composites made of a piezoelectric material with a non-piezoelectric material the global behavior of the effective coefficients of the structure is similar to a non-piezoelectric material. The dependence of the effective coefficients C, e, and κ on the local functions dN 1 dy , dΨ 1 dy , dΘ 1 dy , dΠ 1 dy is shown in Equations (43)- (45). From Equations (21), (22), the Appendix Equations (A1), (A3), and the solutions of the systems (30), (A5), the non-influence of the µ parameter in the local functions is derived. Therefore, the effective coefficients C, e, and κ are not affected by the flexoelectric component in the case of non-flexoelectric materials. From the homogenized problem (47)-(50), it can be followed that the influence of the effective flexoelectric component in the effective stress and electric displacement cannot be underestimated.
The investigation of the flexoelectric effect goes beyond the theoretical analysis of the materials properties to include real measurements [43]. Recent papers have established detailed descriptions of their experiments to determine the flexoelectric properties of PVDF. More precisely, References [44] and [45] describe experiments to find the effective flexoelectric components µ 1123 and µ 2312 of the fourth order tensor µ ijkl for polyvinylidene fluoride. On account of the great difficulty in obtaining some experimental measurements, the relationship between the strain gradient and torque is deduced theoretically and further verified with finite element analysis. The approach is applied to test responses in bars machined from bulk polyvinylidene fluoride [45]. On the other hand, in Reference [46], the flexoelectricity of prototypical semicrystalline polymer, α-phase PVDF, films are investigated. The paper presents a step-by-step description of the experiment highlighting the direct flexoelectric effect in the α-phase PVDF films.

Solution of the Heterogeneous and Homogenized Problems
In this section, the homogenized problem (47)-(50) is solved for a flexoelectric one-dimensional structure. In order to compare the approximation between the solution of the heterogeneous problem and the homogenized problem, the finite element method is used. Consider Ω a two-element structure made of barium titanate (BaTiO 3 ) and strontium titanate (SrTiO 3 ) periodically distributed along Ω. In order to illustrate the example, the prescribed conditions for the heterogeneous (5)-(8) and the homogeneous (47)-(50) problems are given as follows: u 0 = 1, w 0 = 1, ϕ 0 = 0, S 1 = 0, r 1 = 0, and τ 1 = 0. Two different volume fractions for BaTiO 3 are considered V 1 = 0.5, case 1, and V 1 = 0.8, case 2. The external forces are taken as f (x) = e −x for case 1 and f (x) = sin(x) for case 2. The body charge density ρ is considered identically zero along Ω. The values of the effective coefficients of the homogenized problem are computed using the methodology described in Section 3 and are shown in Table 3.
In Figure 5, a comparison of the stress obtained from the solution of the heterogeneous problem (5)-(8) and the homogeneous equivalent problem (47)-(50) for ε = 1/5 is shown. As a result of the continuity of the coefficients of the homogeneous problem, the stress function has a smooth behavior along Ω. On the other hand, the discontinuity of the coefficients in the heterogeneous problem leads to the obstacles in the solution, see Reference [47].

Conclusions
In this manuscript, the equilibrium equation for a flexoelectric composite has been studied. The two-scales asymptotic homogenization method is used to find the homogenized problem formulation. The work presents the expressions of the effective coefficients for the case of a flexoelectric structure. The local problems are derived along with their corresponding variational formulations. The manuscript offers a detailed description to solve the local problems for one-dimensional structures. The effective coefficients reported in References [19,30] for piezoelectric composite are computed here as particular cases of both the present model (Section 3) and the finite element method. The methodology is used to find the effective properties of random structures. The results illustrate the convergence of the effective coefficients for random structures to the case of a composite with ideal periodicity. In addition, the influence of the flexoelectric parameter on the global behavior of composites has been studied. The effective properties of three bi-materials (composites) made of combinations of BaTiO 3 , SrTiO 3 , PVDF, a soft non-piezoelectric polymer, and PZT-7A were obtained. The influence of the flexoelectric parameter (µ) on the effective stress and electric displacement was shown. As an example, a comparison between the solutions of the heterogeneous and homogeneous problem for a flexoelectric structure is presented. A good approximation between the solutions can be appreciated when ε is zero. The results presented in this manuscript can be extended to more complex structures with wavy behavior and three-dimensional composites.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Solutions of the (PR) System
The solutions of the (PR) system (17)- (20) are where ∆ = e 2 + Cκ and λ i for i = 5, 6, 7, 8, satisfy the following system of equations