Liquid metal phase change materials for thermal management of electronics

ABSTRACT

With the rapid development of various chips towards high performance and integration, the ‘thermal barrier’ difficulty facing electronic devices has become ever increasingly challenging, while at the same time the space size left for thermal management system already becomes much narrower. In view of the intermittent operation or short-time large power consumption of some electronic devices, the thermal management technology based on phase change materials (PCMs) has received more and more attention. Among those cutting edge PCMs, the liquid metal phase change materials (LMPCMs) especially have aroused much interest due to their outstanding merits in thermal conductivity, energy storage density and stability. In this article, the representative works on LMPCMs are comprehensively reviewed. First, the thermophysical parameters of LMPCMs characterized by high thermal conductivity and large density are summarized. Then, some basic research cases on the thermophysical properties of metallic PCMs that are worthy of further investigation are briefly introduced. In view of the serious supercooling problem of the LMPCMs, the factors affecting the supercooling degree of liquid metal (LM), particularly the thermal history effect, based on the nucleation theory, are presented. Strategies for suppressing supercooling using nucleating agents and crystal species as well as external fields are also introduced. In addition, the mathematical models and related numerical algorithms for interpreting the phase change heat transfer mechanisms of LM are illustrated. Finally, the performance advantages of LMPCMs over other PCMs, as well as the technical approaches to further improve various properties of LMPCM heat sinks, are demonstrated by showing typical application cases. Based on the above, the current research status of LMPCMs and promising future research hotspots are prospected.This review is expected to help researchers interested in LMPCMs to quickly grasp the performance, heat transfer evaluation methods, and research status of LMPCMs, and also to provide a guidance for engineers in related fields to conduct further design and practices.

Keywords:

• Liquid metal
• phase change materials
• thermal management
• electronics

1. Introduction

In this day and age, the rapidly increasing integration and miniaturization of various electronic device poses ever severe challenges to thermal management technology. Due to the operational limitations, the temperature inside the electronic devices should not exceed 85°C [1]. For every 1% of temperature exceeding the maximum limit, the life of the devices is reduced by 5% [2]. According to statistics, more than 55% of electronic equipment failures are caused by high temperatures at local hot spots where heat is not dissipated in time [3]. The problem of heat dissipation has become one of the key issues restricting the continuous development of electronic devices. The current evolution of electronic devices has brought some new challenges to thermal management technology: (1) Ultra-high heat flux density; (2) Small and portable electronic devices have limited space and are not suitable for installing active cooling system; (3) The power of continuously operating electronic device may also fluctuate wildly. If a continuously working cooling system is designed based on the power peak of the thermal management object, it will inevitably cause redundancy in power consumption, volume and cost. These challenges are not independent of each other and may occur simultaneously in some application scenarios. For example, active devices and pulsed laser devices require a thermal management system that can not only maintain the normal temperature of the equipment during continuous operation, but also meet the heat dissipation requirements during high-power operation in a short period of time within a limited volume and mass [4].

The current techniques applied to thermal management in electronic devices are: air cooling, liquid cooling and phase change material (PCM) cooling [5]. Air cooling technology was first proposed. Although the technology is mature and the structure is simple, air cooling has limited heat transfer capacity and slow cooling speed, so it is not suitable for handling heat flux densities greater than 100W/cm² [6]. Liquid cooling has larger convective heat transfer coefficient than air cooling, which can better control the temperature of electronic devices. However, as an active cooling technology, it is usually bulky, noisy, consumes extra pump power [7,8], and it cannot cope well with short-time ultra-high thermal impulse.

PCM cooling, as a passive thermal management technique, does not require additional drive components or additional power consumption [9]. PCM can absorb or release vast amounts of latent heat during phase change processes while maintaining a constant temperature, providing a unique opportunity to develop thermal management systems with high cooling capacity [4]. At the same time, PCM can effectively cope with thermal shocks over short periods of time [10,11].

Generally, PCM can be categorized into organic PCM and inorganic PCM based on their chemical composition, where organic PCMs include paraffins, fatty acids and polyols, and inorganic PCMs cover-hydrated salts, inorganic salts and metal [12,13]. Paraffins are currently the most widely used organic PCM with a range of advantages such as high latent heat, good chemical stability, non-toxicity, no phase separation and almost no supercooling [14,15]. However, their low thermal conductivity (~0.2 W/(m·K)) [16] severely limits thermal management capability. As inorganic PCM, hydrated salts have been attracting much attention due to their reasonable price, wide availability, good thermal conductivity and high energy storage density [17], but defects such as corrosivity, supercooling and phase separation have limited their practical applications [18,19]. Liquid metal (LM), is a class of metal PCM with melting points near room temperature (<100°C in this review). Due to a series of properties such as high thermal conductivity (20–60 W/(m·K)) [20], high energy storage density, non-toxicity and good stability, LM has received more and more attention in the field of thermal management [21,22]. In addition to direct application in thermal management, liquid metal phase change materials (LMPCM) can be further enhanced with built-in fins [23,24], carbon-based [25] or metallic porous materials [22]. LMPCM can be encapsulated into phase change capsules to stabilize their shape as well as to prevent leakage while further enhancing their thermal management performance [26]. LM phase change capsule can also be added to liquid carrier to prepare phase change slurries [27], utilizing their high thermal conductivity and latent heat to achieve a better thermal management performance.

It is worth noting that LM have a very serious supercooling effect, and the maximum supercooling degree of bulk gallium at atmospheric pressure can reach about 55–70°C [28,29]. If the phase change material does not solidify in time due to overcooling during the working cycle, it will not be able to accommodate the heat generated by the heat source during the melting process in the next working cycle [30], causing the heat source temperature to continue to rise out of control. LM represented by gallium is corrosive to common metal materials such as copper and aluminum [31,32], and as the temperature increases, the corrosion phenomenon becomes more and more serious [33]. In addition, since LMs are prone to leakage and have good electrical conductivity [5,34], leakage of LM can cause serious problems such as short circuits in electronic devices.

In order to deepen the understanding of LMPCMs and promote the further application on LMPCMs, the properties, heat transfer models and representative applications of LMPCM are reviewed in this article. Firstly, the thermophysical properties of LM as well as the supercooling effect are introduced, and the supercooling phenomenon is analyzed based on the nucleation theory. The Methods of suppressing supercooling with the help of nucleating agents and external fields are also summarized. Then the LM phase change heat transfer models are described, and two numerical simulation methods, the Enthalpy-porosity method and the Lattice Boltzmann Method (LBM), are illustrated. Typical application studies on LMPCM, as well as cases of improving the comprehensive performance of LMPCM through phase change composite (PCC) and combined heat transfer methods, are also introduced. Finally, the challenges and prospectives of LMPCM for future electronic thermal management applications are summarized.

2. Properties of liquid metal phase change materials

2.1. Thermophysical properties of liquid metal phase change materials

Table 1 lists the thermophysical properties of some PCMs and PCCs. The high thermal conductivity of LMPCMs brings huge advantages in the field of thermal management, which is difficult for other PCMs to achieve even after enhancing heat transfer through various methods. Within the same melting point range, although the latent heat of LMPCMs is relatively low, due to the high density, its energy storage density is considerable. This feature makes LMPCM particularly suitable for working in confined space like the inside of electronic devices. In addition, LMPCMs also have advantages in thermal stability and phase change volume change. Unlike most materials, the density of gallium increases during melting, which means there is no need to design too much space margin to avoid the leakage of melting gallium caused by expansion. As a result, higher filling amount of LM can be achieved, which is beneficial for the performance of heat sink. It is also possible to utilize the abnormal volume change of LM during the phase change process to combine with ordinary PCM to form zero volume change PCM. Generally, the thermal conductivity and latent heat of gallium-based LM are higher, and the melting point and density of bismuth-based LM are higher.

Table 1. The thermophysical properties of some phase change material and phase change composite [47–50].

PCM	Melting Point (℃)	Latent Heat of Fusion （kJ/kg）	Thermal Conductivity (W/(m·K))	Specific Heat Capacity (J/(kg·K))	Density (kg/m3)	Volumetric Latent Heat of Fusion (MJ/m3)
Ga67In20.5Sn12.5	10.7	67.2	23.4	-	6170	415
Eutectic Ga-In (EGaIn)	15.7	69.7	26.58(l)	404(l)	6197	432
Gallium (Ga)	29.8–30	80.1–80.3	33.5–40.6(s) 28.1–33.7(l)	340–370(s) 390–400(l)	5903–5907(s) 5907–6096(l)	473
Bi44.7Pb22.6In19.1Sn8.3Cd5.3	47–48	28.5–36.8	16.7(s)	150–163(s)	8860–9330	253
Bi47.5Pb25.4In5Sn12.6Cd9.5	57–65	36	15	160(s)	9470	341
Bi49In21Pb18Sn12	58–61	23–29.9	7.143–33.2(s) 10.6(l)	167–323(s) 211(l)	9307(s) 8220(l)	218
Bi31.6In48.8Sn19.6	60.2	27.9	19.2(s) 14.5(l)	270(s) 297(l)	8043	224
Bi32In51Sn17	60.1–60.7	27.7	19(s)	197	8057	223
Bi33In51Cd16	61	25	-	-	8000–10040	200
Bi50Pb26.7Sn13.3Cd10 (Wood’s metal)	70–72	29.6–45.8	16.7–23	146–340(s) 170–184(l)	8974–9620	266
Eutectic Bi-In	72–72.8	20.4–25	-	-	7990–10000	163
Bi57In26Sn17	77–80.7	32.5–47.7	35.8(s) 28.8(l)	400(s) 880(l)	8780	285
Bi54In30Sn16	90.9	38	-	-	8470	322
Bi50Pb31Sn19	94.5	26.4	16.5	138	9590	253
Octadecane	29	244	0.358(s) 0.152(l)	2150(s) 2180(l)	814(s) 724(l)	177
Heneicosane	41	294.9	0.145(l)	2386(l)	773(l)	228
Tricosane	48.4	302.5	0.124(l)	2181(l)	777.6(l)	235
Tetracosane	51.5	207.7	0.137(l)	2914(l)	773.6(l)	161
Palmitic Acid	64	185	0.162(l)	2200(s) 2480(l)	989(s) 850(l)	157
Stearic Acid	69	202	0.172(l)	2830(s) 2380(l)	965(s) 848(l)	171
CaCl2·6H2O	29	170–192	1.008(s) 0.561(l)	-	1802(s) 1562(l)	266
Na2SO4·10H2O	32	251	0.544	-	1485(s)	373
Na2HPO4·12H2O	35	265	0.514(s) 0.476(l)	-	1552(l)	411
Ba(OH)2·8H2O	78	265.7	1.255(s) 0.653(l)	-	2180(s) 1937(l)	515
Paraffin+EG	37	160.53–196.65	0.324–3.355	1900–2260	-	-
PEG+BN+EG	40.42	152.7	0.289	-	-	-
Paraffin+GR+Ni Foam	-	99	4.6	1600(s) 1800(l)	-	-

Note: s refers to solid state, l refers to liquid state, respectively.

EG: Expanded Graphite; BN: Boron Nitride; GR: Graphene.

As PCM, LM also has many shortcomings. Although the high density of LM brings advantages in thermal storage density, it is a disadvantageous factor for its application in weight-sensitive situations such as vehicles, portable and wearable devices. Another major drawback of LMPCMs is that their cost is much higher than that of other organic and inorganic PCMs, which seriously limits the practical engineering applications. In addition, the compatibility between LMPCMs and metal containers is also one of the limiting factors.

The thermophysical properties of LMPCMs are affected by many factors, including composition, size, pressure and even casting methods. Understanding these factors is very important for the fundamental research and practical application of LMPCMs. So far, there is still much less effort on investigating the thermophysical properties of LMPCMs compared to those of high melting point metal PCMs. Generally, the property laws of the same type of substance are common, so in certain extent one can find out the clues from the existing endeavors on the thermophysical properties of high melting point metals.

The research on metal PCMs started early in the field of high-temperature heat storage, and researchers also found some rules in the process of studying the thermophysical properties of high melting point alloys. In the study of alloys composed of elements such as Al, Mg, Cu, Si, Zn, Birchenall pointed out that alloys containing high melting point elements usually have higher latent heat [35,36]. Fang et al. found that the thermal conductivity of Mg-Bi alloy decreased with the increase of Bi content, and the higher the eutectic phase ratio of the alloy, the higher the melting latent heat of the alloy [37].

In the study of alloys whose melting point is closer to the thermal management temperature zone of electronic equipment, Wang et al. found that the alloy containing Zn has more micro-interfaces in the fractography when researching Sn-Bi-Pb-Zn low melting point alloy (melting point: 90.78°C–171.75°C) [38]. These interfaces have a great influence on the thermophysical properties of the alloy. More micro-interfaces mean that more metal bonds need to be destroyed in the melting process, so the addition of Zn to Sn-Bi-Pb alloy can improve the latent heat of fusion. In addition to the latent heat, the specific heat capacity of the alloy increases with the increase of Zn content. The saturated vapor pressure of Zn is higher than that of other metal component, so the thermal stability of the alloy decreases with the increase of Zn, but it can still remain stable below 900°C. Besides, the higher the content of Zn, the higher the coefficient of thermal expansion. The addition of Pb can form an intermetallic compound (Pb₇Bi₃) with the high melting point component Bi, which reduces the content of Bi in the alloy and leads to the decrease of the melting point of the alloy. The density of the alloy increases with the increase of the content of high-density components Pb and Bi. Due to the existence of Bi, the density of Sn-Bi-Sn alloy increases after melting, but with the addition of Pb, the formation of intermetallic compound eliminates this phenomenon. The thermal conductivity decreases with the increase of low thermal conductivity component Bi and increases with the increase of high thermal conductivity component Zn. It is not difficult to see that, the properties of the alloy are highly related to its composition.

Zhou et al. characterized the microstructure, melting point, latent heat and phase composition of 13 kinds of ternary, quaternary and quaternary eutectic alloys composed of Bi, Gd, Zn, Pb with melting point between 48.3°C and 124°C [39]. It is found that the micro-morphology of these low melting point metals is very complex, including lamellar, rod, irregular net and other structures, which may relate to the different mixing enthalpy between component elements. According to the test results of Zhang et al., the melting point and melting latent heat of quaternary and quaternary alloys are generally lower than those of ternary alloys.

For the thermophysical properties of LMs with lower melting point, the latent heat of gallium-based alloy and bismuth-based alloy increases with the increase of gallium and bismuth content, respectively. Duan et al. mixed Ga₆₇In_20.5Sn_12.5 (melting point: 10.5°C enthalpy: 60 J/g density: 6.37 g/cm³) and Bi_33.1In_51.3Sn_15.6 alloy (melting point: 60°C enthalpy: 27.9 J/g density: 7.9 g/cm³) according to different mass ratios to form Ga-Bi-In-Sn quaternary alloy system [40]. The quaternary alloy is in the form of paste at room temperature, and the melting point and density are between the two raw material alloys. The melting point and density decrease with the increase of gallium content.

Fu et al. reviewed the research on the melting point of LM and concluded that the melting point of LMs is affected by alloy composition, pressure and size effect [41]. Generally speaking, the melting point of materials increases with the increase of pressure, and the relationship between melting point and pressure can be described by Clausius-Clapeyron equation [42]:

(1) $\frac{\mathit{dP}}{\mathit{dT}}=\frac{\mathrm{\Delta }{H}_f}{\mathit{T}\Delta \mathit{V}}=\frac{\mathrm{\Delta S}}{\mathrm{\Delta V}}$(1)

where, P is the pressure, T is the temperature, H_f is the latent heat, V is the volume, and S is entropy. Jayaraman et al. measured the melting point variation of gallium with pressure, and found that the melting point does not change linearly with pressure [43], which is related to the different stable state of gallium at different pressure levels: the melting point of α-Ga decreases abnormally with the increase of pressure, while the melting point of Ga (II) and Ga (III) increases with the increase of pressure. The abnormal change of melting point of α-Ga with pressure may be related to its complex crystal structure.

When the size of gallium is reduced to the nano level, its melting point is obviously affected by the size. The gallium melting point with 0.45 μm diameter is 27.9°C, which is about 2 K lower than that of bulk gallium, while the gallium melting point with 35 nm particle size is as low as −14.2°C [44].

The main mechanisms of metal heat conduction include electronic heat conduction and lattice vibration, and the former is dominant. The thermal conductivity of electrons can be expressed by the following formula:

(2) ${\lambda }_e=\frac{1}{3}{c}_e{v}_e{l}_e$(2)

where λ_e is the electron thermal conductivity, c_e is the unit volume-specific heat capacity of electron, v_e is the average electron velocity, and l_e is the electron average free path. Equation (2)(2) ${\lambda }_e=\frac{1}{3}{c}_e{v}_e{l}_e$(2) is a general expression for the electron thermal conductivity of metals which is also applicable to liquid metal phase change materials when they are in solid state. In addition to the melting point, according to Huang et al., the thermal conductivity of metals also has size effect [45]: the electron thermal conductivity and phonon thermal conductivity of metal decrease with the decrease of particle diameter, and the former is more sensitive to the size change. When the size of metal nanoparticles is close to the order of average electron-free path, the thermal conductivity decreases. The supercooling of PCMs also has size effect, which will increase with the decrease of size. Therefore, LMPCMs may not be suitable for making into nano-capsules like other PCMs.

In the study of high melting point alloy, it is found that the casting method of the alloy has an effect on the thermophysical properties of the alloy. For example, in the study of 7075 aluminum alloy, Li et al. found that liquid squeeze casting and semi-solid squeeze casting will lead to the increase of grain boundaries, which lead to the scattering of electrons through the grain boundaries, the shortening of the average electron-free path and the formation of thermal resistance [46]. As a result, the thermal conductivity of squeeze casting alloy is about 10% lower than that of mold casting alloy. For LMPCMs which often undergo solid–liquid phase change, although it is impossible to design the so-called casting process, however, in theory, it is possible to control the cooling conditions of LM to make the alloy produce fewer defects in the solidification process to strive for higher solid-state thermal conductivity, but there is a lack of related research so far.

2.2. Supercooling effect of liquid metal

The supercooling effect means that a material does not solidify when its temperature is cooled below its melting point T_m. The difference between the temperature T of the molten liquid and its melting point T_m under a certain pressure is the degree of supercooling ΔT. Since gallium has a typical supercooling effect, many scholars have chosen gallium as the research object in the studies of crystal nucleation. Turnbull et al. [51] investigated the effects of size effect, thermal history effect, inter-atomic interfacial energy, and crystal structure on LM nucleation, and preliminarily constructed a theory of LM nucleation. Kumar et al. [52] used DSC to study the thermal behavior of gallium particles and bulk gallium, and found that the supercooling effect of gallium particles was more pronounced than that of metal blocks. Yamaguchi et al. [53] used ultrasonication for the reversible control of gallium nanoparticles, and experimentally measured the solidification point of the GaNPs prepared at 20°C to be −128.3°C, and the melting point of the GaNPs to be −14.2°C. A study by Cicco et al. [54] found that gallium encapsulation in micron to nano-sized particles lowers the freezing and melting points, and sub-micron Ga droplets in epoxy resin in the experiments could be cryogenically cooled down to 150 K without undergoing a phase change.

In this Section 2.2, we initially provide a concise overview of the LM nucleation theory and then analyses the supercooling phenomenon of gallium based on the nucleation theory. Subsequently, Section 2.3 details various methods to inhibit the effects of supercooling.

2.2.1. Nucleation theory

The classical nucleation theory (CNT) was proposed by Turnbull [55,56], Volmer [57], Farkas [58] and Becker [59], and has been improved and developed by many subsequent scholars. Two important thermodynamic quantities in CNT are the critical nucleation radius r_c and the critical potential barrier ΔG_c. Homogeneous nucleation occurs in liquids when the nuclei come from the growth of crystal embryo within the liquid. Assuming the presence of a spherical crystal embryo of radius r in a supercooled liquid, the free energy of the system changes during the growth of the crystal embryo:

(3) $\mathrm{\Delta G}=\frac{-4}{3}\mathit{\pi }{r}^3\mathrm{\Delta }{G}_P+4\mathit{\pi }{r}^2{\sigma }_{\mathit{LC}}$(3)

where, ΔG_P is the absolute value of the difference in free energy per unit volume between the solid and liquid phases and σ_LC is the interfacial energy per unit area between the solid and liquid phases. Assuming that the melting point of the substance is T_m, its supercooling ΔT = T_m - T. The relationship between ΔG_P and ΔT was given by Aleksandrov et al. [60]:

(4) $\mathrm{\Delta }{G}_P=\frac{\mathrm{\Delta }{H}_P}{T_m}\mathrm{\Delta }\mathit{T}$(4)

where, ΔH_P is the enthalpy of melting per unit volume of the substance. Substituting Equation (4)(4) $\mathrm{\Delta }{G}_P=\frac{\mathrm{\Delta }{H}_P}{T_m}\mathrm{\Delta }\mathit{T}$(4) into Equation (3)(3) $\mathrm{\Delta G}=\frac{-4}{3}\mathit{\pi }{r}^3\mathrm{\Delta }{G}_P+4\mathit{\pi }{r}^2{\sigma }_{\mathit{LC}}$(3) and deriving it, the critical nucleation radius r_c for homogeneous nucleation and the critical potential barrier ΔG_c can be obtained:

(5) $r_c=\frac{2{\mathit{\sigma }}_{\mathit{LC}}{T}_m}{\mathrm{\Delta }{H}_P\mathrm{\Delta }\mathit{T}},\mathrm{\Delta }{G}_c=\frac{16\pi {\sigma }_{LC}^3{T}_m^2}{3\mathrm{\Delta }{H}_P^2\mathrm{\Delta }{T}^2}$(5)

Heterogeneous nucleation occurs when the nuclei are created with the help of foreign materials such as impurity particles and vessel walls. A schematic diagram of heterogeneous nucleation occurring on a planar substrate is shown in Figure 1(a). There is a substrate S in the liquid L, to which a crystal embryo C is attached because of interfacial tension. θ is the three-phase contact angle between the liquid, the crystal embryo and the substrate. The critical nucleation radius $r_c^{{}^{\prime }}$ and critical potential barrier $\mathrm{\Delta }{G}_c^{{}^{\prime }}$ for heterogeneous nucleation are as follows [55,62]:

(6) $r_c^{{}^{\prime }}=\frac{2{\mathit{\sigma }}_{\mathit{LC}}{T}_m}{\mathrm{\Delta }{H}_P\mathrm{\Delta }\mathit{T}},\mathrm{\Delta }{G}_c^{{}^{\prime }}=\frac{16\pi {\sigma }_{LC}^3{T}_m^2}{3\mathrm{\Delta }{H}_P^2\mathrm{\Delta }{T}^2}\left(\frac{2-3\mathrm{cos\theta }+co{s}^3\mathit{\theta }}{4}\right)$(6)

Figure 1. Schematic representation of several different nucleation theories (a) schematic diagram of heterogeneous nucleation on a planar substrate (b) Step-function change of the free-energy density through the interfacial region between the liquid and solid (c) schematic comparison of one-step versus two-step nucleation. Reprinted with permission from ref [61]. Copyright 2016 from American Chemical Society.

The critical nucleation radius r_c and the critical potential barrier ΔG_c symbolize the difficulty of nucleation. The larger r_c and ΔG_c are, the more difficult it is for nucleation to occur, i.e. the more susceptible it is to supercooling. From Equation (5)(5) $r_c=\frac{2{\mathit{\sigma }}_{\mathit{LC}}{T}_m}{\mathrm{\Delta }{H}_P\mathrm{\Delta }\mathit{T}},\mathrm{\Delta }{G}_c=\frac{16\pi {\sigma }_{LC}^3{T}_m^2}{3\mathrm{\Delta }{H}_P^2\mathrm{\Delta }{T}^2}$(5) , when the temperature T of the liquid is lowered, both r_c and ΔG_c are lowered, and crystal nucleation is easier to form, which is the physical mechanism of classical nucleation theory that lowering the temperature can induce solidification of the liquid.

Crystal nucleation and growth are not only related to thermodynamics, but also to kinetics. The probability of an atom crossing from the solid phase to the liquid phase in response to an interfacial driving force is a kinetic factor J, expressed in the form of the Arrhenius equation [63]:

(7) $\mathit{J}=\mathit{C}{e}^{\left(\frac{-\mathrm{\Delta }{G}_A}{\mathit{kT}}\right)}$(7)

where, C is a constant, ΔG_A is the activation energy for diffusion of atoms across phase boundaries, and k is Boltzmann constant with a value of 1.38 × 10⁻²³J/K. With the above thermodynamic and kinetic parameters, Becker et al. went on to propose the nucleation rate I [9]:

(8) $\mathit{I}=\mathit{C}{e}^{\left(\frac{-\mathrm{\Delta }{G}_A}{\mathit{kT}}\right)}{e}^{\left(\frac{-\mathrm{\Delta }{G}_C}{\mathit{kT}}\right)}$(8)

Turnbull et al. [56] derived the nucleation rate I for nucleation processes that do not require long-range diffusion by combining Equation (8)(8) $\mathit{I}=\mathit{C}{e}^{\left(\frac{-\mathrm{\Delta }{G}_A}{\mathit{kT}}\right)}{e}^{\left(\frac{-\mathrm{\Delta }{G}_C}{\mathit{kT}}\right)}$(8) and absolute reaction rate theory:

(9) $\mathit{I}=\frac{\mathit{NkT}}{h}{e}^{\left(\frac{-\mathrm{\Delta }{G}_A}{\mathit{kT}}\right)}{e}^{\left(\frac{-\mathrm{\Delta }{G}_C}{\mathit{kT}}\right)}$(9)

where, N is the total number of atoms or molecules per unit volume of the liquid, and h is Planck constant with a value of 6.63 × 10⁻³⁴J∙s. The nucleation rate is a very important physical quantity in the nucleation process, and its unit is nuclei/(s∙m³). Its physical meaning is the number of nucleation events per unit time per unit volume, which is a good response to the difficulty of nucleation of a system. The CNT has been widely used since it was proposed. It has a simple form and can describe various nucleation behaviors. However, there are still some defects in its theory, such as the simpler capillary assumption adopted in the CNT to describe the physical characteristics of clusters, and its reliability when applied to small-size clusters has been controversial [64]. In order to obtain more accurate results, many scholars have proposed different nucleation theory models, such as diffusion interface theory, two-step nucleation theory, density functional theory and so on.

The diffuse interface theory was developed independently by Gránásy [65,66] and Spaepen [67]. Diffuse interface theory suggests that there exists a region with a certain thickness between the solid–liquid phases, rather than a clear boundary at the solid–liquid phase interface as suggested by CNT [68]. As shown in Figure 1(b), the free energy per unit volume of the diffusion interface is G_i, and the width is δ; the free energies per unit volume on both sides of the diffusion interface are G_s and G_l, respectively. Based on this free energy distribution, the free energy change of the system ΔG can be obtained, and then the critical nucleation radius r_c and the critical potential barrier ΔG_c can be obtained by ∂∆G/∂r = 0.

While the CNT suggests that the system only needs to overcome the critical potential barrier corresponding to a crystal nucleus of a certain critical size, the two-step nucleation theory [69–71] suggests that the process of crystal formation by molecules in solution usually follows a two-step nucleation mechanism, i.e. it needs to overcome two free energy barriers. As shown in Figure 1(c), crystal nucleation must overcome the first-free energy barrier $\mathrm{\Delta }{G}_{n_{\rho },\mathit{twostep}}^{}$ through density fluctuations of the solute, resulting in the formation of connected molecular clusters of size, at which point the object does not yet possess any crystalline order. Subsequently, the system must cross a second free-energy barrier $\mathrm{\Delta }{G}_{n,tw\mathrm{ostep}}^{}$ in order for the molecules within the dense clusters to be ordered in a manner similar to crystallization. The two-step nucleation theory has been demonstrated in the formation of sodium chloride crystals [72], protein crystallization [73], and a host of others.

Density functional theory (DFT) was introduced by Hohenberg-Kohn in 1964 [74]. The theory suggests that for non-simplified ground state molecules, the energy of the system, the electronic wave function, and the nature of the electrons are determined only by the electron density. As a broader theory for electronic system, it is one of the most commonly used methods in condensed matter physics and computational chemistry. By choosing a realistic and simple enough-free energy functional, Oxtoby et al. proposed a systematic theory of nucleation based on DFT that has been applied to several different models of liquid-gas and liquid-solid nucleation [75,76]. It can be emphasized that the theory introduced here has a specific application of DFT for nucleation in liquids. Their group pointed out the existence of a rather thick interfacial layer between the crystal nucleus and the supercooled liquid, which is not sufficiently described by CNT.

These non-classical nucleation models mentioned above tend to improve some of the accuracies while creating new deficiencies in the simplicity or generalizability of the models. The current various nucleation theories are broad in scope, and there is no complete set of theories specifically designed to provide a reasonable explanation for the nucleation behavior of LMs. The high surface tension and easy oxidation of LMs, which are different from other materials, will inevitably have an impact on their nucleation behavior, and some nucleation theories may not be applicable to them. Despite a certain degree of deviation from experimental data [77], the CNT is still the most classical and commonly used theoretical model for describing homogeneous nucleation processes due to its simplicity and predictive power [61]. Therefore, the supercooling phenomenon of LM is still analyzed based on the classical nucleation theory in the next Section 2.2.2.

2.2.2. Analysis of liquid metal supercooling

Based on the thermophysical parameters provided by Costa and Turnbull et al. [49,51], we can calculate the critical nucleation radius and critical potential barrier for different liquids according to Equation (5)(5) $r_c=\frac{2{\mathit{\sigma }}_{\mathit{LC}}{T}_m}{\mathrm{\Delta }{H}_P\mathrm{\Delta }\mathit{T}},\mathrm{\Delta }{G}_c=\frac{16\pi {\sigma }_{LC}^3{T}_m^2}{3\mathrm{\Delta }{H}_P^2\mathrm{\Delta }{T}^2}$(5) . Table 2 shows the critical nucleation radius and nucleation potential barriers for five PCMs at supercooling degrees of 10°C and 30°C, respectively. As can be seen from the table, the critical nucleation radius and nucleation potential barriers of the four LMs are higher than those of water. In the case of gallium, for example, the critical nucleation radius is 7.14 nm at a supercooling degree of 10°C, which is about 24% larger than that of water and similar to that of mercury. However, its nucleation potential barrier is 1.19 × 10⁻¹⁷J, much larger than that of mercury and water, which indicates that gallium is more difficult to nucleate and has a larger degree of supercooling. Compare the two cases of 10°C and 30°C, it can be found that the critical nucleation radius and nucleation barrier of gallium at a supercooling degree of 30°C are only 2.38 nm and 1.32 × 10⁻¹⁸ J, which are decreased by 66.6% and 88.9%, respectively. This indicates that temperature reduction can significantly promote the nucleation of gallium. The supercooling phenomenon of gallium is then explained in conjunction with the nucleation rate. In the context of a supercooled gallium liquid system, Turnbull posited that an exponent of exp(-ΔGA/kBT) = 10⁻¹² is sensible [51]. Applying this principle to Equation (9)(9) $\mathit{I}=\frac{\mathit{NkT}}{h}{e}^{\left(\frac{-\mathrm{\Delta }{G}_A}{\mathit{kT}}\right)}{e}^{\left(\frac{-\mathrm{\Delta }{G}_C}{\mathit{kT}}\right)}$(9) , it becomes apparent that at a supercooling of 10°C, the nucleation rate of gallium is virtually negligible, implying that spontaneous nucleation of pure gallium at room temperature is exceedingly improbable. At a supercooling of 30°C the nucleation rate increases to 2 × 10⁻¹¹⁹, which is still a very small value. The influence of size effects on nucleation was studied at the same time as the classical nucleation theory was proposed by Turnbull et al. [51]. The homogeneous nucleation frequency I_v for a droplet of volume v is as follows:

Table 2. Critical nucleation radius and nucleation potential barrier of five typical phase change materials at a supercooling of 10°C [49,51].

					rc(10^−9m)		ΔGc(10^−18J)
Materials	σ(10^−3J·m^−2)	ΔH(kJ·kg^−1)	Tm(K)	ρ(kg·m^−3)	ΔT=10°C	ΔT=30°C	ΔT=10°C	ΔT=30°C
Water	32.1	333.0	273.2	917	5.74	1.91	4.43	0.66
Mercury	24.4	11.4	234.3	13546	7.04	2.47	5.60	0.49
Gallium	55.9	80.3	303.0	5907	7.14	2.38	11.94	1.32
Tin	54.5	60.5	505.7	7280	12.52	4.17	35.76	3.97
Bismuth	54.4	53.3	544.0	10028	11.07	3.69	27.94	3.10

(10) $I_v=\mathit{Iv}$(10)

The equation illustrates that macroscopically large volumes of liquid gallium are more likely to nucleate at the same degree of supercooling. Probabilistically, the probability of finding a crystal embryo that effectively promotes nucleation is much higher in bulk gallium than in microdroplets. This probability increases as the droplet volume or surface area increases, which is a good explanation for the size effect in LM nucleation [18,19].

Comparison of Equations (5)(5) $r_c=\frac{2{\mathit{\sigma }}_{\mathit{LC}}{T}_m}{\mathrm{\Delta }{H}_P\mathrm{\Delta }\mathit{T}},\mathrm{\Delta }{G}_c=\frac{16\pi {\sigma }_{LC}^3{T}_m^2}{3\mathrm{\Delta }{H}_P^2\mathrm{\Delta }{T}^2}$(5) and (6(6) $r_c^{{}^{\prime }}=\frac{2{\mathit{\sigma }}_{\mathit{LC}}{T}_m}{\mathrm{\Delta }{H}_P\mathrm{\Delta }\mathit{T}},\mathrm{\Delta }{G}_c^{{}^{\prime }}=\frac{16\pi {\sigma }_{LC}^3{T}_m^2}{3\mathrm{\Delta }{H}_P^2\mathrm{\Delta }{T}^2}\left(\frac{2-3\mathrm{cos\theta }+co{s}^3\mathit{\theta }}{4}\right)$(6) ) reveals the relationship between the critical nucleation radius and the nucleation potential barrier for heterogeneous nucleation and homogeneous nucleation:

(11) ${r_c}^{{}^{\prime }}=r_c,\mathrm{\Delta }{G}_c^{{}^{\prime }}=\mathrm{\Delta }{G}_c\left(\frac{2-3\mathrm{cos\theta }+co{s}^3\mathit{\theta }}{4}\right)$(11)

From Equation (9)(9) $\mathit{I}=\frac{\mathit{NkT}}{h}{e}^{\left(\frac{-\mathrm{\Delta }{G}_A}{\mathit{kT}}\right)}{e}^{\left(\frac{-\mathrm{\Delta }{G}_C}{\mathit{kT}}\right)}$(9) , it is known that for $\forall \mathit{\theta f}\left(\mathit{\theta }\right)=\left(2-3\mathit{cos\theta }+\mathit{co}{s}^3\mathit{\theta }/4\right)⩽1$, which means the nucleation barrier of heterogeneous nucleation is always smaller than that of homogeneous nucleation, so heterogeneous nucleation is more easily realized. When the crystal C can achieve complete attachment to the substrate S, i.e. $\mathit{\theta }=0^{\circ }$, then $\mathit{f}\left(\mathit{\theta }\right)=0$，$\mathrm{\Delta }{G}_c^{{}^{\prime }}=0$, indicating that the crystal is extremely easy to nucleate at this time; when the crystal C lacks any attachment to substrate S, i.e. $\mathit{\theta }={180}^{\circ }$, then $\mathit{f}\left(\mathit{\theta }\right)=1$，$\mathrm{\Delta }{G}_c^{{}^{\prime }}=\mathrm{\Delta }{G}_c$, indicating that the difficulty of heterogeneous nucleation at this moment is equivalent to that of homogeneous nucleation. Based on the above analysis, it becomes evident that introducing impurity particles into the liquid can reduce the nucleation potential barrier. Furthermore, the efficacy of this process hinges on the adhesive properties between the supercooled liquid and the impurity particles, with a stronger adhesion effect facilitating easier crystal nucleation.

2.2.3. Thermal history effect of liquid metal supercooling

The previous Section 2.2.2 analyzed various factors affecting LM supercooling, such as size effect and heterogeneous nucleation. In this Section, the thermal history effect is introduced as another important factor affecting LM supercooling. The thermal history effect means that the maximum supercooling of liquid increases with increasing heating temperature over a certain temperature range. Othmer [78] first observed this effect and summarized a large number of experiments to conclude that anisotropic molecules may still exist in liquids formed from melting solids. Their number naturally decreases as the temperature difference between the liquid temperature and the melting point increases, and eventually they seem to disappear altogether. Richards et al. [79] carried out a comprehensive study of the thermal history effect of gallium, concluding that the irrelevant structure retains a ‘crystalline adsorbate’ above the melting temperature. This adsorbate is thermally stabilized by negative adsorption on the foreign structure in the stable region of the liquid. The stability of the crystalline adsorbate decreases as the temperature of the liquid increases, and eventually disappears completely at some higher temperature. Webster [80] observed that the solidification of gallium is strongly related to the thermal history effect, and that gallium crystals grow very rapidly in experiments, at least 1 cm/s at 0°C. Turnbull et al. [55] gave a theoretical explanation for the thermal history effect: a he presence of retained crystalline embryos in cavities or foreign particles in the vessel wall at a certain heating temperature, as shown in Figure 2(a). These embryos may become nucleation centers during solidification, thus increasing the nucleation rate. And at each temperature, there is a critical value r above which no embryos are retained in the pores.

Figure 2. Theoretical explanation and experimental results of the thermal history effect. (a) Retention of embryos on conical cavity and cylindrical cavity [55]. (b) DSC curves with different scanning ranges at a scanning rate of 15°C/min [62]. (c) The relationship of the maximum temperature T_max and supercooling ΔT [62].

Due to the existence of the thermal history effect, it is necessary to determine the appropriate heating temperature in order to maximize the possibility of eliminating the influence of the thermal history effect when studying the degree of supercooling of metal PCMs represented by gallium. In the study of Aleksandrovsk et al. [60], the degree of supercooling during solidification of liquid gallium increases continuously with the superheating temperature at the time of melting, and is independent of any further increase in superheating after a certain limiting value (about 34 K) has been reached. At the same time, as long as the liquid gallium is superheated above the melting temperature by 20–50 K, it will have almost no effect on the solidification temperature of gallium. Zhang et al. [62] investigated the effect of the cooling rate on the supercooling degree of gallium by varying the scanning rate of DSC, and the results of the study is shown in Figure 2(b,c). When the heating temperature is above 70°C, the maximum supercooling degree of gallium is no longer changed, with the maximum supercooling degree of 68°C.

2.3. Approaches to inhibit the supercooling effect of liquid metal

2.3.1. Adding nucleating agents or crystal seeds

According to the previous analysis, an important means to inhibit the LM supercooling effect is to add impurity particles to promote heterogeneous nucleation. Impurity particles are categorized into crystal seeds and nucleating agents, the former being particles with the same composition as the supercooled liquid, and the latter having a composition different from that of the supercooled liquid.

The nucleation effect of crystal seed is excellent and is the easiest way to induce nucleation of the supercooled liquid [81]. Yu et al. [82] placed silicon crystals, graphite, and solid gallium on the surface of supercooled gallium, respectively. As shown in Figure 3(a), when a piece of solid gallium with tiny size was placed on the surface of supercooled liquid gallium, liquid gallium started to crystallize rapidly around the gallium seed. And after some time, the gallium droplet was transformed into a solid ball. Pan et al. [85] simulated the solidification of liquid gallium by molecular dynamics calculations. The simulation results show that the addition of gallium seeds can effectively control the supercooling of liquid gallium, and the crystal structure also has a certain effect on the supercooling. Wang et al. [86] overcame the supercooling effect by using solid gallium wire as nucleating agents in contact with liquid gallium, and the solidification point of gallium containing metal particles (~20°C) was measured by DSC to be several tens of degrees higher than that of pure gallium at a cooling rate of 5°C/min. In addition, its enthalpy of phase change was much lower than that of pure gallium, implying that the addition of metal particles could make gallium easier to solidify. Afterwards, Wang et al. [83] systematically investigated the effect of four different fine wires with a diameter of 300 µm Ga/Cu/Ni/Fe as nucleating agents on the solidification point of liquid gallium. Figure 3(b) illustrates that using solid gallium as a nucleating agent can increase the solidification point from −36.3°C to 16.2°C. In contrast, when using iron or copper as the nucleating agent, the freezing point of gallium would only show a small increase, to −20.7°C and −23.5°C, respectively. In contrast, nickel was completely unable to reduce the supercooling of liquid gallium. The explanation given in the article for the different effects produced by the four wires is that particles with a similar structure to the liquid gallium crystal structure would be beneficial in promoting solidification. Although the nucleation of crystalline seeds is effective, they have the disadvantage that they need to be re-added each time they solidify, so they are not suitable for applications such as PCMs that require multiple melting and solidification.

Figure 3. (a) Crystallization process of supercooled gallium of supercooled gallium liquid metal (LM) after a solid gallium seed was placed on the surface of the LM gallium [82]. (b) Freezing point of liquid Ga without nucleating agent (red region) or with nucleating agent [83]. (c) Relationship between supercooling of Ga and equilibration temperature [84]. (d) Evolution of the supercooling reduction for different proportion and types of nucleating agents [62].

The application of nucleating agents is more common in the application of PCMs [87], where nucleating agents with a higher melting point than the PCM should be selected to ensure that the nucleating agent remains solid during the melting process. Otherwise, the nucleating agent cannot have the effect of promoting heterogeneous nucleation during solidification of the PCM. The results of Turnbull et al. [88] showed that the most efficient nucleation catalysts have a symmetric structure and lattice spacing similar to that of the forming crystals. Telks also suggests that the nucleating agent should have a cell size variation of less than 15% compared to that of the PCM in order to allow the PCM to grow on the structure of the nucleating agent [89]. Shamberger’s group [84] used an epitaxial lattice-matching technique to find two cubic carbides and nitrides with lattice mismatches of less than 0.05 relative to (010) Ga. Figure 3(c) shows that, at all equilibrium temperatures, the supercooling of gallium with the addition of cubic rock salt structural carbides and nitrides was less than 20°C and less than 10°C, respectively, the smallest supercooling reported to date. Other cubic carbide and nitride phases with higher lattice mismatches were found to lead to an increase in the degree of supercooling of Ga. Ma et al. [90] found experimentally that the primary surface oxides of BiInSnZn alloys can act as surface nucleation sites and thus inhibit the supercooling effect of their particles, and the inhibition effect increases with the amount of primary surface oxides. Zhang et al. [62] selected five particles, TeO₂, CaO, MgO, Fe and Cu, as nucleating agents for their experiments based on the selection rules proposed by Telks and the infiltration between nucleating agents and gallium. As shown in Figure 3(d), the lattice constant is one of the most important parameters when selecting the nucleating agent. The addition of TeO₂, which has a lattice constant most similar to that of α-Ga, has the best effect on reducing the supercooling of gallium, and the maximum supercooling of gallium is the lowest at about 38.2°C with the addition of 0.5 wt% TeO₂. In addition, better wettability of the nucleating agent with gallium also helped to suppress the supercooling.

2.3.2. Applying external field

Besides introducing impurity particles into the supercooled liquid to promote heterogeneous nucleation and inhibit supercooling, the external field can also inhibit the supercooling of the liquid.

Magnetic field can change the intermolecular force and molecular barrier, and then change the macroscopic properties of the liquid, which affects the crystallization of supercooled liquid. Zhuang et al. [91] studied the effect of high magnetic field on the crystallization process of amorphous alloy Fe₈₃B₁₀C₆Cu. It was found that a magnetic field of 12T greatly promotes the precipitation kinetics of α-Fe and thus enhances the crystallization process of the amorphous alloy. Stiller et al. [92] concluded that the moving magnetic field can exert an influence on the convection patterns within the liquid phase during the solidification of molten metals, consequently affecting the resulting crystal structures.

The application of an electric field to a solution or a crystalline material changes the chemical potential and thus the chemical potential difference between the liquid and the solid. Therefore, electric field can change the nucleation work by changing the chemical potential difference [93]. Li et al. [94] investigated the effect of electric field on the structure, point defects and Raman spectra of ZnO. The results showed that under high-voltage static electric field conditions, the complete crystallization time of zinc oxide was significantly prolonged and the temperature was significantly increased. Dong et al. [95] developed a novel electro-nucleation technique that utilizes a tiny DC voltage to trigger the crystallization of a supercooled sodium acetate trihydrate solution. It was investigated that increasing the voltage can shorten the induced nucleation time within a certain voltage range. Xin et al. [96] used the voltage to induce rapid nucleation of supercooled gallium metal droplets in a salt solution, resulting in a realistic reversible stiffness transformation of LM-polymer materials. This reversible stiffness transition is easily triggered by a small voltage (5 V) in the electrolyte solution.

Mechanical vibration can intensify the energy rise and fall of molecules in the supercooled liquid thereby accelerating the formation of nuclei, thus favoring the formation of nuclei. Zhou et al. [97] investigated the effect of shock vibration on supercooled sodium acetate trihydrate supercooled liquid. In the experiment, steel balls with different radius were used to impact a heat storage unit filled with sodium acetate trihydrate, and the effects of different drop heights and impact positions on the degree of supercooling were investigated. The results show that with the increase of the diameter and drop height of the steel ball can shorten the time of induced nucleation. Meanwhile, the solidification of supercooled PCM is more easily triggered by the strike position near the lid or edge than near the center part. Zhang et al. [98] analyzed the solidification conditions of supercooled gallium droplets after impact. The effects of droplet Weber number, substrate temperature and substrate material on the solidification of droplets were investigated by experiments under both isothermal and non-isothermal impacts. The experimental results show that solidification of supercooled gallium droplets is more likely to occur when the droplets hit a cold copper substrate with a high Weber number. The effect of ultrasound on the nucleation of supercooled liquids was investigated by Cogné et al. [99] When the sound pressure was increased from 150 kPa to about 220 kPa, the number of nuclei generated increased with increasing supercooling and sound pressure. When an acoustic pressure amplitude of about 220 kPa was applied, bubbles with an initial radius of 5 μm could start nucleating from supercooling degrees as low as 5 K.

3. Modeling of phase change heat transfer in liquid metal

3.1. Mathematical-physical model of the Stefan problem

The physical modeling of phase change heat transfer is generally simplified as a heat transfer problem coupled with phase change and natural convection as well as a phase-interface mobility problem, also called the Stefan problem [100]. Figure 4 shows a schematic diagram of the phase change process for a rectangular PCM, which is initially solid. The temperature of the left wall is suddenly raised to T_H at t = 0, after which the PCM starts to melt, the liquid-phase region starts to form, and the solid–liquid interface starts to move to the right. Since there is a temperature gradient in the liquid-phase region, natural convection is bound to occur. Let the melting temperature of this PCM be T_m, then the control equations are as follows, described in the form of Cartesian tensor:

Figure 4. Schematics of solid–liquid phase change problem.

Liquid phase region:

Continuity equation:

(12) $\frac{\mathit{\partial \rho }}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i\right)}{\mathit{\partial }{x}_i}=0$(12)

where, ρ is the density of the PCM, u is the velocity in the liquid-phase region of the PCM, and the subscript i (i=x, y, z) denotes the component of this vector in the i direction.

Momentum equation:

(13) $\frac{\partial \left(\rho u_i\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i{u}_j\right)}{\mathit{\partial x}}=\frac{-\mathit{\partial p}}{\partial x_i}+\mathit{\mu }\frac{{\partial }^2{u}_i}{\mathit{\partial }{x_i}^2}+\mathit{\rho }{g}_i$(13)

where, p is the pressure of the PCM, g is the gravitational acceleration, and μ is the dynamic viscosity of the PCM.

Energy Equation:

(14) $\frac{\partial \left(\rho T_l\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i{T}_l\right)}{\mathit{\partial }{x}_i}=\frac{\partial }{\partial x_i}\left(\frac{k}{{\mathit{c}}_{\mathit{p},\mathit{l}}}\frac{\partial T_l}{\mathit{\partial }{x}_i}\right)$(14)

where, k_l and c_p,l are the thermal conductivity and specific heat capacity of the liquid phase PCM, respectively.

No flow occurs in the solid-phase region, which is controlled only by the energy equation:

(15) $\frac{\partial \left(\rho T_s\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i{T}_s\right)}{\mathit{\partial }{x}_i}=\frac{\partial }{\partial x_i}\left(\frac{k}{{\mathit{c}}_{\mathit{p},\mathit{s}}}\frac{\partial T_s}{\mathit{\partial }{x}_i}\right)$(15)

where, k_s and c_p,s are the thermal conductivity and specific heat capacity of the solid phase PCM.

The governing equations for the solid–liquid interface:

(16) $-k_l\frac{\partial T_l}{\mathit{\partial x}}+k_s\frac{\partial T_s}{\mathit{\partial x}}=\mathit{\rho }\Delta \mathit{H}\frac{\mathit{ds}\left(\mathit{t}\right)}{\mathit{dt}}$(16)

(17) $T_l=T_s=T_m$(17)

where, ΔH is the enthalpy of phase change of the PCM, s(t) is the interface position, and ds(t)/dt is the moving speed of the phase interface.

Equations (16)(16) $-k_l\frac{\partial T_l}{\mathit{\partial x}}+k_s\frac{\partial T_s}{\mathit{\partial x}}=\mathit{\rho }\Delta \mathit{H}\frac{\mathit{ds}\left(\mathit{t}\right)}{\mathit{dt}}$(16) and (17(17) $T_l=T_s=T_m$(17) ) show that the existence of the solid–liquid interface requires an additional physical quantity s(t) to be solved for the interface position, which brings a lot of difficulties in solving the equations. It is worth mentioning that, in the external force-enhanced close contact phase change problem, the solid PCM which has not yet melted is pressed tightly on the surface of the heat source, and the phase change interface is always maintained on the surface of the heat source [101], then the solid–liquid interface control equation will be simplified away.

In order to solve the Stefan problem, two currently used numerical models are presented next: the Enthalpy–porosity Method and the Lattice Boltzmann Method.

3.2. Numerical simulation methods

3.2.1. Enthalpy–porosity method

In order to solve the momentum coupling problem of solid–liquid two-phase, Voller et al. [102] proposed the Enthalpy–porosity Method in 1988 by borrowing the method of porous media flow. The Enthalpy–porosity method describes the energy equation by replacing the temperature with the total enthalpy, so that the solid–liquid two-phase can be controlled by the same equation. At the same time, the paste-like region in the solid–liquid two-phase region is regarded as the porous medium region, and the porosity is expressed by the local liquid-phase fraction, so the Enthalpy–porosity Method can automatically track the interface position. It has been recognized by a wide range of scholars since it was proposed, and now has become a mainstream method to study and solve the solid–liquid phase change model [103–105]. The volume change of PCM during solid–liquid phase change is neglected, and the Boussinesq approximation is used to consider the effect of natural convection and simplify the equations. The governing equations are as follows:

The continuity Equation (12)(12) $\frac{\mathit{\partial \rho }}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i\right)}{\mathit{\partial }{x}_i}=0$(12) remains unchanged and the momentum Equation (13)(13) $\frac{\partial \left(\rho u_i\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i{u}_j\right)}{\mathit{\partial x}}=\frac{-\mathit{\partial p}}{\partial x_i}+\mathit{\mu }\frac{{\partial }^2{u}_i}{\mathit{\partial }{x_i}^2}+\mathit{\rho }{g}_i$(13) becomes:

(18) $\frac{\partial \left(\rho u_i\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i{u}_j\right)}{\mathit{\partial x}}=\frac{-\mathit{\partial p}}{\partial x_i}+\mathit{\mu }\frac{{\partial }^2{u}_i}{\mathit{\partial }{x_i}^2}-\mathit{\rho }{g}_i\mathit{\beta }\left(\mathit{T}-T_{\mathit{ref}}\right)+\mathit{A}{u}_i$(18)

where μ, β and T_ref are the kinetic viscosity, the coefficient of volume expansion and the reference temperature of the PCM, respectively. Au_i is the source term of the kinetic loss of this vector in each direction. The source term A is defined as follows:

(19) $\mathit{A}=\frac{-\mathit{C}{\left(1-\mathit{\gamma }\right)}^2}{\mathit{\epsilon }+{\gamma }^2}$(19)

where, γ is the local liquid phase fraction, C is a large number used to suppress momentum (usually taken to be 10⁵, and ε is a very small amount to prevent the denominator from being zero (usually taken to be 0.001).

The energy Equations (14)(14) $\frac{\partial \left(\rho T_l\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i{T}_l\right)}{\mathit{\partial }{x}_i}=\frac{\partial }{\partial x_i}\left(\frac{k}{{\mathit{c}}_{\mathit{p},\mathit{l}}}\frac{\partial T_l}{\mathit{\partial }{x}_i}\right)$(14) and (15(15) $\frac{\partial \left(\rho T_s\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i{T}_s\right)}{\mathit{\partial }{x}_i}=\frac{\partial }{\partial x_i}\left(\frac{k}{{\mathit{c}}_{\mathit{p},\mathit{s}}}\frac{\partial T_s}{\mathit{\partial }{x}_i}\right)$(15) ) become:

(20) $\frac{\mathit{\partial }\left(\mathit{\rho H}\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i\mathit{H}\right)}{\mathit{\partial }{x}_i}=\frac{\partial }{\partial x_i}\left(\frac{k}{c_p}\frac{\mathit{\partial H}}{\partial x_i}\right)$(20)

where, H is the total enthalpy of the PCM and the two components included, latent heat H_f and sensible heat h of the PCM, namely:

(21) $\mathit{H}=\mathit{h}+\mathit{\gamma }{H}_f$(21)

(22) $\mathit{h}=h_{\mathit{ref}}+{\int }_{{\mathit{T}}_{\mathit{ref}}}^{\mathit{T}}{c}_p\mathit{dT}$(22)

(23) $\mathit{\gamma }=\{\begin{array}{l}0,\mathit{T}<T_m\\ \frac{T-T_m}{T_l-T_m},T_m<\mathit{T}<T_l\\ 1,T_l<\mathit{T}\end{array}$(23)

where, T_m is the melting onset temperature and T_l is the temperature at the end of the melting process. h_ref is the reference enthalpy corresponding to the reference temperature T_ref.

If the natural convection during the phase change can be neglected, the convection terms in Equation (20)(20) $\frac{\mathit{\partial }\left(\mathit{\rho H}\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i\mathit{H}\right)}{\mathit{\partial }{x}_i}=\frac{\partial }{\partial x_i}\left(\frac{k}{c_p}\frac{\mathit{\partial H}}{\partial x_i}\right)$(20) can be neglected, and the whole problem becomes a simple heat conduction model for the phase change. Equation (20)(20) $\frac{\mathit{\partial }\left(\mathit{\rho H}\right)}{\mathit{\partial t}}+\frac{\partial \left(\rho u_i\mathit{H}\right)}{\mathit{\partial }{x}_i}=\frac{\partial }{\partial x_i}\left(\frac{k}{c_p}\frac{\mathit{\partial H}}{\partial x_i}\right)$(20) then degenerates into:

(24) $\frac{\mathit{\partial }\left(\mathit{\rho H}\right)}{\mathit{\partial t}}=\frac{\partial }{\partial x_i}\left(\frac{k}{c_p}\frac{\mathit{\partial H}}{\partial x_i}\right)$(24)

The Enthalpy–porosity Method perfectly unifies the solid-liquid two-phase equations, which can not only solve the solid–liquid phase change problem with a fixed grid, but also automatically capture the solid–liquid interface, which greatly simplifies the solution process of the Stefan problem. And at the same time, it is simple to operate and has small computational volume. The method has been widely recognized and used since it was proposed, and has now become one of the mainstream methods for solving solid–liquid phase change heat transfer problems. Another fixed-region algorithm similar to the Enthalpy–porosity Method is the Effective Heat Capacity Method. It considers the latent heat of phase change L of the PCM as a part of the heat capacity of the PCM, while the derivative of the enthalpy H change with temperature is defined as the equivalent heat capacity. Both of them equivalently represent the energy change during a phase change to the change of another physical quantity and are essentially the same.

3.2.2. Lattice Boltzmann method

Unlike the Enthalpy–porosity Method for solving the Navier–Stokes equations based on the continuous model introduced above, the Lattice Boltzmann Method (LBM) is a mesoscopic simulation method based on molecular dynamics theory proposed to solve the Boltzmann equation in a discrete form. It was born in the late 1980s [106], initially to solve the statistical noise problem that plagued its predecessor, the Lattice Gas Automata(LGA) [107,108]. Although proposed in a relatively short period of time, it has attracted widespread attention, particularly in complex fluid domains such as multiphase flow, viscoelastic fluids, and phase changes [109,110]. This is due to its advantages of easy boundary treatment, parallelized computational properties and high computational efficiency [111,112].

Denote the physical space infinitesimal volume by $\mathit{d}\overrightarrow{r}$ and the velocity space infinitesimal volume by $\mathit{d}\overrightarrow{c}$. Then, the statistical distribution function of molecules is defined as the number of molecules within the infinitesimal volume $\overrightarrow{r}$ to $\overrightarrow{r}+\mathit{d}\overrightarrow{r}$ at time t with velocities within the range $\overrightarrow{c}$ to $\overrightarrow{c}+\mathit{d}\overrightarrow{c}$. When the velocity distribution function of molecules is known, various macroscopic quantities in dynamics can be obtained by taking moments of f. To take moments means to multiply f by a certain function of molecular velocity and then integrate over the entire velocity space. For example, the density of a fluid ρ in physical space is:

(25) $\mathit{\rho }=\mathit{m}\int \int \int \mathit{fd}\overrightarrow{c}$(25)

The macroscopic velocity of the fluid $\overrightarrow{u}$ is:

(26) $\overrightarrow{u}=\frac{m}{\rho }\int \int \int \overrightarrow{c}\mathit{fd}\overrightarrow{c}$(26)

The change in the number of molecules in physical space $\mathit{d}\overrightarrow{r}$ with velocity in the vicinity of $\overrightarrow{c}$ within the velocity space infinitesimal volume $\mathit{d}\overrightarrow{c}$ over time $\mathit{dt}$ is $\frac{\mathit{\partial f}}{\mathit{\partial t}}\mathit{d}\overrightarrow{r}\mathit{d}\overrightarrow{c}\mathit{dt}$. The term consists of three components, the convective term $-\overrightarrow{c}\frac{\mathit{\partial f}}{\partial \overrightarrow{r}}\mathit{d}\overrightarrow{r}\mathit{d}\overrightarrow{c}\mathit{dt}$ of molecules through the surface of $\mathit{d}\overrightarrow{r}$, the convective term $-\overrightarrow{a}\frac{\mathit{\partial f}}{\partial \overrightarrow{c}}\mathit{d}\overrightarrow{r}\mathit{d}\overrightarrow{c}\mathit{dt}$ of the molecule through the surface of $\mathit{d}\overrightarrow{c}$ and the intermolecular collision term $\mathrm{\Omega }\left(\mathit{f}\right)\overrightarrow{r}\mathit{d}\overrightarrow{c}\mathit{dt}$, namely:

(27) $\frac{\mathit{\partial f}}{\mathit{\partial t}}\mathit{d}\overrightarrow{r}\mathit{d}\overrightarrow{c}\mathit{dt}=-\overrightarrow{c}\frac{\mathit{\partial f}}{\partial \overrightarrow{r}}\mathit{d}\overrightarrow{r}\mathit{d}\overrightarrow{c}\mathit{dt}-\overrightarrow{a}\frac{\mathit{\partial f}}{\partial \overrightarrow{c}}\mathit{d}\overrightarrow{r}\mathit{d}\overrightarrow{c}\mathit{dt}+\mathrm{\Omega }\left(\mathit{f}\right)\overrightarrow{r}\mathit{d}\overrightarrow{c}\mathit{dt}$(27)

Organizing the above equation gives:

(28) $\frac{\mathit{\partial f}}{\mathit{\partial t}}+\overrightarrow{c}\frac{\mathit{\partial f}}{\partial \overrightarrow{r}}+\overrightarrow{a}\frac{\mathit{\partial f}}{\partial \overrightarrow{c}}=\mathrm{\Omega }\left(\mathit{f}\right)$(28)

This equation is the Boltzman equation, where $\overrightarrow{a}$ is the acceleration produced by the external field force and $\mathrm{\Omega }\left(\mathit{f}\right)$ is the collision function. To solve the Boltzman equation, Bhatnagar, Gross and Krook introduced a simple approximate collision operator model to simplify the equation, called as the BGK operator model [113]:

(29) $\mathrm{\Omega }\left(\mathit{f}\right)=\frac{-1}{\tau }\left(\mathit{f}-f^{\mathit{eq}}\right)$(29)

where, $\mathit{\tau }$ is the relaxation time and $f^{\mathit{eq}}$ is the equilibrium distribution function, also known as Maxwell-Boltzman distribution [114]:

(30) $f^{\mathit{eq}}=\frac{\rho }{{\left(\frac{2\mathit{\pi kT}}{m}\right)}^{\frac{D}{2}}}\mathit{exp}\left(\frac{-{\left(\overrightarrow{c}-\overrightarrow{u}\right)}^2}{\frac{2\mathit{kT}}{m}}\right)$(30)

where, k is the Boltzman constant, m is the molecular mass, and u is the macroscopic flow rate of the fluid. The LBM requires discretizing the Boltzman equation and then solving for the equilibrium distribution function of the particles in each lattice cell due to migration and collisions. The equation is discretized by replacing the continuous velocity $\overrightarrow{c}$ with the discrete velocity $\overrightarrow{e_i}$ and assigning the discrete velocity distribution function the corresponding weight $w_i$. Substituting Equation (29)(29) $\mathrm{\Omega }\left(\mathit{f}\right)=\frac{-1}{\tau }\left(\mathit{f}-f^{\mathit{eq}}\right)$(29) into (28(28) $\frac{\mathit{\partial f}}{\mathit{\partial t}}+\overrightarrow{c}\frac{\mathit{\partial f}}{\partial \overrightarrow{r}}+\overrightarrow{a}\frac{\mathit{\partial f}}{\partial \overrightarrow{c}}=\mathrm{\Omega }\left(\mathit{f}\right)$(28) ) and discretizing the equation gives:

(31) $\frac{\partial f_i}{\mathit{\partial t}}+{\overrightarrow{e}}_{\mathit{i}}\frac{\partial f_i}{\mathit{\partial }\overrightarrow{r}}=\frac{-1}{\tau }\left(f_i-f_i^{eq}\right)$(31)

Integrating Equation (31)(31) $\frac{\partial f_i}{\mathit{\partial t}}+{\overrightarrow{e}}_{\mathit{i}}\frac{\partial f_i}{\mathit{\partial }\overrightarrow{r}}=\frac{-1}{\tau }\left(f_i-f_i^{eq}\right)$(31) in a small time interval Δt gives the lattice Boltzman equation [115,116]:

(32) $f_i\left(\overrightarrow{r}+{\overrightarrow{e}}_{\mathit{i}}\mathrm{\Delta t},\mathit{t}+\mathrm{\Delta t}\right)-f_i\left(\overrightarrow{r},\mathit{t}\right)=\frac{-\mathrm{\Delta t}}{\tau }\left[f_i\left(\overrightarrow{r},\mathit{t}\right)-f_i^{eq}\left(\overrightarrow{r},\mathit{t}\right)\right]$(32)

The equation expresses that the particle $f_i\left(\overrightarrow{r},\mathit{t}\right)$ moves to the neighboring point $\overrightarrow{r}+{\overrightarrow{e}}_{\mathit{i}}\mathrm{\Delta t}$ at the next time step $\mathit{t}+\mathrm{\Delta }\mathit{t}$ at the velocity ${\overrightarrow{e}}_{\mathit{i}}$. In numerical calculations, the equilibrium distribution function is generally taken to be its second order Taylor expansion [114]:

(33) $f_i^{eq}(\overrightarrow{r},\mathit{t})=\mathit{\rho }{w}_i[1+\frac{\overrightarrow{e_i}\cdot \overrightarrow{u}}{c_s^2}+\frac{{(\overrightarrow{e_i}\cdot \overrightarrow{u})}^2}{2c_s^4}-\frac{\overrightarrow{u}\cdot \overrightarrow{u}}{2c_s^2}]$(33)

where, c_s is the sound speed at the lattice reference temperature.

The relaxation time and the equilibrium distribution function in the LBM determine the model of the solved problem, which can be applied to the solution of different problems by imposing different equilibrium distribution functions and source terms. Compared with the Enthalpy method, it is suitable for complex fluids due to its easy boundary treatment and parallelization of the computation. This method can naturally be used for the calculation of phase change heat transfer in LMs. Yang et al. [117] proposed a LBM with two distribution functions and multiple relaxation factors based on the Enthalpy-based LBM of Huang et al. [118] and used it to study the Rayleigh-Bénard convective solid–liquid phase change problem in LMs. The computational results are in good agreement with the Enthalpy–porosity Method and the speed of computation is much faster, indicating the high accuracy and computational efficiency of the method.

4. Applications of liquid metal phase change materials

4.1. Pure liquid metal phase change materials

As shown in Figure 5, LMPCMs have great potential in thermal management of various electronic devices such as consumer electronics, high power density LEDs, data centers, IGBT modules, high power lasers, power batteries, etc. A series of work by Ge and Liu et al. from 2012 to 2013 demonstrated the great potential of LM, represented by gallium, in the field of passive thermal management of mobile storage devices and smartphones [20,119,120]. As shown in Figure 6(a), Zhao et al. proved that even compared with stearic acid (SA)/copper foam PCC, the thermal control performance of pure LM has obvious advantages [121]. The experimental results of Xu et al. also show that the performance of the heat sink with pure gallium as PCM is better than that of paraffin/copper foam composite heat sink (Figure 6(b)) [122]. The average effective thermal conductivity of gallium heat sink reaches 71.97 W/(m·K). Given the same heat storage capacity, compared with the paraffin/copper foam composite heat sink, the volume and weight of the gallium heat sink are reduced by 67% and 18% respectively. As shown in Figure 6(c), Alipanah and Li numerically compared the cooling capacity of gallium and octadecane used for thermal management of lithium-ion battery packs [123]. The results show that gallium can lower the battery temperature, improve the temperature uniformity, and lengthen the limit discharge time by 4.7 times longer than that of octadecane. Fan et al. tested a kind of Pb-Sn-In-Bi alloy heat sink whose thermal performance is much better than that of organic PCM (Figure 6(d)) [124]. When using the same volume of 80 mL, the performance of LM is obviously better than that of organic PCM at various heating power (up to 105.3 W/cm²). During the heating period, the presence of LM can significantly prolong the effective protection time by nearly twice and reduce the maximum superheating temperature by up to 50°C. Omari et al. studied the vibration-enhanced direct contact heat exchange between hot water and a heat sink composed of gallium [125]. The temperature rise in the melting zone is effectively limited by the high latent heat and low melting temperature of gallium. The results show that vibration enhances the heat transfer performance of gallium.

Figure 5. Potential application scenarios of LMPCM in the field of thermal management of electronic equipment.

Figure 6. Properties of single liquid metal phase change materials (LMPCMs). (a) Performance comparison of LMPCMs and SA, as well as SA/copper foam phase change composite (PCC) [121] (b) performance comparison of LMPCMs and paraffin/copper PCC [122] (c) performance comparison of LMPCMs and octadecane [123] (d) performance comparison of LMPCMs and octadecanol [124].

In theoretical research, gallium has been widely used as the research object of all kinds of melting models [126–128], while Yang et al. have constructed a numerical model of phase change heat transfer from the point of view of phase change heat storage or thermal control [129]. The phase change processes of gallium and eicosane are simulated and compared, and the results show that gallium has excellent melting characteristics that the paraffin does not have, and the proposed dimensionless relation provides a valuable reference for the design of the thermal management system. The effects of inclination and heating mode on the melting process of gallium were studied by Rivero [130]. The numerical simulation results show that the greater the inclination of the container is, the faster the melting process is. Positive dip angle can improve the melting process, negative dip angle can restrain natural convection, both of them can change the hydrodynamics of gallium during the melting process. It is found that only LM can restrain the temperature rise of local heat source, because it has thermal stability in a wide temperature range, which is not available for organic PCM.

Compared with other PCMs, LM usually exhibits significant advantages, but they still have some shortcomings. For example, the thermal conductivity of LM is still relatively low compared with other high thermal conductivity metals, leaving a space for further improvement; LM has low specific heat capacity and phase change enthalpy per unit mass, so it is not suitable to be used in situations sensitive to weight. The high cost of LM limits its practical application in engineering. LM is corrosive to copper and aluminum, and it is also a kind of conductor. Once leakage occurs, it will do great harm to electronic equipment. In order to conquer the above defects of LMPCMs, the performance of LMPCMs is optimized by integrating LM and other materials into PCC, or by combining LM phase change heat transfer process with other heat transfer methods.

4.2. Liquid metal phase change composites

4.2.1. Weight and cost optimization of liquid metal phase change materials

In order to control the weight and material cost of the system, the common strategy is to mix LM with other low-density and low-cost PCMs. By rationally designing the mixing method, mixing ratio and melting point, the weight and cost of the system can be reduced while the thermal management performance can even be further enhanced.

In the study of Shu [131], Bi_32.55In₅₁Sn_16.5 powder and paraffin were mixed by ultrasonic concussion according to the mass ratio of 1:1, 2:1 and 3:1. The low melting point alloy not only increased the thermal conductivity of the material, but also contributed part of the heat storage capacity. The introduction of a large number of alloys with small phase change volume change (abnormal expansion) also weakens the phase change volume change of the material of whole composite. There is a huge density gap between LM and organic PCMs, so the problem of phase separation after direct mechanical mixing is inevitable. Using expanded graphite as medium, Yao et al. prepared octadecane/LM PCC whose properties can be adjusted in a large range as shown in Figure 7(a) [132]. Octadecane and gallium are stored in the microstructure of expanded graphite, which solves the problem that it is difficult to mix organic PCMs and LM directly. At the same time, the microstructure of expanded graphite also plays a role in storing nucleated embryos, and the inhibition rate of supercooling is more than 86.82%. The latent heat of the prepared PCCs can vary from 211 to 70.19 J/g, and the corresponding thermal conductivity can vary from 2.639 to 39.158 W/(m·K).

Figure 7. The effect of combination of LMPCMs with other PCM. (a) Performance range of gallium/octadecane PCC [132]. (b) Comparison of the average bottom surface temperature of paraffin heat sink (PHS), wood alloy heat sink (WHS) and dual-PCM heat sink (DHS) [133]. (c) Temperature variation of solid gallium and PCM [134]. (d) Temperature rising under different gallium filling ratios in the 2-stage phase change heat sink [135].

McCann et al. stirred paraffin phase change capsules into Field’s metal [136]. Although the long-term stability of this direct mixing method remains to be investigated, the mixed use of organic and inorganic PCMs can bring more flexibility to the design of heat storage or thermal management systems. In order to reduce the weight and material cost of LM phase change heat sink, as shown in Figure 7 (b), Huang et al. proposed and analyzed a dual PCM heat sink based on numerical simulation, which can even achieve the effect close to that of pure LM heat sink by adjusting the ratio of LM and paraffin in the container [133].

Omari et al. proposed a new strategy to improve the applicability of gallium as a phase change heat sink [134]. Some cavities are prefabricated in the solid gallium and filled with organic PCMs with a melting point slightly lower than gallium. As shown in Figure 7 (c), when the hot fluid comes into contact with gallium, the gallium surface becomes molten. Although the internal temperature of gallium does not reach the melting point, it has reached the melting point of organic PCMs. Organic PCMs begin to absorb heat to slow down the growth of molten gallium and restrain the overheating of the melting layer with low heat capacity. As a result, the deterioration of the endothermic condition of gallium is prevented.

Ding et al. designed a two-stage PCM heat sink to alleviate the problem of photovoltaic cell overheating [135]. The pure gallium and RT25 (organic PCM, melting point 26.58°C) in the heat sink are isolated by aluminum plate, and the gallium is on the side near the heat source. The numerical analysis results show that the temperature control time of the PCC is the longest when the proportion of gallium is 56.9%, which is 20.2% and 78.8% longer than that of pure gallium and pure RT25 (Figure 7(d)). This result is of great significance for reducing the cost of PCMs and lightening the heat sink.

4.2.2. Heat conduction enhancement of liquid metal phase change materials

In the research of enhancing the thermal conductivity of LMPCMs, the most common reinforcement materials are built-in fins, carbon-based and metal porous materials.

As shown in Figure 8(a), Yang et al. demonstrated through numerical simulation that gallium with built-in copper fins as phase change heat sink can cope with thermal shock (100 W/cm²) which cannot be effectively handled by traditional PCMs [11]. Mozafari et al. designed a built-in crossed copper fin heat sink based on LM. As exhibited in Figure 8(b) [24], the effects of fin thickness, height and PCM thickness on the cooling performance were studied, and the numerical model was verified by experiments. Five cycles with a period of 45 s are monitored, and the results show that the optimal design can successfully control the bottom temperature below 48°C. Omari et al. integrated high thermal conductivity fins into LM heat sink [137]. The experimental results show that the fin can not only enhance the heat conduction capacity, but also enhance the natural convection in the phase change region and effectively improve the heat dissipation performance. Especially at high heat flux, the fins play a decisive role.

Figure 8. Effect of fins on performance of LMPCM heat sink. (a) The effect of the size of plate fins, cross fins and pin fins on the heat sink performance [11]. (b) The effect of the size of the cross fins on the heat sink desirability [24].

Zhao et al. infiltrated low melting point alloy into poco foam (a kind of carbon foam) by vacuum heating, improving the thermal conductivity of the composite significantly [25]. When the phase change occurred, the slope of the surface temperature response curve of the heater was greatly reduced, and the temperature difference was reduced from 21K to 3K (Figure 9(a)). Huang et al. compounded Wood’s metal into expanded graphite by heating and stirring method to form a form-stable PCM with high thermal conductivity [138]. And the distribution of wood alloy in the expanded graphite skeleton is uniform. As shown in Figure 9(b), the measurement results of thermophysical properties show that increasing the compaction density can reduce the internal voids of the composites, thus improving the heat transfer performance of the composites. At the same compaction density, increasing the content of expanded graphite can increase the thermal conductivity, but it will also decrease the latent heat.

Figure 9. Effect of thermal conductivity enhancing porous matrix on the performance of LMPCM heat sinks. (a) Temperature variations with time for different composite PCM at the heat flux of 7000 W/m² [25]. (b) Thermal conductivity and latent heat storage density of the composite material with different mass percentage of Wood’s alloy and compacting density [138]. (c) IR images of the thermal buffer [4]. (d) Top-cooled and bottom-cooled GaN device junction temperatures [22]. (e) Velocity distribution of copper foam/47 alloy and temperature curves of different PCM [139].

As exhibited in Figure 9(c), Yang et al. use a tube furnace and make the Field’s metal melt and infiltrate into the copper foam in a H2/Ar environment to form a PCC thermal buffer [4]. At the same time, the finite element method simulation was used to verify the test results and the effects of different thermal buffer thickness and copper volume fraction on the heat transfer performance were studied. After that, the PCC was integrated into the heat sink of the GaN device, and compared with the copper fin array, and the composite thermal buffer with higher cooling capacity can reduce the junction temperature of the device by 10–20°C (Figure 9(d)) [22]. Hou et al. used numerical methods to analyze the thermal control performance of 47 alloy (melting point: 51°C)/copper foam PCC [139]. It is found that copper foam weakens the convective heat transfer during the phase transformation process, effectively improves the response speed of the alloy, alleviates the phase transformation temperature rise (Figure 9(e)). And the higher the porosity of copper foam, the better the thermal management effect. However, the introduction of copper foam reduces the proportion of the alloy, so when the thermal control temperature limit is more than 80°C, the advantage of temperature control time is not obvious.

4.2.3. Encapsulation of liquid metal phase change materials

In the above research on the use of porous materials to enhance the thermal conductivity of LM, they can not only enhance the overall thermal conductivity of materials, but also play a certain role in packaging liquid PCMs. however, the leakage of LM cannot be completely avoided in this way. In order to avoid leakage, it is an important strategy to encapsulate LM into phase change capsules.

Yu et al. [26] mixed micron EBiInSn particles and expanded graphite into silicone rubber to make a flexible shell material with high thermal conductivity and phase change heat storage capacity, and phase change macro capsule was made by pouring with octadecyl alcohol/expanded graphite composite core. Compared with the shell material mixed with copper powder, researchers proved that EBiInSn particles can also significantly improve the thermal conductivity of the shell, and can provide heat storage function that copper powder cannot achieve (Figure 10(a)). As shown in Figure 10(b), Yao et al. have prepared a potential LM phase change capsule wrapped with self-assembled copper film [140]. The capsule not only has considerable energy storage density, but also can withstand the stress impact caused by the volume change of LM core in the phase change cycle. Raj et al. added 5 wt% nano-encapsulated liquid eutectic Ga-In alloy exhibited in Figure 10(c) to organic solid–solid PCM (SS-PCM) [141]. Because the LM has been packed in nano-capsules and limited in SS-PCM, the problems of agglomeration and phase separation of LMs have been solved. As shown in Figure 10(d), Kashiyama et al. prepared a kind of low temperature phase change capsule by chemical and oxidation treatment of dispersed gallium droplets [142]. The capsule is completely covered by the gallium oxide shell and the latent heat reaches 56% of pure gallium. Wang et al. prepared phase change microcapsule composed of LM core and PMMA shell (Figure 10(e)), and then mixed microcapsules with thermal gel to be used in the thermal management of LED chips [143]. The experimental results show that the composite with LM microcapsules can control the temperature of the chip at 71.1°C, which is 11.2°C lower than using single thermal gel. As exhibited in Figure 10(f), Gao et al. made phase change macro-capsule by wrapping 3 mm Field’s metal spheres with shell consist of Field’s metal and silicone rubber composites [144]. When the temperature reaches 110°C, the shape of the capsule can remain stable. The numerical simulation results show that the packed bed heat sink made of the phase change macro-capsules can reduce the temperature rise of pulse heat source by 73.4% and exhibits excellent thermal management performance.

Figure 10. Encapsulation of LMPCMs (a) the schematic diagram of preparation process of the spherical PCM microcapsule and thermal conductivities of the silicone composites with different mass fraction of EBiInSn [26]. (b) Schematic diagram of self-assembly mechanisms with images showing the state transition and mechanical testing and DSC measurement results of EGaIn@Cu. Reprinted with permission from ref [140]. Copyright 2021 American Chemical Society. (c) Nano-encapsulation of liquid gallium-indium eutectic metal alloy with in-situ polymerization method [141]. (d) Experimental procedure for Ga-based micro capsules synthesis [142]. (e) Forming schematic diagram of LM micro capsules [143]. (f) The preparation process of the spherical LMPM-PCM macro-capsules and the effect of spherical LMPM-PCM macrocapsules on the temperature response of convective cooling heat sink [144].

4.3. Combined phase change heat transfer of liquid metal

Although phase change thermal management technology has a series of advantages such as good temperature control effect and high temperature uniformity, its duration is limited, especially when dealing with long-term heating devices. PCM may not be able to recover all latent heat during solidification for the next working cycle, or even cause thermal failure after complete melting [145–147]. Because a single thermal management strategy is difficult to meet a variety of thermal management requirements [148], researchers have proposed hybrid cooling technology based on PCMs. Liu et al. designed and optimized a hybrid cooling system combining PCM, liquid cooling and fins, and the heat dissipation performance of the battery pack was studied by numerical simulation [149]. Xie et al. proposed a new hybrid thermal management system combining PCM, micro heat pipe array and liquid cooling, and the cooling performance of this system was verified by experiments [150]. In the field of LMPCMs, in order to meet the increasingly stringent heat dissipation requirements, Yang et al. put forward the concept of LM combined heat transfer, focusing on the composite cooling technology based on LMPCMs, and summarized the series and parallel heat transfer combined modes based on LMPCMs (Figure 11(a-c)) [151]. Yang et al. have proposed a quiet heat sink that combines LMPCMs with natural air convection, and established a transient heat transfer model [154]. The results of numerical simulation show that the combined heat sink can effectively control the heat source temperature when dealing with a 40W thermal shock lasting 100 s. By building heat pipes and fins inside the Field’s metal, with the aid of air cooling, Yang et al. have also realized a kind of heat sink which can cope with 1000 W periodic thermal shock produced by large electronic equipment (Figure 11(d)) [152]. As shown in Figure 11(e), Li et al. combined gallium with flat heat pipe to design a phase change thermal control module for electronic equipment with high packaging density, which can maintain the highest surface temperature below 45°C for several times longer than that of a single copper plate [153]. Even if the PCM is completely melted under the condition of transient heating, the module can keep the surface temperature of electronic equipment below 45°C in 600s, which indicates that the combined cooling technology based on LMPCMs is a feasible solution for thermal management of high-power portable electronic devices.

Figure 11. LMPCM combined series and parallel cooling system and application. (a) LMPM PCM as primary heat sink in series cooling system [151]. (b) LMPM PCM as secondary heat sink in series cooling system [151]. (c) Parallel cooling system [151]. (d) Thermal responses of the low melting point metal/finned heat pipe module under 200 W heating power [152]. (e) Module temperatures as a function of time at 3.65 W [153].

In addition to combining LMPCMs with other efficient heat transfer methods, there is also an application form of using carrier fluid to pump PCMs circularly between cold and heat sources, namely microencapsulated phase change slurry. Microencapsulated phase change slurry is a kind of fluid for convective heat transfer formed by dispersing phase change capsules into the liquid carrier, which is beneficial to improve the latent heat utilization efficiency of PCMs, as well as the effective heat carrying capacity of the fluid. LM phase change capsules have been used to improve the latent heat and thermal conductivity of the heat transfer fluid. As shown in Figure 12(a), Wu et al. dispersed the nano-Field’s metal particles treated with silane monolayers in HFE7100 dielectric fluid to form a stable suspension [155]. The suspension is used in jet impingement heat transfer test by absorbing heat through solid–liquid phase change of alloy particles and evaporation of dielectric liquid. In the range of alloy melting point, the heat transfer coefficient of suspension containing 30% alloy particles is 70% higher than that of pure HFE7100. Huang et al. proposed a method to prepare native or surface modified Field’s metal nanoparticles in one step based on the nano-scale emulsification process, and dispersed the nanoparticles in PAO to form a thermal functional fluid with good dispersion stability and thermal properties (Figure 12(b)) [156]. As shown in Figure 12(c), Mingear et al. prepared nano-sized gallium-indium alloy particles by ultrasonication, and the particles were wrapped by oxide film to form a capsule structure [27]. These LM phase change capsules were dispersed in AR-200 silicone oil to make a low-temperature microencapsulated phase change slurry. When the volume fraction of the capsules was 10%, the thermal conductivity of the phase change slurry at 24°C was 50% higher than that of pure silicone oil. Because the volumetric specific heat capacity of AR-200 and Ga-In capsules was similar, the volumetric specific heat capacity of the phase change slurry did not decrease after adding the microcapsule.

Figure 12. Liquid metal phase change slurry. (a) Schematic of synthesis process of nanoparticles and HFE7100 jet impingement (56–68°C) heat flux at the flow rate of 7.15 cc/s [155]. (b) DSC curves of bulk Field’s alloy pellets and corresponding nanoparticles with or without surfactant of ethyl carbamate synthesized at 180°C for 3 h [156]. (c) Ga–in nanoparticles with a schematic representation of a core/shell particle fabricated through ultrasonication process and the cumulative stored thermal energy starting at −40°C, including the increased energy from phase change [27].

5. Conclusion and perspective

In this review, the properties, heat transfer models and application research progress of LMPCMs are summarized. The main conclusions and future research perspectives are outlined as follows.

Compared to other PCMs, LMs have significant advantages in properties such as thermal conductivity and latent heat per unit volume, making them more suitable for working in the compact spaces of electronic devices. According to the results of fundamental studies on the thermophysical properties of metallic PCMs, the properties of LMs are related to their composition. For example, bismuth-based alloys have higher melting point and density, while gallium-based alloys have higher thermal conductivity and latent heat. The corresponding properties of LMPCMs can be influenced by adjusting the content of components with properties such as high melting points and high coefficients of expansion. The melting point of LMPCM is also pressure dependent. When the size of LMPCM is reduced to the nanometer scale, the melting point also decreases significantly. Microstructural defects in LMPCMs can also adversely affect their solid-state thermal conductivity.

Gallium-based LMs have significant supercooling effects. According to the nucleation theory, the size effect, heterogeneous nucleation, and thermal history all affect the supercooling degree of LM, so it is necessary to reduce the influence of these effects when studying the supercooling of LM. Although the crystal species suppresses supercooling better than the nucleating agent, a nucleating agent that needs to be added only once is clearly more suitable for practical applications for PCM heat sinks that are repeatedly frozen and thawed. In addition, the application of external fields such as magnetic fields, electric fields and mechanical vibrations can also partially inhibit supercooling in LMs.

In order to address the challenges associated with the coupling of phase change and natural convection as well as the motion of the phase interface in the classical Stefan’s problem, various numerical simulation methods have been developed. Among them, the Enthalpy–porosity Method based on continuous model solution and the LBM based on molecular dynamics have significant advantages in solving the above problem. The enthalpy method can effectively avoid the complex iterative calculations associated with tracking the moving interfaces, while the lattice Boltzmann method is particularly suitable for complex fluid systems similar to LMs due to its easy boundary treatment, computational parallelization, and higher computational efficiency.

LMPCMs offer significant advantages in thermal control over other PCMs. However, LMPCMs also have shortcomings, such as lower latent heat per unit mass and higher price. To overcome the shortcomings of LMPCMs and further expand their thermal management capabilities, researchers have developed two routes: liquid metal composite phase change materials and combined heat transfer. LM PCCs have been proposed to combine the benefits of multiple materials. In order to overcome the low latent heat per unit mass and high cost of LMPCMs, researchers combined them with low-density and inexpensive PCMs. If properly designed, the performance of heat sink with cascade design is even better than that of pure LM. By combining well-designed fins and carbon-based or metal porous matrices, the effective thermal conductivity of LMPCMs can be further improved. The packaging of LMs in the form of phase change capsules can effectively prevent the risk of corrosion and short circuit caused by the leakage of molten LMPCM.

Under stronger thermal shock, the LMPCMs will fail if the absorbed latent heat cannot be released in time. In order to solve this problem, researchers have proposed the LM combined heat transfer technology, which combines LM phase change heat sinks with other efficient heat transfer or heat dissipation modalities. LM phase change microcapsules can also be directly dispersed into fluids to participate in the convective heat transfer process, absorbing heat and melting at the heat-generating device, while releasing heat and returning to the solid state at the heat dissipator, thus increasing the convective heat transfer coefficient. At the same time, LM microcapsules also help to improve the thermal conductivity and apparent specific heat of the carrier fluid.

As the research on LMPCMs continues, several scientific and technical challenges need to be further considered. Currently, there are few non-toxic low-melting-point alloys available for phase change thermal management of electronic devices that do not contain lead and cadmium. Therefore, there is a need to discover more LMPCMs with melting points in the range of 30 ~ 60°C. The combination of calculating material properties from first principle and predicting new materials using machine learning has emerged as a highly effective approach for discovering novel materials. This approach is expected to efficiently improve the material system of LM, establish a comprehensive database of the corresponding physicochemical properties, and ultimately construct a genome of LM materials. Current research on LMPCMs is still dominated by eutectic alloys, while non-eutectic LMs are preferred for applications requiring a wider range of phase change temperature control. This is because they have a broader phase change temperature region and can fill the gap of non-toxic eutectic alloys with melting points in the range of 30 ~ 60°C.

LMs are highly susceptible to supercooling, which in turn can greatly affect the thermal control performance of LMPCMs. Although there are various methods, such as adding nucleating agents or applying external fields, to prevent the supercooling of LMs, it is necessary to find more effective methods to prevent LM supercooling since LM are generally significantly supercooled, for example, the minimum supercooling degree of gallium is still between 10°C and 20°C.

Close contact melting has been proved to be helpful to improve the performance of LMPCMs [101]. However, there is still a lack of systematic and in-depth research on the close contact melting of LM. In addition, there is no practical design scheme of LM phase change heat sink based on close contact melting that can work continuously.

Most of the current studies based on LM PCC heat transfer have been combined with fins and heat pipes, and fewer studies have been combined with other cooling methods. In future research work, more LM composite heat transfer methods should be tried, such as phase change combined with convection cooling, phase change combined with thermal radiation. LM phase change-based composite heat transfer technology may become the best solution to face the increasingly complex heat dissipation requirements for thermal management of electronic devices.

LMs are good carriers of both heat and electricity. Due to their conductivity, the leakage of LMPCMs in electronic devices can lead to short circuits and consequent serious damage. In addition, a series of LMs represented by gallium and gallium-based alloys have different degrees of corrosiveness to copper, aluminum and other commonly used metal materials, and the corrosion phenomenon is more serious with the increase in temperature. In practice, it is necessary to consider the compatibility of LMPCMs with the substrate and take good anti-corrosion measures.

Due to their unique physicochemical properties, LMs not only have outstanding performance in thermal control, but also have many other characteristics such as flexibility, electromagnetic shielding and good biocompatibility. Therefore, it is possible to combine these characteristics to create multifunctional PCMs, such as LMPCMs with electromagnetic shielding function, which can be used in aerospace and medical fields where electromagnetic shielding is required for thermal control, and flexible PCMs, which can be used in the manufacture of thermally controlled flexible robotic skins. It is expected that LM phase change multifunctional materials will greatly broaden the related research and application fields in the coming time.

Author contribution

Jing Liu, and Zhongshan Deng conceived the idea. Yuchen Yao and Wei Li composed the first version of the manuscript. All authors participated in interpretation of the results and revised the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [51890893].