Jiyoung
Jung
^{a},
Seung Hee
Jeong
*^{b},
Klas
Hjort
^{b} and
Seunghwa
Ryu
*^{a}
^{a}Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology (KAIST), 291 Daehak-ro, Yuseong-gu, Daejeon 34141, Republic of Korea. E-mail: ryush@kaist.ac.kr
^{b}Division of Microsystems Technology, Department of Engineering Sciences, The Angstrom Laboratory, Uppsala University, Box 534, SE 751 21, Uppsala, Sweden. E-mail: seunghee.jeong@angstrom.uu.se

Received
17th February 2020
, Accepted 2nd May 2020

First published on 6th May 2020

For the facile use of liquid metal composites (LMCs) for soft, stretchable and thermal systems, it is crucial to understand and predict the thermal conductivity of the composites as a function of liquid metal (LM) volume fraction and applied strain. In this study, we investigated the effective thermal conductivity of LMCs based on various mean-field homogenization frameworks including Eshelby, Mori–Tanaka, differential and double inclusion methods. The double inclusion model turned out to make the prediction closest to the experimental results in a wide range of LM volume fractions. Interestingly, we found that the theoretical models based on the assumption of ideal LM dispersion and zero interfacial resistance underestimated the thermal conductivity compared to the experimental results in a low volume fraction regime. By considering the accompanied variations in the LM inclusion's aspect ratios under a typical size distribution of inclusions (∼μm), the change of effective thermal conductivity was predicted under a uniaxial 300% tensile strain. Our study will deepen the understanding of the thermal properties of LMCs and support the designs of stretchable thermal interfaces and packaging with LMCs in the future.

For the facile use of LMCs, it is important to understand and predict their thermal conductivity as a function of the LM volume fraction and applied strain, which would allow one to design LMCs for a specific application. One can employ numerical tools, such as finite element analysis, but multiple calculations are required to obtain statistically meaningful data by varying the positions and the number of inclusions within the representative volume elements (RVEs). As an alternative, theoretical models with cost-efficient analysis have been proposed for predicting the effective properties of composites. Mean field schemes such as the Maxwell–Garnett's model^{7} and Bruggeman's theory,^{8–10} differential schemes such as the Tavangar's model,^{11,12} and fast Fourier transform schemes^{13,14} have been applied to predict the effective thermal conductivity of composites such as LMCs. In this paper, we have considered the mean-field homogenization scheme based on the micromechanics theory and Eshelby's solution because of its convenience of handling a composite involving multiple types of inclusions (here, liquid inclusions have different sizes and different aspect ratios). Because of the mathematical analogy between the governing equations of elasticity and thermal conductivity, the Eshelby's theory on micromechanics has been adapted and applied to study the effective thermal properties.^{15,16}

In this work, we have employed four different homogenization frameworks, namely, Eshelby,^{17,18} Mori–Tanaka,^{15,16,19} differential^{20–22} and double inclusion methods^{20,23} to predict the thermal conductivity of the LMC as a function of the LM volume fraction and an applied uniaxial strain. These methods are classified according to how the interaction between inclusions is (partially) incorporated. In the first half, we show that the double inclusion method matches the experimental results with the smallest error in a wide range of LM volume fraction and discuss why it is so. We then investigated the unusual discrepancy between the mean-field theoretical prediction and experiments, namely, the underestimation of thermal conductivity even in the idealized LM droplet dispersion and zero interfacial thermal resistance assumptions. We hypothesize the origin of the discrepancy in terms of the LM aspect ratio as a candidate cause of the discrepancy and a few plausible explanations on the discrepancy are discussed. In the second half of the paper, we predict the change of the thermal conductivity of LMC under tension, by predicting the aspect ratio change of LM inclusions. Interestingly, the thermal conductivity of LMC is shown to change by almost a factor of three. We conclude the paper by discussing the implication and future direction of the present study.

∇·q(x) − g(x) = 0, | (1) |

q(x) = Ke(x) where e(x) = −∇T(x), | (2) |

K_{1}(e^{0} + e^{pt}) = K_{0}(e^{0} + e^{pt} − e*), | (3) |

(4) |

e^{pt} ≡ Se*. | (5) |

Obtaining the proper Eshelby tensor is the most important in applying the homogenization framework based on the Eshelby's solution. The Eshelby tensor is a function of material properties of a matrix and the shape of an inclusion. In heat conduction problems, the Eshelby tensor with an anisotropic matrix and ellipsoidal inclusion is derived as follows,^{16}

(6) |

S_{11} = 1 − 2S_{22}, | (7) |

(8) |

(9) |

When a single inclusion exists in the infinite matrix, substituting eqn (5) into (3), and the local concentration tensor T, which maps the exterior intensity field to the intensity field within the inclusion e_{1}, can be obtained as follows,

T = [I + SK_{0}^{−1}(K_{1} − K_{0})]^{−1}, | (10) |

e_{1} = e^{0} + e^{pt} = Te^{0} | (11) |

By extending Eshelby's solution on a single inhomogeneity to multiple inhomogeneity problems with the appropriate approximations, various homogenization techniques have been developed, such as the Eshelby method, Mori–Tanaka method, the differential method, the double inclusion method and others. Considering the volume average over composites, the volume averaged heat flux and intensity field satisfy the following equations,

〈q〉 = c_{0}〈q_{0}〉 + c_{1}〈q_{1}〉 = c_{0}K_{0}〈e_{0}〉 + c_{1}K_{1}〈e_{1}〉, | (12) |

〈q〉 = K_{eff}〈e〉, | (13) |

K_{eff} = K_{0} + c_{1}(K_{1} − K_{0})A where 〈e_{1}〉 = A〈e〉, | (14) |

K^{Eshelby}_{eff} = K_{0} + c_{1}(K_{1} − K_{0})T. | (15) |

Since the Eshelby method assumes an isolated inclusion, the theoretical prediction by the Eshelby method is known to be valid in a very dilute volume fraction of inclusion.

As an alternative to the Eshelby method, the Mori–Tanaka method reflects interactions between inclusions by assuming that an intensity field within a matrix around each inclusion is equal to the volume averaged intensity field of the matrix (e_{0}(x) = 〈e_{0}〉). The global concentration tensor and the effective thermal conductivity tensor in the Mori–Tanaka method are derived as follows,

A ≈ T(c_{0}I + c_{1}T)^{−1}, | (16) |

K^{MT}_{eff} = K_{0} + c_{1}(K_{1} − K_{0})T(c_{0}I + c_{1}T)^{−1} = (c_{0}K_{0} + c_{1}K_{1}T)(c_{0}I + c_{1}T)^{−1}. | (17) |

The Mori–Tanaka method is known to be valid for up to 20% of the volume fraction.

The differential method was developed to predict the effective properties of composites with moderate finite volume fractions. The differential method is proposed with the idea that the prediction from the Eshelby/Mori–Tanaka methods is reliable at low volume fractions. The differential scheme homogenization proceeds as follows. First, a small amount (Δc_{1}) of inclusions is added to the matrix with K_{0} and homogenization is applied, resulting in a homogenized phase with K^{(1)}_{eff}. Second, a small amount (Δc_{1}) of inclusions is added to the homogenized phase and homogenization is applied again. Third, the previous step is repeated until a desired volume fraction is reached and an effective thermal conductivity by the differential method K^{Diff}_{eff} is obtained. Here, a small enough incremental volume fraction Δc_{1} is chosen by testing the convergence over Δc_{1} reduction. Applying the differential method based on the Mori–Tanaka approach, the effective thermal conductivity tensor is derived as follows,

(18) |

(19) |

The double inclusion method was developed to predict the effective properties of composites with arbitrary volume fractions. When inclusions with low volume fractions are sparsely distributed as shown in Fig. 1(b), the effective properties can be predicted by applying the Mori–Tanaka method. Similarly, with eqn (10), a local concentration tensor at a low volume fraction is derived as follows,

T_{low} = [I + SK_{0}^{−1}(K_{1} − K_{0})]^{−1}. | (20) |

On the other hand, when inclusions with a very high-volume fraction are densely distributed as shown in Fig. 1(c), the matrix phase seems to be surrounded by inclusions. Considering the configuration of the matrix phase wrapped by inclusions, the local concentration tensor with a high-volume fraction is reduced as follows,

T_{high} = [I + SK_{1}^{−1}(K_{0} − K_{1})]. | (21) |

In the double inclusion scheme, the local concentration tensor T_{D} is obtained by smoothly interpolating T_{low} and T_{high’} using an appropriate interpolation function ζ. The interpolation function ζ should satisfy the following conditions,^{20}

(22) |

We set the volume fraction of inclusions as an interpolation function, and the local concentration tensor in the double inclusion method is represented as follows,

ζ(c_{1}) = c_{1}, | (23) |

T_{D} = (1 − ζ)T_{low} + ζT_{high}. | (24) |

Then, the effective thermal conductivity tensor in the double inclusion method K^{DI}_{eff} is obtained as follows,

K^{DI}_{eff} = (c_{0}K_{0} + c_{1}K_{1}T_{D})(c_{0}I + c_{1}T_{D})^{−1}. | (25) |

Fig. 2 Different homogenization models for effective thermal conductivities of LMCs. (a) Effective thermal conductivities compared to the experimental data.^{6} The thermal conductivities are normalized by the thermal conductivity of PDMS. (b) A detailed plot from 0 to 20% of volume fraction. |

LM volume fraction | 5.0% | 13.7% | 32.2% | 58.7% | 66.1% |
---|---|---|---|---|---|

Eshelby | 11.49% | 20.74% | 43.23% | 67.12% | 77.41% |

Mori–Tanaka | 10.91% | 17.28% | 30.76% | 39.72% | 50.88% |

Differential | 9.98% | 10.54% | 9.06% | — | — |

Double inclusion | 10.32% | 13.39% | 13.41% | 16.19% | 10.16% |

Interestingly, we find that the theoretical models based on ideal LM dispersion and zero interfacial resistance assumptions underestimate the thermal conductivity compared to the experimental results in a low volume fraction regime, as shown in Fig. 2(b). A similar phenomenon was observed in the previous studies predicting the Young's modulus of LMC.^{30,31} One can predict a few plausible causes of the discrepancy, such as (i) the non-spherical shape of the inclusions, (ii) the change of inherent property of PDMS upon fabrication or measurement variation, (iii) the effect of a thin oxide layer or interfacial thermal resistance and (iv) the limitation of mean-field homogenization.

First, we can exclude the third hypothesis because an interfacial thermal resistance would lower the effective thermal conductivity prediction even more and the conductivity of a thin oxide layer is comparable with that of the LM droplet. It is obvious that the mean-field homogenization framework cannot account for the effect of the local agglomeration of inclusions, but the agglomeration effect would lower the thermal conductivity. Because the theoretical model assumes the ideal dispersion of inclusions, the underestimation from the mean-field homogenization theory cannot be explained this way.

Second, regarding to the first hypothesis, we tested the effect of realistic inclusion shapes that deviate from a perfect sphere. According to the CT scan measurements on LMCs, due to the influence of a solid oxide layer, the LM inclusions actually have an ellipsoidal shape with different aspect ratios.^{9} The average aspect ratio was considered as the reported 1.5.^{9} We predicted the effective thermal conductivities of the LMCs with randomly oriented ellipsoidal inclusions in order to test whether the experimental value can be explained by accounting for realistic inclusion shape distribution. Following the previous studies of the orientation averaged Mori–Tanaka method,^{32,33} we can predict the thermal conductivity of LMCs with randomly oriented LM droplets as follows,

(26) |

(27) |

T_{ij}′ = r_{ip}r_{jq}T_{pq}, | (28) |

(29) |

LM volume fraction | 5.0% | 13.7% |
---|---|---|

ρ = 1 | 10.91% | 17.28% |

ρ = 1.5 | 10.29% | 15.92% |

To confirm our prediction, we conducted a finite element analysis (FEA) of the LMC involving randomly oriented ellipsoidal inclusions with an aspect ratio of 1.5 by generating ten different representative volume elements (RVEs) by Digimat-FE software, as shown in Fig. 3(b). The numerical results from the FEA show good agreement with the orientation averaged Mori–Tanaka method below 20% of the volume fraction, which implies that the theoretical framework is valid and also that the experimental data cannot be explained despite accounting for a realistic aspect ratio.

Finally, we tested the last possible cause of the discrepancy, which is the change of the inherent thermal conductivity of PDMS upon the fabrication of LMCs. We found that if the thermal conductivity of PDMS within LMCs increases by 10% (from 0.17 W m^{−} K^{−} to 0.187 W m^{−} K^{−}), the theoretical prediction matches better with the experimental data (see the SI.3, ESI†). In that sense, we suspect that the thermal conductivity of PDMS might have changed after fabrication, or varied from the measurements, especially the PDMS conductivity on its own and in LMCs. This can be correlated with the measured thermal property variations of PDMS. This would also explain the observation in the previous studies, which show that the experimentally measured stiffness of LMC is higher than the mean-field homogenization theory prediction.^{30,31}

(30) |

Considering the volume conservation of an axisymmetric ellipsoidal inclusion , the aspect ratio ρ can be obtained as follows,

(31) |

We considered the typical LM inclusion size distribution, ranging from 2 μm to 50 nm of the diameters of inclusions, as reported in a previous study.^{6} We predict the effective thermal conductivity of the LMC with the 20% volume fraction, based on the double inclusion method because the double inclusion method showed the best match with experiments among four different methods in the previous section. Because we consider the relatively small volume fraction, the aspect ratio prediction for a single inclusion, eqn (30), can be applied to estimate the aspect ratio change of LM inclusions within LMC, which is initially spherical inclusions. The Eshelby tensor S for an ellipsoidal inclusion is obtained by substituting eqn (31) into (7) and (8). Then, the Eshelby tensor S is substituted into eqn (20) and (21). The effective thermal conductivity tensor from the double inclusion method is obtained by substituting eqn (20) and (21) to eqn (24) and (25). The effective thermal conductivities of LMCs with 2 μm and 50 nm inclusions along the tensile direction are shown in Fig. 4(a) (red and yellow lines, respectively). The results are normalized by the thermal conductivity at zero strain, κ_{ε=0}. The effective thermal conductivity for a large inclusion (2 μm) increases significantly with applied tensile strain. On the other hand, the variation of the effective thermal conductivity for smaller inclusions (50 nm) is small since the surface tension tends to keep the inclusion spherical. The dashed lines in Fig. 4(a) are limiting case results for zero surface tension and rigid inclusion, respectively. The thermal conductivity of the LMC converges to the upper limit, that of a laminated composite at a very high strain limit (see SI.4, ESI†).

Fig. 4 Tensile strain effect on effective thermal conductivities of LMCs. (a) Effective thermal conductivities in the tensile direction under the far-field strain at 20% of volume fraction. Thermal conductivities are normalized by thermal conductivity at zero strain. The upper limit in the high strain regime is represented as the blue dotted line (see SI.4, ESI† for details). (b) The assumed log-normal distribution of LM inclusions. (c) Normalized effective thermal conductivities in a transverse direction to the tensile strain under the far-field strain at 20% of volume fraction. (d) Aspect ratios of stretched inclusions of different sizes as a function of the far-field tensile strain. |

According to the SEM images of the previous study, inclusion sizes have a log-normal distribution.^{9,38} A log-normal distribution of inclusions is assumed as shown in Fig. 4(b). The blue histogram shows the fraction of inclusions within each bin and the yellow line shows the corresponding volume fraction. When the distribution of inclusions has a histogram of N bins, the effective thermal conductivity tensor of the LMC is represented as follows,

(32) |

(33) |

We note that the contribution of the oxide layer is not accounted for in our theoretical framework. Because the oxide layer observed on the surface of LM acts like a thin and rigid shell,^{9,39,40} it would give rise to a decreased surface tension effect on the properties of the LMCs. Quantitative analysis on the complex interplay between the oxide layer and the liquid surface tension, and the resulting effective properties of LMCs will remain an important task for the LMC community. Also, because eqn (30) is derived for a single inclusion embedded in a matrix subjected to a far-field strain, the accuracy of our model is expected to be limited in a high strain regime when the distance between inclusions significantly reduces due to the Poisson effect.

We also predicted that the effective thermal conductivity of typical LMCs, having micrometer scale inclusions with a log-normal distribution, increases by almost a factor of three along the tensile direction at 300% strain. The sensitivity of thermal conductivity to strain can be tuned by controlling the inclusion size distributions because it originates from the elastocapillary effect. Our study deepens the understanding on the thermal conductivity of LMCs and provides guidelines for designing and implementing thermal interfaces and packaging for applications. Still, the rigorous considerations of modeling the parameters of LMCs as a complex system, including the oxide skin layers between inclusions and a matrix, thermal interface resistance, and orientation distribution, will challenge future research.

- J. A. Rogers, T. Someya and Y. Huang, Science, 2010, 327, 1603–1607 CrossRef CAS PubMed.
- N. Lu and D.-H. Kim, Soft Robot., 2014, 1, 53–62 CrossRef.
- M. Amjadi, K. U. Kyung, I. Park and M. Sitti, Adv. Funct. Mater., 2016, 26, 1678–1698 CrossRef CAS.
- M. H. Malakooti, N. Kazem, J. Yan, C. Pan, E. J. Markvicka, K. Matyjaszewski and C. Majidi, Adv. Funct. Mater., 2019, 29, 1906098 CrossRef CAS.
- C. Pan, E. J. Markvicka, M. H. Malakooti, J. Yan, L. Hu, K. Matyjaszewski and C. Majidi, Adv. Mater., 2019, 31, 1900663 CrossRef PubMed.
- S. H. Jeong, S. Chen, J. Huo, E. K. Gamstedt, J. Liu, S.-L. Zhang, Z.-B. Zhang, K. Hjort and Z. Wu, Sci. Rep., 2015, 5, 18257 CrossRef CAS PubMed.
- P. Fan, Z. Sun, Y. Wang, H. Chang, P. Zhang, S. Yao, C. Lu, W. Rao and J. Liu, RSC Adv., 2018, 8, 16232–16242 RSC.
- R. Tutika, S. H. Zhou, R. E. Napolitano and M. D. Bartlett, Adv. Funct. Mater., 2018, 28, 1804336 CrossRef.
- M. D. Bartlett, A. Fassler, N. Kazem, E. J. Markvicka, P. Mandal and C. Majidi, Adv. Mater., 2016, 28, 3726–3731 CrossRef CAS PubMed.
- M. D. Bartlett, N. Kazem, M. J. Powell-Palm, X. Huang, W. Sun, J. A. Malen and C. Majidi, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 2143–2148 CrossRef CAS PubMed.
- M. I. Ralphs, N. Kemme, P. B. Vartak, E. Joseph, S. Tipnis, S. Turnage, K. N. Solanki, R. Y. Wang and K. Rykaczewski, ACS Appl. Mater. Interfaces, 2018, 10, 2083–2092 CrossRef CAS PubMed.
- R. Tavangar, J. M. Molina and L. Weber, Scr. Mater., 2007, 56, 357–360 CrossRef CAS.
- F. Nosouhi Dehnavi, M. Safdari, K. Abrinia, A. Sheidaei and M. Baniassadi, Mater. Today Commun., 2020, 23, 100878 CrossRef CAS.
- D. J. Eyre and G. W. Milton, Eur. Phys. J.: Appl. Phys., 1999, 6, 41–47 CrossRef.
- H. L. Quang, Q.-C. He and G. Bonnet, Philos. Mag., 2011, 91, 3358–3392 CrossRef CAS.
- S. Lee, J. Lee, B. Ryu and S. Ryu, Sci. Rep., 2018, 8, 1–11 CrossRef PubMed.
- J. D. Eshelby, Proc. R. Soc. London, Ser. A, 1957, 241, 376–396 Search PubMed.
- J. Qu and M. Cherkaoui, Fundamentals of micromechanics of solids, Wiley, Hoboken, 2006 Search PubMed.
- J. Jung, S. Lee, B. Ryu and S. Ryu, Int. J. Heat Mass Transfer, 2019, 144, 118620 CrossRef.
- I. Doghri and A. Ouaar, Int. J. Solids Struct., 2003, 40, 1681–1712 CrossRef.
- R. McLaughlin, Int. J. Eng. Sci., 1977, 15, 237–244 CrossRef.
- A. Norris, Mech. Mater., 1985, 4, 1–16 CrossRef.
- M. Hori and S. Nemat-Nasser, Mech. Mater., 1993, 14, 189–206 CrossRef.
- H. Hiroshi and T. Minoru, Int. J. Eng. Sci., 1986, 24, 1159–1172 CrossRef.
- M. Hodes, R. Zhang, L. S. Lam, R. Wilcoxon and N. Lower, IEEE Trans. Compon., Packag., Manuf. Technol., 2013, 4, 46–56 Search PubMed.
- J. E. Mark, Polymer data handbook, Oxford University Press, 2009 Search PubMed.
- Z. Guo, A. Verma, X. Wu, F. Sun, A. Hickman, T. Masui, A. Kuramata, M. Higashiwaki, D. Jena and T. Luo, Appl. Phys. Lett., 2015, 106, 111909 CrossRef.
- G. J. Weng, Mech. Mater., 2010, 42, 886–893 CrossRef.
- S. Torquato and H. Haslach Jr, Appl. Mech. Rev., 2002, 55, B62–B63 CrossRef.
- R. W. Style, R. Boltyanskiy, B. Allen, K. E. Jensen, H. P. Foote, J. S. Wettlaufer and E. R. Dufresne, Nat. Phys., 2015, 11, 82–87 Search PubMed.
- F. Mancarella, R. W. Style and J. S. Wettlaufer, Proc. R. Soc. London, Ser. A, 2016, 472, 20150853 CrossRef PubMed.
- S. Lee and S. Ryu, Eur. J. Mech. A:Solids, 2018, 72, 79–87 CrossRef.
- G. M. Odegard, T. Gates, K. Wise, C. Park and E. Siochi, Compos. Sci. Technol., 2003, 63, 1671–1687 CrossRef CAS.
- R. W. Style, A. Jagota, C.-Y. Hui and E. R. Dufresne, Annu. Rev. Condens. Matter Phys., 2017, 8, 99–118 CrossRef CAS.
- R. W. Style, J. S. Wettlaufer and E. R. Dufresne, Soft Matter, 2015, 11, 672–679 RSC.
- S. H. Jeong, A. Hagman, K. Hjort, M. Jobs, J. Sundqvist and Z. Wu, Lab Chip, 2012, 12, 4657–4664 RSC.
- V. Kocourek, C. Karcher, M. Conrath and D. Schulze, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 026303 CrossRef CAS PubMed.
- Z. J. Farrell and C. Tabor, Langmuir, 2018, 34, 234–240 CrossRef CAS PubMed.
- H. Liang, Y. Zeng, K. Zuo, X. Xia, D. Yao and J. Yin, Ceram. Int., 2017, 43, 5517–5523 CrossRef CAS.
- E. Ricci, E. Arato, A. Passerone and P. Costa, Adv. Colloid Interface Sci., 2005, 117, 15–32 CrossRef CAS PubMed.

## Footnote |

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sm00279h |

This journal is © The Royal Society of Chemistry 2020 |