Effect of topological imperfections on the electro‐mechanical properties of structured piezoelectric particulate composites

A novel method based on the Virtual Particle Mori–Tanaka (VPMT) is developed to predict the effective electro‐elastic properties, d33 and g33, of structured piezoelectric particulate composites with improved accuracy by means of a single parameter related to the spatial distribution of imperfectly aligned rod‐like PZT particles. The VPMT method is found to have excellent prediction capabilities for idealized particle configurations. Several new correction functions are presented to capture the drop in piezoelectric composite’s electro‐elastic properties as a function of topological imperfections. These imperfections are related to longitudinal and lateral inter‐particle spacings and the topology of the chain like structures themselves. The functions are evaluated in detail and show physically consistent behaviour.


Introduction
Piezoelectric ceramics distinguish themselves from other materials through their electromechanical coupling which makes them suitable for sensor, actuator, and energy harvesting applications [1][2][3]. The inherent brittleness and difficult manufacturability of these materials limits their use for applications in which mechanical flexibility is required. Currently available flexible piezoelectric materials, such as PVDF, show poor piezoelectric coupling and high susceptibility to piezoelectric degradation due to environmental influences, whereas piezoelectric ceramics are resilient and chemically inert.
To bypass the inherent drawbacks of piezoelectric ceramics, piezo-polymer composite materials (PCMs) have been developed to combine the superior piezoelectric properties of piezoelectric ceramics with the mechanical properties of a polymer matrix allowing for high piezoelectric coupling characteristics with increased flexibility and toughness. These composites feature a particulate piezo phase embedded in a polymer matrix. Newnham et al [4] distinguished 9 different types of connectivity of which the random particle distribution, 0-3 composites with particles mutually unconnected and the matrix fully connected in 3 dimensions, and a structured particle distribution, 1-3 composites with fibrous particles fully aligned and connected in one dimensions, are the most relevant. Perfect 1-3 connectivity offers rather superior piezoelectric properties in the direction of 2 alignment as compared to 0-3 composites. However, production of 1-3 composites is a complex process which often results in composites with a significant amount of microstructural imperfections, making it hard to predict their effective electromechanical properties using existing analytical models [5][6][7][8][9]. As the current models do not incorporate microstructural imperfections, they generally overestimate the effective electromechanical properties. Better predictions may be obtained by using numerical schemes, such as finite element modelling, or the self-consistent and generalized self-consistent schemes Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. [6,[10][11][12][13]. FEM based methods use numerical techniques that take into account electrical and mechanical boundary conditions in a composite system that is represented by a unit cell. An advantage of FEM methods is no restrictions on connectivity patterns and inclusion sizes. However these are time consuming and computationally expensive due to the large amount of iterations involved.
Recently many theoretical micromechanics models of heterogeneous piezoelectric systems, such as porous piezo-ceramics and ceramic-polymer composites, have emerged among which the self-consistent and Mori-Tanaka (MT) methods have received the most attention and use [14][15][16][17][18]. The generalized self-consistent method is based on the three-phase model, i.e. an inclusion and a surrounding ring matrix constitute a representative unit cell, which, in turn, is embedded in an infinite composite [19,20]. MT is based on the two-phase model of inclusion/matrix. It has been shown that the generalized self-consistent method and MT method can verify each other's results [21][22][23].
In this paper a novel computational method is proposed to predict the electromechanical properties of structured piezoelectric ceramic-polymer composites with improved accuracy, taking into account several (quantified) topological imperfections in the arrangement of the piezo electric particles. The method is essentially based on the established MT model reformulated by Benveniste [23] and adopted by Dunn and Taya [24] who expanded the Eshelby tensor to incorporate piezoelectric behaviour [14]. The MT theory considers a single inclusion in an infinite matrix aligned in a specified direction and has been proven to be an accurate and effective model for the homogenization of particulate composites. For this study a numerical approximation for the Eshelby tensor is used which has been extensively described and validated in a FEM analysis by Andrews et al [25].

Homogenization of 0-3 composites
To approximate the full set of electro-elastic properties of a 0-3 composite with significant accuracy, the MT method can be applied to homogenize the effective properties. It is assumed that the particles are aligned in the poling direction and that the inter-particle distance is large. By itself, the MT method cannot be used for piezoelectric composites, which was addressed by Dunn and Taya [24] who derived the following relation for the effective electro-elastic modulus, E c , of a particulate electro-active composite.
where E r , c r , N and A r are the electro-elastic moduli of the rth phase, the volume fraction of the rth phase, the total number of phases (including the matrix material) and the concentration tensor for the rth phase, respectively. The first phase is the matrix phase and it is denoted with subscript 1. Of all available methods that have relatively short computation times, the MT method has proven to give the most accurate results [25][26][27][28]. The MT method uses a concentration tensor which is represented by A r dil is the dilute concentration tensor which is expressed as where S r is the Eshelby or constraint tensor for the rth phase. The Eshelby tensor is a function of the geometry of the ellipsoidal particle or inclusion with semi-axes a, b and c, and the electro-elastic properties, E 1, of the matrix phase.
For an ellipsoidal inclusion the following relations define three distinct particle shapes.

( )
A numerical approximation of the constraint and Eshelby tensors is provided by Andrews et al [25].

Randomly oriented particles
To approximate the effective properties of composites with randomly oriented particles the MT method can be applied in combination with the rule of mixtures. The first step is to calculate the effective electro-elastic properties for composites with aligned slender cylindrical particles, oriented in discretized Euler rotations, which replicate all possible directions. It should be noted that the poling axis remains in the 3-direction for all these composites. The previously discussed Eshelby tensor can be rotated with tensor transformation laws which allows the generation of electro-elastic properties for composites with slender cylindrical particles with different orientations. Because the electro-elastic Eshelby tensor is an extension of the Voigt notation and essentially not a true tensor the transformation laws cannot be applied directly. Therefore, the transformation laws should first be applied to the separate parts of the tensor before it is recompiled in the Voigt extension. The modelling of 0-3 composites with randomly oriented particles using the rule of mixtures is elaborated in more detail by Andrews et al [25].

Homogenization of 1-3 composites
In this section a new method is presented for the prediction of the effective properties of structured piezoelectric composite materials, also referred to as (quasi) 1-3 composites, with the piezoelectric charge constant d 33 and the composite permittivity κ in particular since they can be used to determine the piezoelectric voltage coefficient, g 33 , and the energy efficiency. The theory behind the model will be elucidated and an experimental validation is performed to demonstrate the model's capabilities.
Columnisation of particles and its effect on the properties of PCMs has been described by Bowen et al [29] and Van den Ende et al [5,7,8,30] who showed that columnisation can significantly improve the effective electroelastic properties of PCMs in the poling direction. They proposed models for the prediction of the composite permittivity and the d 33 , respectively. Bowen's model for the permittivity assumes perfectly stacked particles with a discrete inter-particle spacing in the poling direction, also referred to as axial particle misalignment. Van den Ende's model for the d 33 and g 33 expands on this by additionally taking angular misalignment of the cylindrical particles into account. These models give good preliminary prediction results for 1-3 composites with a limited amount of imperfections and act as the upper limit for the newly developed model in which additional topological imperfections will cause a further reduction of the composite's properties. Benveniste [20] concluded that the MT method already takes particle interaction into account up to a certain extent which implies that it can be used for homogenization of composites where particle interaction is present even though an infinite particle spacing is assumed initially.
In the current approach the MT method with Eshelby tensors is used twice in succession to determine the properties of composites with columnised particles by redefining a column as a single virtual particle (VP) which is in essence a composite in its own right [5]. The particulate columns or chains have a high aspect ratio and the proposed method originates from the fact that for these high aspect ratios the values of the composite performance parameters as a function of the particle volume fraction converge towards a common limit [25]. The contour of the limit for high aspect ratios is traversed by varying the volume fraction of particles within the column.
For construction of a VP it is assumed that the particles in the column are close to each other within a section that can be captured by a rectangle and that the formed chains or columns have an aspect ratio of at least 20, but preferably 100 and up because the numerical error is smaller for higher aspect ratios. To have the geometry of the VP comply with the MT method requirements an ellipse is fitted over the rectangle which incidentally also reduces the very high particle volume fraction within the rectangle, increasing calculation accuracy. A schematic representation of a VP is shown in figure 1.
The effective composite properties depend on several nonlinear functions, which are mainly influenced by the volume fraction and aspect ratio of the piezoceramic phase [5,7,[25][26][27][28][29]31]. In order to correct for the high VP volume fraction [26,31] and incorporate columnisation and imperfections through manipulation of the volume fraction of the VP the novel Virtual Particle Mori-Tanaka (VPMT) method is proposed.
The VPMT method is a derivative of a non-dilute MT method proposed by Schjødt-Thomsen and Pyrz [32] where an intermediate step with a comparison medium is used to account for high volume fraction deviations. For structured composites this comparison medium is defined as being a two phase particulate piezoelectric composite with the same initial phases as the basic composite. The volume fraction, c VP , of this composite can be adjusted to account for the high-volume fraction error and nonlinear conversion effects when using the VP with the MT method in a similar way that Schjødt-Thomsen and Pyrz determine the non-dilute solution for the MT method.
where I is the identity tensor, S p is the Eshelby or constraint tensor for the piezoelectric phase, E m is the electroelastic tensor of the matrix and E p is the electro-elastic tensor of the piezoelectric material. A r cm dil , is used in the MT method to calculate the comparison medium's electro-elastic tensor E cm with c VP . In the successive step, the dilute concentration tensor of the VP, A , is used with a the VP's aspect ratio and a corrected composite volume fraction, c VPMT , which is expressed as a function of the piezoelectric material's initial volume fraction in the VP, c 0 VP , and its original volume fraction, c p .
The volume fraction that is used in the second step of the VPMT method is larger than the original volume fraction because there is also matrix material included in the VP. This method distinguishes itself from the nondilute MT method through the change of the particle geometry and volume fraction by means of a VP which accounts for the chain formation. The bounds for c VP are 0<c VP <=c .

VP 0
For c VP =c VP 0 =1 the chains in the composite are uninterrupted and perfectly oriented and the results should coincide with the results of the MT method when the particle geometry is taken to be a continuous wire. By substituting equation (7) in (2), the dilute concentration tensor of the VPMT is obtained as A :

Experimental determination of c VP
The value of c VP can be determined experimentally by comparing test results of 0-3 and 1-3 composites and selecting a value for c VP that captures the differences in performance. The accuracy of the prediction of the electro-elastic properties of the structured composite relies on an accurate determination of the properties of the 0-3 composite.
To demonstrate the prediction capabilities of the newly proposed method is given by fitting the MT method to the experimental results of measured d 33 and permittivity of 0-3 composite samples published elsewhere [5]. Some of the material parameters that are used as input for the MT model may have to be adjusted slightly with respect to the nominal values because the properties of the composites' constituents change due processing of the material. Piezoceramics in particular experience significant changes in piezoelectric and dielectric properties when ground after sintering [33][34][35][36][37] which results in changes of the effective electromechanical properties of the composite in which they are embedded, different particle-matrix interaction [34][35][36], and different processing behaviour [7,38,39,41]. Because the piezoceramic powder is too small to conduct reliable property measurements and because the particle size and shape are not uniform the material properties were altered in such a way that the MT method prediction is accurate. Since the MT method is considered to be very reliable it is expected that the result in this case is similar to the result that would have been calculated had the right material properties been available.  It is shown that when the MT method is applied well, the newly proposed method can yield exceptionally well matching predictions for structured composites, here shown for c VP =0.06. The consistency of the analysis is demonstrated by the optimal c VP fit value is the same for d 33 as g 33 . The fit value of c VP of 0.06 is considerably lower than the value of c VP =0.5, which was approximately the initial volume fraction in the VP where c 0 VP =0.55. Clearly the actual spatial arrangement of the particles deviates significantly from an equi-spaced perfectly aligned arrangement.
The c VP as derived here can serve as a simple yet informative figure of merit to gauge the quality of real dielectrophoretically or otherwise structured PZT-polymer composites.

Theoretical prediction of c VP using imperfection parameters
While the c VP parameter gives an attractive single parameter to cover the effects of all topological imperfections, it is desirable to analyse the effects of specific types of imperfections in more detail. Topological defects can be distinguished at a level of a single strand or chain of quasi-aligned particles and at the level of aggregated chains. These will now be analysed separately.

Single chain imperfections
The three most commonly observed imperfections in isolated particulate chains are angular misalignment, axial misalignment and lateral misalignment. Single chain formation is most commonly found in composites with low volume fraction of the piezoceramic phase.

Angular misalignment
The misalignment of particles with respect to the direction of the applied electric field causes the effective properties of the composite in the field direction to be reduced. This is due to the fact that the particles are poled in the field direction. A misalignment of the particles will cause the effective particle length in the poling direction to be reduced, therefore the particle aspect ratio as experienced by the composite is smaller as well. The effects of the average orientation angle of chains of particles and the average orientation angle of particles within the chains have been described by Van den Ende et al, as seen in figure 3 [5].
The effective AR of particles is expressed by Van The effective aspect ratios can directly be used with the VPMT method, leading to a decrease in effective properties when misorientation is added, as is to be expected ( figure 4). In the figure it can be seen that the calculated results for the VPMT method are lower than those of the MT method for structured composites as long as the misorientation angle is less than 10 0 the conventional MT model describes the properties adequately. For higher misorientation angles the MT method rather over-estimates the properties.

Axial misfit
The average axial misfit, l vs , see figure 5, can directly be used with the VPMT method or alternatively be translated into an efficiency factor for c VP which allows for crossing the AR-limit boundary by reducing c VP to mimic the increase of the inter-particle spacing.
A decrease of the volume fraction due to an increase in particle spacing will lead to exponentially decaying effective properties [5,7,8,25]. The axial misfit, shown in figure 5, is averaged over the inter-particle spacing between vertically adjacent particles. The proposed function to describe this is l cos tan . 11 ips vs h = ( ( )) ( ) In figure 6 it is shown that multiplying c VP by η ips drives the effective properties towards those of the MT model, which is as expected since the MT model assumes an infinite spacing between particles.
The simulations on the effect of axial misfit were made for AR=5. Reducing the aspect ratio decreases the initial values of the VPMT and MT methods for zero spacing and the value for which the VPMT method attains the same value as the MT method and vice versa, which is also found in the literature [5,7,8,25].

Lateral misalignment
Another imperfection that can be defined is the average lateral misalignment, l lat , between the ends of adjacent particles in a chain. This parameter will have a large effect in particular for particles with a small aspect ratio, the reason being that the force may not travel efficiently through the particle but can partly bypass it through the matrix. For large particles however, the influence will be less since the force can travel from one particle to another by means of shear and a higher matrix loading in the vicinity of the particle edges, causing a larger portion of the strain energy to travel through the particle instead of the adjacent matrix. The average lateral misalignment is defined as the average of lateral distances between the ends of adjacent particles within a chain ( figure 7). It is assumed that the adjacent particles are part of the same chain. To be part of a chain at least one end of a particle should have a lateral misalignment that is smaller than the particle width.
The average lateral misalignment is incorporated in the VPMT model by using it to define an efficiency factor that gauges how longitudinal stress is transferred between particles. For this case the stiffnesses levels of the composites' constituents are deemed unimportant since the force through the thickness of the composite will have an affinity for the stiffest path and these effects are taken into account intrinsically by the MT method.
It can be imagined that if the particles are perfectly stacked the load path will travel uninterrupted through the chain. If edge misalignment is introduced, it is logical that an increased amount of energy will be lost to the matrix in order to transfer the load from one particle to another. The amount of energy that is lost will be a  function of the particle aspect ratio. The efficiency function describes the change in effective properties, such that its impact is low for high aspect ratios and vice versa. It also includes factors that describe the drop in the performance parameters due to increased lateral spacing and hindered stress propagation. The proposed expression that complies with these demands is c VP is multiplied with the efficiency factor to obtain a new c VP value. Even though the actual behaviour is highly nonlinear the approximation was found to be relatively accurate for lateral spacings that are smaller than twice the particle thickness. Figure 8 shows the decay of the d 33 and dielectric constant for different particle aspect ratios. As shown the effect of the lateral displacement is significant and more so for shorter fibres.

Aggregated chain imperfections
The imperfections that occur in single chains also occur in situations where the single chains in a structured composite have agglomerated into multi-chains or columns, a behaviour that is found commonly for high volume fractions. Additionally, such aggregated chains show two additional topological imperfections: particle overlap and the bridging of gaps by particles. Even though these imperfections may occur in single chains as well, for dielectophoretically structured composites such forms of topological imperfections are unlikely given the physical force field causing the alignment of the particles.

Overlap
Particles that overlap cause the surrounding stress-field to be augmented and therefore the strain in each involved particle changes. Because the connection between the particles is through matrix material more energy will be absorbed by the matrix resulting in less strain in the particles and lower effective composite properties. The average overlap length, l ol , is defined as the average fraction of the mean particle length that overlaps with any adjacent particles ( figure 9).
In addition to the average overlap length, the efficiency function for the overlap should contain a factor that describes the energy transfer as a function of the material stiffnesses. If the matrix material is considerably less stiff than the piezoelectric material, the piezoelectric material will attract a large portion of stress. However, if the materials have approximately the same stiffness the force will flow straight from the top of the composite to the bottom without alterations in the load path and without being diverted to the piezoelectric particles. The force transfer between particles in a matrix with large differences in stiffness and dependency on particle shape and spacing is highly nonlinear. In order to incorporate the energy loss in an efficiency function, the expected physical behaviour should be broken down into comprehensible steps.
Consider an uninterrupted (ideal) chain with locally a particle that has twice the thickness of other particles and a matrix material with negligible stiffness. In the thicker particle the strain will be half that of the rest of the chain. The total output of this section remains the same since the strain is halved, but the area is doubled. If the particle is split into two equal parts across its length and gradually moved away from the chain, it can be imagined that the strain in the particle in the chain will increase and that the strain in the separated particle will decrease. Because the matrix deforms easily the force will be distributed equally over the parts of the split particle for zero lateral spacing and it will be only on the in-chain particle for a large lateral spacing. This is due to the matrix deforming freely causing the strain to be dissipated before it reaches the separated particle. Because the longitudinal strain in a piezoelectric particle is commonly small with respect to the particle's dimensions, it is expected that the force that is transferred to the separated particle diminishes quickly for increased lateral  distance. When the matrix material is stiffer, it will absorb more energy upon deformation, but its ability to transfer force between laterally spaced particles becomes more prevalent.
The energy that is lost through overlap is proposed to be captured with the following efficiency factor l l where l olat is the average lateral distance between the surfaces of the overlapping areas, where Y m and Y p is the stiffness of the matrix and particle respectively. The three parameters that influence the overlap efficiency function show different kinds of behaviour as shown in figures 10-12. Contrary to other imperfection parameters the particle overlap does not converge toward the MT results since it is assumed that there still are columns present in the material. Since the proposed efficiency function does not converge towards the MT model, it is more difficult to fine tune the function and to assess if the impact has the right magnitude. It can however be concluded that the function is very versatile and therefore it offers good prospects for fine tuning capabilities.

Gap-bridging
The average gap-bridging length, l b , is defined as the average fraction of the mean particle aspect ratio that is adjacent to a vertical gap between two stacked particles and overlapping both particles' ends ( figure 13).
The gap-bridging effect is a positive effect which increases the effective properties of the VP. The spacing between the bridged particles is considered to be axial misalignment. The gap-bridging function converges to its neutral point (which is 1) for a large lateral spacing and high matrix stiffness. The efficiency function for the gap-  bridging parameter contains the average fraction of the mean particle length that is bridging a gap, l b , the average lateral distance of the particle from the gap, l olat , the average number of bridges per particle, n b , and the stiffnesses of the composite constituents, Y m and Y p , for reasons mentioned in the previous section on particle overlap. It is proposed to use an efficiency function similar to the function that was used for the particle overlap. The gapbridging factor η b that is proposed is Since the formula is quite similar to that for particle overlap only the effects of the number of bridges per particle and the bridged length are shown in figure 14 to illustrate that it is a beneficial function and that an optimum can be found for certain combinations of parameters in which the energy is distributed more efficiently over the piezoelectric phase.
The impact of gap-bridging has not been researched thoroughly for particulate composites and as for the other imperfections more research is needed to check the validity of the proposed efficiency function. It can however be stated that the trends are as desired. The gap-bridging efficiency function showed trends similar to that of the overlap function and therefore it has the same difficulties with assessment of the magnitude of the impact and fine-tuning capabilities. However, the gap-bridging function is very versatile and offer good prospects for fine tuning when more insight is gained in the mechanics behind this imperfection.

Overall Imperfection parameter
The efficiency parameters defined in the previous sections are all multiplication factors to adjust the c VP value and therefore they can be combined using the following equation for a correction factor that takes into account multiple topological imperfections.   The imperfection volume fraction is used in the calculation as input for the VPMT model instead of c VP and shows a lot of potential to obtain accurate property predictions. The range of η is theoretically 0<η<=1. For empirically optimized dielectrophoretically structured composites the imperfection factor is about 0.8 [5,9,16] showing that some significant improvement in properties can be obtained.

Conclusions
A new VPMT method for the prediction of the properties of particulate composites with imperfections was proposed. The model has the MT model as a lower bound and it coincides with the MT method for high aspect ratios, complying with the desired behaviour which ranges between the dilute limit and the saturated limit.
Efficiency functions were proposed for angular, axial, and lateral misalignment and particle overlap. Angular misalignment is taken into account in a very direct manner and it has the largest impact on the effective properties, yielding a 10% decrease in performance parameters for an 8 degree misalignment angle. Axial misfit causes the performance of the composite to converge towards the results of the dilute MT model rapidly depending mainly on the particle aspect ratio. The other imperfection parameters' efficiency functions have a moderate influence but together they can reduce the effective composite properties by a further to 20%. The only beneficial topological defect leading to an increase in properties over a perfectly structured composite is particle gap-bridging which shows an optimum value for which the strain energy is most efficiently distributed over the surrounding piezoelectric particles. However, the effects of gap-bridging are low compared to the influence of other, clearly detrimental, parameters. Taking realistic degrees of topological imperfections in quasi 1-3 composites it is predicted that the effective dielectric and piezoelectric properties of currently produced quasi 1-3 PZT polymer composites are 20%-50% of the values for perfectly aligned and spaced particulate composites.