Transverse Thermoelectricity in Fibrous Composite Materials

Transverse thermoelectric elements have the potential to decouple the electric current and the heat flow, which could lead to new designs of thermoelectric devices. While many theoretical and experimental studies of transverse thermoelectricity have focused on layered structures, this work examines composite materials with aligned fibrous inclusions. A simplified mathematical model was derived based on the Kirchhoff Circuit Laws (KCL), which were used to calculate the equivalent transport properties of the composite structures. These equivalent properties, including Seebeck coefficient, electrical conductivity, and thermal conductivity, compared well with finite element analysis (FEA) results. Peltier cooling performance was also examined using FEA, which exhibited good agreement to KCL model predictions. In addition, a survey was conducted on selected combinations of thermoelectric materials and metals to rank their transverse thermoelectricity with respect to the dimensionless figure of merit.


Introduction
Thermoelectric (TE) devices, which can achieve direct conversion between thermal and electric energies, are usually constructed from n-type and p-type semiconducting elements that are operated under the conventional longitudinal mode, where the electric current and the heat flux flow parallel to each other [1,2].In contrast, TE devices may also be constructed from elements with transverse thermoelectricity, i.e., electric current can induce a heat flow in the perpendicular direction or vice versa [3,4].Therefore, transverse TE can potentially decouple the electric current and heat flux, and thus enable new device designs such as thin-film coolers [5,6] and cascading transverse TE devices [7,8].
While most known thermoelectric materials have insignificant transverse Seebeck coefficients (the off-diagonal terms in the Seebeck coefficient matrix) [9][10][11][12], appreciable transverse TE effect could be achieved in composite structures with engineered microstructural anisotropy, such as lamella structures or a matrix material embedded with aligned fibers [1,3].In the former case, a lamella structure can be constructed from alternating layers of two materials with different TE properties, e.g., a semiconductor and a metal, and may exhibit transverse Seebeck effect in the directions neither parallel nor perpendicular to the layer planes [4].
A framework describing the anisotropic transverse TE properties of the lamella configuration was initially developed by Babin et al. [3], and was further discussed by many other researchers [13][14][15].This framework first utilized the Kirchhoff Circuit Laws (KCL) to calculate the equivalent properties in the directions parallel and perpendicular to the layer plane, and then employed a tensor transformation to obtain the properties in an arbitrary direction.The validity of this simplified mathematical model has been demonstrated by multiple studies using either finite element simulations [16,17] or experimental

Mathematical Modeling
In order to examine the effective properties of a composite material with aligned fibers, a representative volume (unit cell) of the composite structure is studied (Figure 1).It is assumed that the composite structure is infinitely large and the fibers are continuous and uniformly distributed in the matrix.Therefore, the effects of external boundaries and fiber aspect ratio can be neglected.For the sake of simplicity, fibers are assumed to possess identical square cross-sections.The unit cell is further divided into two subunits: C 1 consisting of the fiber phase (F) and part of the matrix (M 1 ) and C 2 consisting of the rest of the matrix (M 2 ).
In the current study, we developed a simple KCL-based mathematical model to describe the transverse TE behavior of fibrous composite materials using an approach similar to that developed by Babin and others for lamella structures [13][14][15].This simple KCL model was based on fibers with square cross-sections.The validity of this model was then confirmed by comparison to finite element analysis (FEA) results, which also indicated little difference between fibers with square and circular cross-sections.FEA was further conducted to show that appreciable cooling might be obtained at room temperature.Finally, we surveyed combinations of some thermoelectric materials and metals, and ranked their performance with respect to their figure of merit.

Mathematical Modeling
In order to examine the effective properties of a composite material with aligned fibers, a representative volume (unit cell) of the composite structure is studied (Figure 1).It is assumed that the composite structure is infinitely large and the fibers are continuous and uniformly distributed in the matrix.Therefore, the effects of external boundaries and fiber aspect ratio can be neglected.For the sake of simplicity, fibers are assumed to possess identical square cross-sections.The unit cell is further divided into two subunits: C1 consisting of the fiber phase (F) and part of the matrix (M1) and C2 consisting of the rest of the matrix (M2).Using the principles of the Kirchhoff Circuit Laws (KCL) and following the steps described in our previous work [16], thermal conductivity, λ, electrical resistivity, ρ, and Seebeck coefficient, S, of subunit C1 in the three principle directions can be calculated as: where x1, y1, z1 are coordinates of C1 as shown in Figure 1.Subscripts F and M denote the properties of the fiber phase and the matrix phase, respectively, and n is defined as: Using the principles of the Kirchhoff Circuit Laws (KCL) and following the steps described in our previous work [16], thermal conductivity, λ, electrical resistivity, ρ, and Seebeck coefficient, S, of subunit C 1 in the three principle directions can be calculated as: where x 1 , y 1 , z 1 are coordinates of C 1 as shown in Figure 1.Subscripts F and M denote the properties of the fiber phase and the matrix phase, respectively, and n is defined as: where f is the volume fraction of the fiber phase in the unit cell.Similarly, the effective properties of the entire unit cell can be calculated from the effective properties of C 1 and the properties of C 2 , such that: Energies 2017, 10, 1006 where x 2 , y 2 , z 2 are parallel to x 1 , y 1 , z 1 , respectively (Figure 1).It is important to note that in this step, subunit C 1 is treated as a homogeneous body with its properties calculated from Equation (1) as previously described.
Using Equation (3), the effective properties in the directions parallel and perpendicular to the fibers are calculated at room temperature (300 K) for a hypothetical composite material fabricated from Bi 0.5 Sb 1.5 Te 3 as the matrix and Cu as the fiber phase (Figure 2).Material properties are listed in Table 1.As shown in Figure 2, both the thermal conductivity and electrical conductivity in the parallel direction increase linearly with increasing volume fraction of the Cu fiber, while the conductivities in the perpendicular direction remain relatively unchanged.On the other hand, the Seebeck coefficient in the parallel direction decreases rapidly from 30 µV/K at f = 0.01 to 5 µV/K at f = 0.1 and remains lower than 5 µV/K with increasing fiber volume fraction.In contrast, the Seebeck coefficient in the perpendicular direction approaches that of the matrix phase (~220 µV/K, Table 1) at low fiber volume fractions, and decreases slightly with increasing f.
where f is the volume fraction of the fiber phase in the unit cell.Similarly, the effective properties of the entire unit cell can be calculated from the effective properties of C1 and the properties of C2, such that: where x2, y2, z2 are parallel to x1, y1, z1, respectively (Figure 1).It is important to note that in this step, subunit C1 is treated as a homogeneous body with its properties calculated from Equation (1) as previously described.
Using Equation ( 3), the effective properties in the directions parallel and perpendicular to the fibers are calculated at room temperature (300 K) for a hypothetical composite material fabricated from Bi0.5Sb1.5Te3as the matrix and Cu as the fiber phase (Figure 2).Material properties are listed in Table 1.As shown in Figure 2, both the thermal conductivity and electrical conductivity in the parallel direction increase linearly with increasing volume fraction of the Cu fiber, while the conductivities in the perpendicular direction remain relatively unchanged.On the other hand, the Seebeck coefficient in the parallel direction decreases rapidly from 30 μV/K at f = 0.01 to 5 μV/K at f = 0.1 and remains lower than 5 μV/K with increasing fiber volume fraction.In contrast, the Seebeck coefficient in the perpendicular direction approaches that of the matrix phase (~220 μV/K, Table 1) at low fiber volume fractions, and decreases slightly with increasing f.In order to confirm the validity of the KCL-based model, finite element analysis (FEA) was conducted using commercial FEA software (COMSOL Multiphysics, COMSOL Inc., Burlington, MA,  In order to confirm the validity of the KCL-based model, finite element analysis (FEA) was conducted using commercial FEA software (COMSOL Multiphysics, COMSOL Inc., Burlington, MA, USA).A unit cell similar to that shown in Figure 1 was constructed.Periodical boundary conditions, including flux and potential continuity for both electrical and thermal properties, were applied to simulate the infinitely large material.Bi 0.5 Sb 1.5 Te 3 and Cu (Table 1) were used as the matrix and fiber materials.Fiber volume fraction, f, was varied between 0.01 and 0.7.In addition to fibers with square cross-sections, fibers with circular cross-sections were also studied.As shown in Figure 2, all FEA results of both square fibers and circular fibers agree fairly well with the KCL model predictions in both directions parallel and perpendicular to the fibers.
Figure 2 also shows the comparison to the effective conductivities estimated using the effective medium approach (EMA).Using the equations developed by Nan et al. [34] for predicting anisotropic thermal conductivity of composites with ellipsoidal inclusions, the effective thermal and electrical conductivities of the Bi 0.5 Sb 1.5 Te 3 matrix-Cu fiber composite were calculated by extending the length-to-diameter ratio of the ellipsoidal inclusion to infinity.As shown in Figure 2, the KCL and EMA results also compare well.
As previously discussed, a transverse TE element can be obtained by rotating the composite structure shown in Figure 1 by an angle of α about the x (or y) axis, where 0 < α < 90 • (Figure 3).If the new coordinates are x, y, z, the TE properties can then be calculated as: mwhere P = λ, ρ, and S. Assuming that the electric current flows in the y direction and the heat flux flows in the z direction, a transverse figure of merit, Z trans , can be defined as [35]: where S yz = 1 2 (S z 2 − S y 2 ) sin 2α is the transverse Seebeck coefficient, λ zz = λ z 2 cos 2 α + λ y 2 sin 2 α is the thermal conductivity in the z direction, and ρ yy = ρ y 2 cos 2 α + ρ z 2 sin 2 α is the electrical resistivity in the y direction.
Energies 2017, 10, 1006 5 of 11 USA).A unit cell similar to that shown in Figure 1 was constructed.Periodical boundary conditions, including flux and potential continuity for both electrical and thermal properties, were applied to simulate the infinitely large material.Bi0.5Sb1.5Te3and Cu (Table 1) were used as the matrix and fiber materials.Fiber volume fraction, f, was varied between 0.01 and 0.7.In addition to fibers with square cross-sections, fibers with circular cross-sections were also studied.As shown in Figure 2, all FEA results of both square fibers and circular fibers agree fairly well with the KCL model predictions in both directions parallel and perpendicular to the fibers.Figure 2 also shows the comparison to the effective conductivities estimated using the effective medium approach (EMA).Using the equations developed by Nan et al. [34] for predicting anisotropic thermal conductivity of composites with ellipsoidal inclusions, the effective thermal and electrical conductivities of the Bi0.5Sb1.5Te3matrix-Cu fiber composite were calculated by extending the length-to-diameter ratio of the ellipsoidal inclusion to infinity.As shown in Figure 2, the KCL and EMA results also compare well.
As previously discussed, a transverse TE element can be obtained by rotating the composite structure shown in Figure 1 by an angle of α about the x (or y) axis, where 0 < α < 90° (Figure 3).If the new coordinates are x, y, z, the TE properties can then be calculated as: where P = λ, ρ, and S. Assuming that the electric current flows in the y direction and the heat flux flows in the z direction, a transverse figure of merit, Ztrans, can be defined as [35]:

Cooling Performance
It can be shown that the coefficient of performance (COP) of a transverse TE device when operated in the refrigeration mode can be calculated as [3,35]:

Cooling Performance
It can be shown that the coefficient of performance (COP) of a transverse TE device when operated in the refrigeration mode can be calculated as [3,35]: Energies 2017, 10, 1006 6 of 11 where j y is the electrical current density in the y direction, T h is the heat sink temperature, ∆T is the temperature difference between the cooling surface and the heat sink, d is the device thickness, and L is the device length (Figure 3).A derivation of Equation ( 6) assumes a one-dimensional electrical current flow in the y direction and a heat flow in the z direction, which requires that the transverse TE element has a large length-to-thickness ratio (L/d) [17].The maximum temperature difference, ∆T max , is achieved when the COP approaches zero [16,35]: The ∆T max and the dimensionless figure of merit (Z trans T) of the Bi 0.5 Sb 1.5 Te 3 matrix-Cu fiber composite were calculated using Equations ( 6) and ( 7) with a fixed hot-side temperature (T h ) of 300 K.As shown in Figure 4, both ∆T max and Z trans T initially increase with increasing α (rotation angle, see Figure 3), and after reaching the peak values start to decrease as α further increases.For a given α, the effect of fiber volume fraction, f, is rather weak, especially for f > 0.1.This result implies that a moderate amount of cooling (∆T max > 30 K) can be theoretically obtained within a wide range of fiber volume fractions (0.1 < f < 0.7).
where jy is the electrical current density in the y direction, Th is the heat sink temperature, ΔT is the temperature difference between the cooling surface and the heat sink, d is the device thickness ,and L is the device length (Figure 3).A derivation of Equation ( 6) assumes a one-dimensional electrical current flow in the y direction and a heat flow in the z direction, which requires that the transverse TE element has a large length-to-thickness ratio (L/d) [17].The maximum temperature difference, ΔTmax, is achieved when the COP approaches zero [16,35]: The ΔTmax and the dimensionless figure of merit (ZtransT) of the Bi0.5Sb1.5Te3matrix-Cu fiber composite were calculated using Equations ( 6) and ( 7) with a fixed hot-side temperature (Th) of 300 K.As shown in Figure 4, both ΔTmax and ZtransT initially increase with increasing α (rotation angle, see Figure 3), and after reaching the peak values start to decrease as α further increases.For a given α, the effect of fiber volume fraction, f, is rather weak, especially for f > 0.1.This result implies that a moderate amount of cooling (ΔTmax > 30 K) can be theoretically obtained within a wide range of fiber volume fractions (0.1 < f < 0.7).The cooling performance, in terms of the ΔTmax, of the Bi0.5Sb1.5Te3matrix-Cu fiber composite was also investigated using FEA.The total volume of the model was 2 mm (x) × 500 mm (y) × 10 mm (z), with periodical boundary conditions applied in the x direction.An aspect ratio (L/d) of 50 was selected based on our previous study [16], which indicated that the end effect of the device on the temperature distribution became relatively small at a large device length-to-thickness ratio (L/d > 30).Five cases with fiber volume fractions ranging from 0.1 to 0.5 were examined.In each case, the diameter of the fibers was allowed to change to obtain different fiber volume fractions, and in all cases, the L/d ratios were greater than 40.The rotation angle was fixed at 10°, the hot-side temperature was fixed at 300 K, and the electrical current density in the y direction (jy, see Figure 3) was varied until ΔTmax in the z direction was found.
As an example, Figure 5 plots the temperature distribution when f = 0.3 for both square fibers (Figure 5a) and circular fibers (Figure 5b), where temperature gradients were established between the top (cold-side) and bottom (hot-side) surfaces.There exist some non-uniform regions near the ends of the device (near y = 0, or 500 mm) due to the device end effect as previously discussed [16].The cooling performance, in terms of the ∆T max , of the Bi 0.5 Sb 1.5 Te 3 matrix-Cu fiber composite was also investigated using FEA.The total volume of the model was 2 mm (x) × 500 mm (y) × 10 mm (z), with periodical boundary conditions applied in the x direction.An aspect ratio (L/d) of 50 was selected based on our previous study [16], which indicated that the end effect of the device on the temperature distribution became relatively small at a large device length-to-thickness ratio (L/d > 30).Five cases with fiber volume fractions ranging from 0.1 to 0.5 were examined.In each case, the diameter of the fibers was allowed to change to obtain different fiber volume fractions, and in all cases, the L/d ratios were greater than 40.The rotation angle was fixed at 10 • , the hot-side temperature was fixed at 300 K, and the electrical current density in the y direction (j y , see Figure 3) was varied until ∆T max in the z direction was found.
As an example, Figure 5 plots the temperature distribution when f = 0.3 for both square fibers (Figure 5a) and circular fibers (Figure 5b), where temperature gradients were established between the top (cold-side) and bottom (hot-side) surfaces.There exist some non-uniform regions near the ends of the device (near y = 0, or 500 mm) due to the device end effect as previously discussed [16].The insets in Figure 5a, b illustrate the temperature distribution on the top surface in the uniform region, where small local temperature variations can be observed between the fiber and the matrix.The insets in Figure 5a, b illustrate the temperature distribution on the top surface in the uniform region, where small local temperature variations can be observed between the fiber and the matrix.Table 2 compares the ΔTmax obtained from FEA simulation and the KCL model predictions for the five cases, where the rotational angle was fixed at 10°.The ΔTmax values in FEA cases were obtained by averaging the cold-side surface temperatures in the uniform region (Figure 5).In general, the two sets of data show good agreement.The ΔTmax decreases with increasing fiber volume fraction, f, except for f = 0.5 where a small increase in ΔTmax was observed.This is because ΔTmax is a nonlinear function of fiber volume fraction at a fixed rotational angle (Figure 4b).A larger discrepancy between the KCL model and the FEA results was observed at f = 0.1, indicating that the KCL approach may be less effective to model composites with low volume fiber fractions.The simplified KCL model assumes the unit cell as a homogeneous body and ignores any localized current at the fiber/matrix boundary.However, the localized current exists due to the dissimilar nature of the two phases and circulates at the phase boundary, which leads to additional Joule heating and reduces the cooling efficiency.Consequently, the ΔTmax values obtained from FEA, Table 2 compares the ∆T max obtained from FEA simulation and the KCL model predictions for the five cases, where the rotational angle was fixed at 10 • .The ∆T max values in FEA cases were obtained by averaging the cold-side surface temperatures in the uniform region (Figure 5).In general, the two sets of data show good agreement.The ∆T max decreases with increasing fiber volume fraction, f, except for f = 0.5 where a small increase in ∆T max was observed.This is because ∆T max is a nonlinear function of fiber volume fraction at a fixed rotational angle (Figure 4b).A larger discrepancy between the KCL model and the FEA results was observed at f = 0.1, indicating that the KCL approach may be less effective to model composites with low volume fiber fractions.The simplified KCL model assumes the unit cell as a homogeneous body and ignores any localized current at the fiber/matrix boundary.However, the localized current exists due to the dissimilar nature of the two phases and circulates at the phase boundary, which leads to additional Joule heating and reduces the cooling efficiency.Consequently, the ∆T max values obtained from FEA, which took into account the localized circulating current, are always lower than the KCL predictions (Table 2).The effect of the circulating current is Energies 2017, 10, 1006 8 of 11 reduced as the composite becomes more homogeneous at higher fiber volume fraction.It is noted that the FEA result of ∆T max on the circular fiber for f = 0.2 is higher than that of the square fiber (Table 2).This is possibly caused by numerical inaccuracies encountered during simulation.Using Equations ( 3)-( 5), the maximum Z trans T values were calculated for various combinations between some state-of-the-art thermoelectric materials and metals (Table 1).As first suggested by Goldsmid when discussing lamella structures [4], a good transverse TE composite may be constructed by selecting two phases with large contrasts in their transport properties.This principle should also be applicable to fibrous composites.In this study, the thermoelectric materials are selected as the matrix phases, and the metals are selected as the fiber inclusions (Table 3).Should the opposite configuration be used, the metals will form a continuous network and the thermoelectric effects will diminish due to the highly conductive nature of the metals.For all three operating temperatures (300, 700 and 1000 K) included in Table 2, a fibrous inclusion of Ni seem to give the highest Z trans T values.Au, Ag, and Cu as the fiber phase show comparable contributions.At 300 K, other metals such as Al, Bi, In, Sn, and Pb also lead to non-trivial transverse figure of merits.Although the maximum Z trans T values listed in Table 2 are lower than the corresponding TE matrix materials, which is also true for lamella structures [4], transverse TE devices may have some practical advantages.For example, a transverse TE element only uses one type of semiconductor, i.e., p-type or n-type, while conventional TE devices require both p-type and n-type materials that possess compatible properties.On the other hand, the TE performance of the transverse elements could be further improved by introducing porosity to the less conducting phase [4], which are the thermoelectric materials in the cases included in Table 3, to sharpen the contrast in conductivities between the matrix and the fiber phases.

Conclusions
In this study, a mathematical model based on the Kirchhoff Circuit Laws (KCL) was developed to describe the effective thermoelectric properties, including thermal and electrical conductivities and Seebeck coefficient, of composite materials with aligned, uniformly distributed fiber inclusions.Comparisons to finite element analysis (FEA) results and EMA models (note: there is no EMA model for the Seebeck coefficient of fiber composites) indicate that the KCL-based model can predict the effective thermoelectric properties of fiber composite materials with relatively good accuracy.Unlike FEA, which requires a significant amount of computing effort, the simplified KCL model may be useful as a quick tool for material screening and prototype device design.
The cooling performance of composite materials utilizing Bi 0.5 Sb 1.5 Te 3 as the matrix and Cu as the fiber phase was examined using both the KCL-based model and FEA, which also exhibited agreement.Moderate cooling (>30 K) near room temperature (300 K) may be theoretically achieved using this composite structure.Furthermore, predictions using the KCL approach indicate that non-trivial values of figure of merit (0.25-0.42) may be obtained when combining common thermoelectric materials and metals.
Acknowledgments: This work is supported by the Temple University Faculty Start-Up Fund.
Author Contributions: Fei Ren and Bosen Qian conceived and designed the project; Bosen Qian performed the mathematical modeling and finite element simulation; Bosen Qian and Fei Ren analyzed the data and wrote the paper.

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

Figure 1 .
Figure 1.Schematic of a unit cell used in the KCL model.

Figure 1 .
Figure 1.Schematic of a unit cell used in the KCL model.

Figure 2 .
Figure 2. Thermal conductivity, electrical resistivity, and Seebeck coefficient as a function of fiber volume fraction in directions parallel and transverse to the fiber inclusions. where resistivity in the y direction.

Figure 3 .
Figure 3. Schematic showing a piece of the fibrous composite tilted on an angle α about the x-axis.

Figure 3 .
Figure 3. Schematic showing a piece of the fibrous composite tilted on an angle α about the x-axis.

Figure 4 .
Figure 4. KCL model results of the Bi0.5Sb1.5Te3matrix-Cu fiber composite as a function of fiber volume fraction and the angle of rotation: (a) ZtransT at 300 K, and (b) the maximum ΔT when the hot-side temperature is fixed at 300 K.

Figure 4 .
Figure 4. KCL model results of the Bi 0.5 Sb 1.5 Te 3 matrix-Cu fiber composite as a function of fiber volume fraction and the angle of rotation: (a) Z trans T at 300 K, and (b) the maximum ∆T when the hot-side temperature is fixed at 300 K.

Table 1 .
Material properties used in this study.

Table 1 .
Material properties used in this study.

Table 2 .
Comparison of ∆T max between the KCL model and FEA simulations.

Table 3 .
Maximum Z trans T values for various combinations of thermoelectric matrices and metal fibers.