Generalized locally-exact homogenization theory for evaluation of electric conductivity and resistance of multiphase materials

The locally-exact homogenization theory is further extended to investigate the homogenized and localized electric behavior of unidirectional composite and porous materials. Distinct from the classical and numerical micromechanics models, the present technique is advantageous by developing exact analytical solutions of repeating unit cells (RUC) with hexagonal and rhomboid geometries that satisfy the internal governing equations and fiber/matrix interfacial continuities in a point-wise manner. A balanced variational principle is proposed to impose the periodic boundary conditions on mirror faces of an RUC, ensuring rapid convergence of homogenized and localized responses. The present simulations are validated against the generalized Eshelby solution with electric capability and the finite-volume direct averaging micromechanics, where excellent agreements are obtained. Several micromechanical parameters are then tested of their effects on the responses of composites, such as the fiber/matrix ratio and RUC geometry. The efficiency of the theory is also proved and only a few seconds are required to generate a full set of properties and concomitant local electric fields in an uncompiled MATLAB environment. Finally, the related programs may be encapsulated with an input/output (I/O) interface such that even nonprofessionals can execute the programs without learning the mathematical details.


Introduction
The past fifty years have witnessed tremendous growth in the development and use of composite materials, the applications of which extend from microelectronic devices to advanced aerospace structural and engine components [1][2][3]. Compared to the monolithic metals, polymers and ceramics, whose properties are usually stable and do not vary spatially, fibrous composites possess some remarkable properties that can be manufactured and designed for specific applications by tailoring the internal microstructure such as fiber shapes, sizes and array patterns, as well as using reinforcement with different properties [3][4][5][6][7]. For example, monolithic piezoelectric materials have a limited range of coupled properties and pronounced directionality. Those problems may be circumvented by using the piezoelectric materials in the form of composites where optimum electromechanical coupling effects are obtained by selecting proper fiber volume fractions and poling orientations [8,9].
Electric conductivity and resistance are very important properties to be considered in the development of smart devices such as sensors and transducers. The potential benefits that may be obtained from the composite technologies have led to an increasing interest in understanding the mechanical, thermal and electric properties of multiphase and multifunctional composite materials. However, most of the papers in the literature focus on the mechanical properties and thermal conductivity [10,11], and the electric conductivity phenomenon of piezoresistive composites is mostly limited in the experimental measurement [12][13][14] or manufacturing techniques [15]. Tchmutin et al. [16] synthesized the conductive composites with undoped conjugated polymers and conductive fillers and investigated their electophysical and electrochemical properties. Yakovenko et al. [17] reports the morphological effect of carbon nanofillers on the electric properties of polymer based composties. Park et al. [18] tested the effect of including micro-scale secondary fillers on the electrical conducitivity of hybrid nanotube composites. Feng et al. [19] prepared Cu-coated SiC whiskers into coppermatrix composties and tested the effect of whisker orientation on the electrical conductivity.
Although numerous experiments are conducted in the effort of improving the electrical performance of complex material systems, only a few theoretical models are proposed to understand the micro-/nano-mechanics behavior of composites. Some early contributions that study the microstructural effects of piezoresistive composites adopted the classical models, such as the Hashin-Shtrikman bounds [20] and Mori-Tanaka model [21], which, however, are mainly based on the mean-field assumption and not capable of considering the fiber-fiber interactions and recovering the localized field distributions. The numerical simulations, mostly finite-element and finite-volume based, can effectively overcome these shortcomings [22][23][24][25][26]. Nevertheless, the detailed mesh discretization and pre-or postprocessing always requires large-scale computational effort and prohibits the numerical approaches' application in multiscale and optimization analysis. Thus, computationally efficient homogenization and localization of the complex material systems still require a sophisticated mathematical micromechanics model.
Despite the widespread use of numerical techniques in the homogenization of composites due to their ability to address complex microstructures and inelastic behavior of constituent phases, interest in elasticity-based approaches has also gain extensive attention within the past decade in light of advances in computational technology, as well as due to the potential advantages offered by these algorithms. The locally-exactly homogenization theory (LEHT) developed by Drago and Pindera [27] is an ideal elasticity-based method for solving micromechanical boundary value problem of unidirectional multiphase materials. The advantages of the LEHT algorithms relative to the numerical approaches are three-fold. Firstly, in contrast with the finite-element and finite-volume numerical approaches that are based on the approximation of the governing equations in the limit when sufficient mesh discretizations are employed, the LEHT satisfies exactly the governing differential equations in a point-wise manner. Secondly, the construction of an input data file in the LEHT is at least an order-of-magnitude faster relative to numerical approaches based on geometric discretization, hence parametric or optimization studies may be carried out with much fewer difficulties. Thirdly, fiber/matrix con-tinuity conditions are exactly satisfied in the interior of the repeating unit cell in the case of LEHT approach, eliminating detailed geometric discretization in the vicinity of the interface to ensure the continuities of electric potentials and electric current densities when variational techniques are employed. As a result, the converged results for the homogenized properties and localized stress field distributions are obtained with much greater efficiency and better accuracy.
The successful applications of the locally-exact homogenization theory with the aforementioned features motivate us to extend the LEHT to investigate the homogenized electric conductivity and local electric field distribution mechanism of a unidirectional composite. The construction of the new LEHT is carried out in three steps. Based on the Trefftz concepts, the electric potential is firstly represented using the Fourier series expansions, which directly satisfy the electric governing equation. The remaining unknown coefficients are then solved by imposing the continuity conditions at the fiber/matrix interface as well as the newly proposed boundary variational principle. Finally, the homogenized constitutive equations are established to obtain the effective conductivity parameters of composite materials. It should be pointed out that relative to the numerical volume integration of current densities of composite constituents adopted by most of the numerical techniques, we present the analytical expressions of homogenized conductivities through introducing the newly defined electric-field concentration matrix. In order to test the accuracy and efficiency of the LEHT, the present technique is first validated against the exact analytical solution in the case when the fiber volume fraction goes to a small value and the solution is thus reduced to the generalized Eshelby problem under far-field loading condition. In the finite-volume fraction case, the homogenized moduli and local electric fields of unidirectional composites with rhomboid and hexagonal arrangements are compared extensively with those generated by the finite-volume direct averaging micromechanics (FVDAM) which has been verified to be of comparable accuracy but with greater efficiency.
The remainder of the present work is organized as follows: Section 2 develops the theoretical framework of the generalized locally-exact micromechanics with electric capability and periodicity implementation, as well as the homogenization relations. Section 3 derives the generalized Eshelby problem and validates against the develop LEHT technique. Parametric studies are conducted to investigate the effects of microstructural parameters and properties of constituent phases on the homogenized properties or localized stress/electric concentrations within periodic 2 Generalized locally-exact homogenization

Unit cell overview
A periodic material with unidirectional inhomogeneities (fibers/porosities) along the x 1 axis is considered in this work. Each inhomogeneity is defined within a repeating unit cell (RUC), which is loaded by effective electric-field components E i (i = 1, 2, 3) and satisfied with periodic boundary conditions. In order to generalize the fiber arrangement in the heterogeneous media, several different geometries of the RUCs are considered, including parallelogram, square and hexagonal arrays, Figure 1. For easy mathematical implementation of the boundary conditions, the electric potential and field are decomposed into averaged and fluctuating components, e.g. and where y j (j = 2, 3) defines the local coordinate system for the microstructural RUCs. According to the Ohm-Kirchhoff law, the current density components of a piezoresistive media are defined as: where i, j = 1, 2, 3. σ ij are the electric conductivity components.
Herein we define the electric current that is the surface integral of the current densities of Eq. (3) along the surface S i : The elastic solution of the present boundary value problem is determined by firstly generating the exact analytical solution that satisfies the interior governing equations within fiber and matrix phases and then imposing the periodic boundary conditions on the RUC's exterior surfaces. The complete solution involves the transformation between the Cartesian and cylindrical coordinates and is hence defined as the inseparable solutions [28]. To solve the exterior problem, a new variational principle is established for the piezoresistive effect of periodic composites [27]: where ϕ = ϕ 0 and I = I 0 are the periodic electric potential and current density constraints imposed on S ϕ and S I , respectively. Taking the first variation of Eq. (5) and considering the periodicities across the opposite surfaces of RUCs yield For a parallelogram microstructure, the explicit periodicities of electric potential and current density take the following form: where d 2 and d 3 are the length and height of the parallelogram array, Figure 1a. It should be noted that Eq. (7) is degenerated to the square array once φ = π/2. Similarly, in the case of hexagonal array, we have where d i (S j ) are the i-th direction projections of the shortest vectors across two opposite faces (j-th surface with its opposite surface) of the unit cell.

The interior problem
In order to establish the continuity conditions at the perimeter of fiber domain, the boundary value problem is initially solved within the cylindrical coordinate system [29]. Thus, we start with the relationship between electrical potential and fields for each constituent of the composites: as well as the relationship between the electrical field and current densities (Ohm-Kirchhoff law): where the superscript i = f , m, and f and m denote the fiber and matrix phases, respectively. σz, σr and denotes the axial, radial and tangential inplane electrical conductivity components, respectively. For a transversely isotropic media, we have σr = σ θ = σ. The governing electrical equations (Maxwell equation) of each constituent can be expressed as: Substitution of Eq. (10) into Eq. (11) yields the Maxwell equation in the cylindrical coordinate in the form: The solution of Eq. (12) can be obtained by assuming: Substitution of Eq. (13) into Eq. (12) yields: where both terms within the bracket are zeros: Solving the Euler equations, Eq. (15), yields the final solution leading to the analytical expressions of current densities: It should be pointed out that the fiber coefficients H f n3 and H f n4 are zeros since the electric potential needs to be bounded in the fiber domain. From the continuity conditions at the perimeter of fiber: the relationship of unknown coefficients between fiber and matrix phases are obtained as where Besides the present internal series expansion functions, we can also employ the traveling wave variable substitution technique to solve the governing partial differential equation. The basic idea is to replace the coupled real variables by an equivalent independent complex variable with a travelling wave variable, helping obtain the solution of simplified ordinary differential equation directly. Direct comparison between the two techniques will be introduced in the future work.

The exterior problem
The remaining unknown coefficients H f n can be determined by implementing the proposed variational principle in Eqs. (5)(6), which imposes the surface potential and current densities on the opposite faces of the unit cell with periodic boundary conditions [27,29]. The two-scale electric-potential representation in Eq. (1) can be reduced to the fluctuating constraints. For instance, we have for parallelogram (square) array and for hexagonal array. Implementation of those conditions into the variational principle, Eq. (6), yields for parallelogram array and for hexagonal array. Both of the variational principles yield the final solution for the unknown fiber coefficients: and H f is a vector comprised of the unknown fiber coefficients. In addition, A and B are the matrices in which the elements are functions of the geometrical and conductivity constants of the microstructures.

Homogenization equations
Distinct from the finite-element-based techniques that usually adopt the volume integration of the fiber and matrix domains in the homogenization, herein we derive the explicit analytical expression for the homogenization constitutive equations. Firstly, we define a generalized Hill electric-field concentration matrix A f [30] to relate the averaged fiber electric field and macroscopic electric field: The averaged fiber electric field can be obtained in a closed form by integrating its local expressions in Eq. (9) over the fiber domain, leading to the explicit expression of Eq. (25) by imposing one unit macroscopic electric field component at a time, with the other components kept zeros, and we have: Finally, the homogenized constitutive equation for the two-phase unidirectional piezoresistive composites is where σ * is the homogenized piezoresistive matrix where the effective piezoresistive conductivity constants are given as: 3 Numerical results

Convergence study
The success of the present technique depends on the fast convergence of the series expansion of the internal solutions as well as the periodic boundary implementation. Thus, we first demonstrate the robustness of the developed  Figure 2, which are normalized by the corresponding converged results when 15 harmonic terms are employed. Two cases with relatively large material property contrasts are considered, namely, σ f ⧸︁ σ m = 100 and 0.001, representing conductive and nonconductive fibers embedded in conductive matrices, respectively. Unit cells with hexagonal, square and parallelogram arrangements are used in the calculation. The angle between the two connected edges of the parallelogram unit cell is taken as 75 ∘ . The ability to generate the homogenized properties and local electric field distributions of/within an arbitrary parallelogram unit cell is a major advantage of the present elasticitybased approach. It should be noted that the discretization of arbitrary parallelogram microstructures cannot be easily achieved based on the finite-element and finitevolume-based numerical techniques that require complicated mesh discretization techniques to accomplish. As observed, the normalized effective properties converge to "1" with around 8 harmonic terms for both cases. No difference is shown further increasing harmonic terms in the electric potential representation, Eq. (16). This indicates that the convergence rate of the locally-exact homogeniza-tion theory with the Fourier series expansion is very rapid and the effective electric conductivities can be calculated accurately with only a few harmonics.
We proceed to illustrate the convergence of localized current density distributions within hexagonal and parallelogram composite microstructures during a unit electric loading E 2 = 1V/m as it is the most demanding. The material system investigated herein is a conductive fiber 1000 S/m embedded in a less conducive matrix 200 S/m. The volume fractions for both unit cells are prescribed as 20%. The angle between two connected edges of the parallelogram unit cell is taken as 75 ∘ . The small fiber diameter relative to the overall unit cell dimensions and relatively large property contrast result in high electric field gradients at the fiber/matrix interface, providing a demanding test of the methods' accuracy. Figure 3 illustrates the transverse local current density distributions J 2 for three different harmonics, N = 2, 5, 10, respectively. For 2 harmonics and within both microstructures, the basic features of the local current density distributions are reasonably well-captured. The periodicity feature of the parallelogram unit cell in y 3 direction is, however, poorly predicted. 5 harmonics yield improved results for the current density distribution that differs from that obtained based on 2 harmonics. Further increasing harmonics to 10 produces converged local current density field but no fundamental differences are observed between 5 and 10 harmonics, as also observed in Figure 2. Thus, the convergence of the LEHT is very fast and doesn't require a large number of harmonic terms.

Generalized Eshelby problem
Next, the locally-exact homogenization theory with electric capability is validated against the generalized Eshelby solution that describes an inclusion embedded into an infinite matrix, which is subjected to far-field current density loading J ∞ 2 ≠ 0. According to the transformation equations between the Cartesian and cylindrical coordinates, the boundary conditions can be re-expressed as: Similar to the analytical solution of Eq. (17), the expression for a transversely isotropic piezoresistive material is given as: while the rest unknown coefficients are obtained by implementing continuity conditions between fiber and matrix phases: which generates systems of equations: from which all of the remaining unknown coefficients are solved.
To simulate the far-field loading condition in the modified generalized Eshelby problem, a square repeating unit cell with a dilute volume fraction, namely, 5% in the present case, is employed in the generalized locally-exact simulation. A macroscopic current density is then J 2 = 1 A/m 2 applied to the unit cell. Figure 4 presents a comparison of the inplane local current density distributions, J 2 (y2, y 3) and J 3 (y2, y 3) , in the Cartesian coordinate system near the fiber for a composite with a conductive fiber, based on the exact solution of an infinite plate under farfield loading J ∞ 2 = 1 A/m 2 and the present LEHT technique. Figure 5 presents similar results in the case of a nonconductive fiber. As observed, the results generated by the locallyexact homogenization theory correlates well with the modified Eshelby solutions almost everywhere in the domain of interest. It should be pointed out that a current-density concentration factor (CDCF) of two is observed in at the top and bottom vicinities of the nonconductive fiber, meaning the current density could be magnified two times under far-field loading condition for a conductive matrix with a nonconductive fiber or porosity. In particular, the transverse local current density within the fiber is uniform while the local current density J 3 (y2, y 3) is zero in the fiber domain. The small departures between the locally-exact and Eshelby solutions that are seen occur in the vicinity of the boundary of the unit cell, which is due to the differences in the applied boundary conditions. While far-field conditions are imposed in the infinity of the Eshelby problem, the periodical boundary conditions are applied to the mirror faces of the unit cell in the LEHT theory. Nonetheless, these differences are absent when a smaller fiber volume fraction is employed in the LEHT solution.

Numerical validations and new results
To assess the accuracy of the locally-exact homogenization theory in calculating the homogenized and localized behavior at moderate and higher volume fractions, the homogenized version of the FVDAM theory is extended to fulfill this demand. The FVDAM is among a few micromechanics theories that yield accurate results for both the effective and localized behavior, the accuracy of which is proved to be comparable to the Q8 finite-element method but with greater efficiency and stability. Herein, we use the FVDAM as a gold standard for comparison with our newly extended locally-exact theory. Figure 6 firstly compares the present simulations agains the experimental measurements for the electric conductivities in the longitudinal and transverse directions of a porous nickel. The conductivity parameter of a homogenenous nickel is σ M = 1.41 × 10 7 Ω −1 m −1 [21]. It is seen that most of the measurements are conducted for nickel with porosity volume fractions between 0.2 and 0.4. Both hexagonal and square microstructures are employed for the predictions where the geometric effect doesn't play an important role when the porosity volume fraction is lower than 0.5. However, good agreement is still obtained against the experimental data, indicating that the present simulation tool is reliable in predicting effective conductivities with high volume fractions that may not be easily conducted in the laboratories. Figure 7 illustrates comparison of the complete set of effective electric conductivities of composites with three different matrix properties over a wide range of volume fractions. The conductive fibers' arrangement is modeled as either square or hexagonal arrays. The specific values for the investigated material systems are listed in Table 1, where the fiber electric conductivities are fixed as 1000 S/m with the matrix properties varied from 10 S/m to 200 S/m. The correlation is seen to be excellent for all the material systems and microstructural geometries using 10 harmonics, whose choice was motivated by the results presented in Figure 2. It should be noted that only 5 seconds and 7 seconds are required to generate a full set of effective properties for composites with fourteen different volume fractions under an uncompiled MATLAB environment on PC Intel(R) Core(TM) i5-3320M CPU @ 2.60GHz.  Figure 8 illustrates comparison of the complete set of homogenized electric conductivities in the case of nonconductive fibers or pores in three different matrices in square and hexagonal arrays over a wide range of volume fractions, simulating porous materials. In the locally-exact micromechanics simulation, the electric conductivities of the matrices are the same as those in the proceeding case while the electric conductivities of the fiber phase are prescribed as one-tenth of the corresponding values of the matrix phases. In the FVDAM simulation, the fiber phase is excluded from the composite such that a real porosity is modeled. Once again, the locally-exact homogenization exhibits remarkable correlation with the finite-volume predictions in the low and intermmediate volume fractions. A cursory examination at higher volume volumes indicates that the locally-exact homogenization generates slightly larger effective conductivities. The small differences are due to the way in which the porosity is treated. The contri-bution from the fiber phase is negligible at low and intermmediate volume fractions but becomes noticeable when the fiber-volume fraction is sufficiently large to contribute some load-bearing capability despite the small conductivities. The locally-exact results are therefore consistent with this observation, illustrating that the method is sufficiently sensitive to correctly capture these small effects.
The localized current density distributions, J 2 (y2, y 3) and J 3 (y2, y 3) , generated by the locally-exact homogenization and FVDAM simulation of a hexagonal unit cell with 20% volume fraction during transverse electric loading  Figure 9. The electric conductivities of the fiber and matrix phases are 1000,200 S/m, respectively. As expected, no visible difference is observed between the two approach's predictions. Similarly, Figure 10 illustrates comparison of local electric field distributions, E 2 (y2, y 3) and E 3 (y2, y 3) , within a square unit cell inserted with a nonconductive fiber in the case of LEHT and a porosity in the case of FVDAM. Given the negligible contributions from the nonconductive fiber to electric field in the vicinity of fiber/matrix interface, the LEHT generates generally good correlation with the FVDAM in the matrix domain.
To quantitatively highlight the accuracy of the locallyexact theory relative to the FVDAM results, Table 2 compares effective electric conductivity for rhombic unit cells with two different characteristic angles, φ = 50 ∘ and 70 ∘ , and three different fiber volume fractions, 0.1, 0.3 and 0.5, thus covering a wide range of reinforcement content. Similar research is conducted in Table 3 for the porous materials when the fiber phase is excluded for the same microstructure. Table 4 presents comparison of the effec-    tive electric conductivity for rhombic unit cells with two large property contrasts, σ f / σm = 10,100, and three different volume fractions, 0.3, 0.5 and 0.7 at a fixed angle φ = arccos(1/4). The electric conductivities predicted by the locally-exact theory with N = 10 converge to the FVDAM results to within approximately less than 0.1% and remain stable thereafter. Despite a few publications appeared in the literature that address the homogenized behavior of a rhombic unit cell, substantially less work was done concerning the localized field distributions within a rhombic microstructure [31,32]. Figure 11 presents local electric current density distributions, J 2 (y2, y 3) and J 3 (y2, y 3) , during a transverse electric loading E 2 = 1 V/m for a rhombic unit cell with a connecting angle of φ = 50 ∘ and 30% volume fraction, generated by the locally-exact homogenization theory. Figure 12 illustrates electric field distributions, E 2 (y2, y 3) and E 3 (y2, y 3) , when the connecting angle increases to φ = 70 ∘ and volume fraction increases to 60%. Not surprisingly, the two approaches predict nearly identical results almost everywhere in the interest of domain. We note that the above two cases are very demanding for the finite-volume and finite-element approaches because of extensive mesh refinements are required. However, no meshes are involved in the present simulations.

Summary and conclusion
The electrical conductivity and resistance of the unidirectional composites are investigated in this work. To improve the computational accuracy and efficiency, the Trefftz concept is adopted, where the explicit analytical solutions are developed for the internal fields and an electrical variational principle is proposed to impose the periodic boundary conditions. By avoiding mesh discretiza-tion in the fiber/matrix domain, the execution of the programs is significantly improved: only a few seconds are required to generate a full set of effective properties. Distinct from the effective response generated in most of the papers in the literature, the present work also emphasizes the recovery of local field distributions, which are normally disturbed by the existence of inhomogeneities. The maximum field concentrations are important in indicating crack initiations. The accuracy of the LEHT in generating the micromechanical responses of unidirectional compos- ites is validated against the independently developed Eshelby problem and FVDAM, guaranteeing its robustness and thus deserving to be a numerical standard.
A few other conclusions are presented below: 1. The effective electrical conductivities can be reinforced with the fiber reinforcement or weakened by the porosity existence. More importantly, the effective conductivity in the fibrous direction varies linearly against the fiber/porosity volume fraction while a nonlinear behavior is observed for the transverse conducitivity component. The geometrical fiber/porosity fiber arrangement has a limited influence on the effective responses, especially for the composites with smaller fiber volume fractions. 2. Besides the effective conductivities, the localized electric/current-density fields are also recovered. Distinct from the porous materials, the fiber reinforcement generates higher electric magnitudes within the fibrous domains. However, the field concentrations still occur at the fiber/matrix interfaces that could initiate possible damages within microstructures that are invisible to naked eyes.
In order to facilitate the present theory's application by professionals or non-professionals alike, the related programs are encapsulated into a "black box" with an input/output (I/O) interface, through which the users only need to type in the input data construction without learning the detailed mathematical derivations in Section 2, Figure 13. The effective or localized results are then automatically generated based on the users' demands.