Comprehensive Overview of the Effective Thermal Conductivity for Hydride Materials: Experimental and Modeling Approaches

Abstract

In metal hydride beds (MHBs), reaction heat transfer often limits the dynamic performance. Heat transfer within the MHB usually involves solid and gas phases. To account for both, an effective thermal conductivity (ETC) is defined. Measuring and predicting the ETC of metal hydride beds is of primary importance when designing hydride-based systems for high dynamics. This review paper presents an integral overview of the experimental and modeling approaches to characterize the ETC in MHBs. The most relevant methods for measuring the ETC of metal hydride beds are described, and the results and scopes are shown. A comprehensive description of the models applied to calculate the ETC of the MHBs under different conditions is developed. Moreover, the effects of operation parameters such as P, T, and composition on the ETC of the presented models are analyzed. Finally, a summary and conclusions about experimental techniques, a historical overview with a classification of the ETC models, a discussion about the needed parameters, and a comparison between ETC experimental and calculated results are provided.

1. Introduction

Fossil fuels mainly cover the ever-increasing global energy demand nowadays, and cause global warming due to greenhouse gas (GHG) emissions [1]. Therefore, there is a need for sustainable and renewable energy sources, which has led to the research of clean and efficient energy carriers [2]. Hydrogen is one of the most promising options, thanks to its abundance and high gravimetric energy density (LGHV (low gravimetric heating value) of hydrogen: 120 MJ/kg or 33.33 kWh/kg, in comparison to the LGHV of methane: 55.6 MJ/kg or 15.44 kWh/kg) [3,4,5]. However, hydrogen is a flammable gas with a very low boiling point (around 20 K at 1 atm), making it hard to store and transport hydrogen efficiently and safely. In addition, hydrogen has a low volumetric energy density (LVHV (low volumetric heating value): 10.8 MJ/Nm³ or 3.0 kWh/Nm³) [3]. Physical storage methods, like high-pressure gaseous hydrogen storage (at 350 bar and 700 bar) or liquid and cryo-compressed hydrogen storage, are characterized by high costs for compression and/or liquefaction (5–20% and 30–40% of the LGHV hydrogen is required for compression and liquefaction, respectively) [4,5,6,7,8,9]. Hydrogen storage in carrier materials and, especially, in solid-state, in which hydrogen is bonded to a host material, is a promising alternative [10,11,12,13,14].

Thanks to the high hydrogen volumetric capacity ranging from 60 to 150 kg/m³, wide pressure and temperature operative ranges (cryogenic temperature to 500 °C, and 1 bar to more than 700 bar), and reversibility of the absorption–desorption process, metal hydrides (MHs) are among the materials with the highest potential for hydrogen storage in solid-state for mobile and stationary applications [15,16,17,18,19,20,21,22]. Furthermore, MHs are also used for many other applications, including heat pumps, thermal energy storage systems, heat transfer systems, and thermal hydrogen compressors [17,18,23,24,25,26].

The reversible hydrogenation/dehydrogenation process in MHs under equilibrium conditions can be described with the following Equation (1) [27]:

$$xMe+{\displaystyle \frac{y}{2}}{H}_2\rightleftharpoons {Me}_x{H}_y+{∆H}_r$$

(1)

where Me is a metal, an alloy, or a compound, Me_xH_y is the hydride compound, and ΔH is the reaction enthalpy.

The equilibrium pressure, P_eq, is described by the Pressure–Composition Isotherms (PCIs) in the alpha–beta (α_ss − β_MH) equilibrium condition, as can be seen in Figure 1A. An ideal Me/Me_xH_y system presents flat plateaus that provide the P_eq for the hydrogenation and dehydrogenation processes at different temperatures (Figure 1A). The P_eqs at different T_eqs, equilibrium temperatures, are used to calculate the enthalpy, ΔH_r, and entropy, ΔS_r, of the hydrogenation/dehydrogenation reaction through the van’t Hoff equation (Figure 1B). The ΔH_r represents the strength of the Me-H bond and is derived from the slope of the van’t Hoff plot. According to Equation (1), the hydrogen absorption reaction is an exothermic process, while the desorption is an endothermic process. The heat absorbed or released varies depending on the hydride. For example, the enthalpy of the reaction for the formation/decomposition of the hydride compound for LaNi₅ is approximately 30 kJ·mol H₂⁻¹, while for Mg, it is around 75 kJ·mol H₂⁻¹ [28,29,30]. The ΔS_r is taken as the hydrogen entropy variation from molecular gas to dissolved solid and vice versa, neglecting the entropy change in the α_ss and β_MH phases; the coordinate to the origin determines a ΔS_r value of about 130 J·K⁻¹·mol⁻¹, Figure 1B [30,31]. According to the van’t Hoff Equation, the equilibrium pressure varies exponentially with the temperature, and the equilibrium conditions of the system are defined just by a P_eq and T_eq condition, as shown in Figure 1C [30,31]. The hydride compound exists on the left side of the equilibrium line, and on the right side, the metal phase is present.

However, the real Me/Me_xH_y system exhibits plateaus with a slope and hysteresis between the hydrogenation and dehydrogenation plateau, as shown in Figure 2. The cause of sloped plateaus can be attributed to the H-H interaction, since the space lattice does not continuously increase with the hydrogen concentration [32]. The hysteresis phenomenon can be ascribed to the material’s strain, which generates a plastic deformation when it absorbs hydrogen [32]. Kinetic limitations near the equilibrium conditions can also lead to an apparent hysteresis in the experiments.

For hydrogen absorption and desorption under dynamic conditions, if the applied hydrogen pressure, P_abs/des, at a given T is higher/lower than the absorption/desorption equilibrium pressure P_eq,abs/des, the Me/Me_xH_y system starts to absorb/desorb hydrogen as described in Figure 2. The difference between P_eq,abs/des and P_abs/des is the driving force for the chemical reaction. Therefore, if the heat associated with the hydrogenation and dehydrogenation reaction is not efficiently taken out or provided from/into the Me/Me_xH_y system, the reaction rate slows down. Such slow rates are associated with the increase/decrease in the hydride bed temperature, which changes the P_eq and reduces the driving force for the hydrogenation and dehydrogenation processes.

In most cases, the hydride-forming materials are in powder form. The hydrogen gas surrounds the solid particles upon hydrogen interaction. The primary heat transport phenomenon upon hydrogen absorption and desorption is conduction. In the case of metal hydride beds (MHBs), the thermal conductivity involving the solid and gas phases is called effective thermal conductivity (ETC) [33,34,35,36,37,38]. The swelling phenomenon also influences the performance of metal hydrides. The hydride-forming material expands and shrinks during cycling due to the metal’s varying crystal parameters, and thus, pulverization occurs during repeated hydrogen absorption and desorption cycles [39,40], usually down to a certain particle size. With a reduction in the particle size, the ETC of MHBs drastically decreases. The ETC of hydrides is relatively low and in the range of ≈1 W m⁻¹ K⁻¹, and therefore limits the hydrogenation and dehydrogenation reaction rate [33,34,35,36,37,38,39,40]. Hence, measuring and predicting the ETC of metal hydride beds is of primary importance when designing hydride-based systems, particularly if high absorption or desorption dynamics are requested. Experimental methods and mathematical models are required to evaluate the ETC of MHBs. In both cases, however, some issues arise due to the complexity of MHBs, owing to their powdered nature, the swelling phenomenon, incomplete hydrogenation/dehydrogenation reactions, and, in some hydride systems, the formation of intermediate phases. Hence, effective thermal conductivity becomes a complex function of various variables over the hydrogenation and dehydrogenation reaction.

Regarding the experimental measurements, conventional methods and commercially available instruments can usually not fully cover the complex nature of an MHB and can only provide ETC values with certain experimental constraints. It is necessary to perform measurements in a hydrogen pressure-controlled atmosphere of porous media made of very small particles. Because of this, innovative approaches have been proposed over time, and several examples are reported in the literature [33,41]. However, challenges remain, and, so far, it is still impossible to identify a general experimental procedure and setup for measuring the effective thermal conductivity of metal hydride beds.

For modeling the ETC of material beds, various equations and methods have been proposed [42,43,44,45]. However, the main limitations related to this class of “general” models lie in the lack of account for some peculiar aspects of MHBs, like the reactivity between the solid and the gas phase, the consequent expansion/compression of solid particles, or the heat release from the reaction, as expressed in Equation (1). Some of these models, like the Zehner–Schlünder model, have been subsequently adapted to metal hydride systems, while others have been directly developed for this purpose [46,47,48,49,50].

This review paper presents an integral overview of the experimental and modeling approaches to characterize the ETC in MHBs. Section 2 highlights the most relevant methods for measuring the ETC of metal hydride beds, describing the results and limitations. Section 3 provides a comprehensive development of the models in the literature employed to calculate the ETC of the MHBs at different conditions. The difficulties in estimating every parameter are also discussed. Most of the parameters regarding the solid phase undergo a relevant variation during operation, as the solid matrix is originally constituted and is first activated and then cycled. Being interested in the hydrogenation–dehydrogenation processes, the parameters should be evaluated after cycling once the material can be considered stable. In addition, how the parameters change during ABS/DES cycles, including the expansion/compression related to the process, should also be considered. The thermal conductivity of the solid phase is also investigated, a critical parameter that has not been covered so far in the literature on the hydrides’ field. Section 4 is about the effects of pressure, temperature, and composition on the ETC of the models presented in Section 3. An analysis of the behavior of the models comparing experimental data for LaNi₅ is also performed. Finally, a summary and conclusions about experimental techniques, an overview with a classification into gas film and particle modification models of the ETC models developed so far, and a discussion about the needed parameters are provided.

2. On the Overview of the Methods Applied to Measure Thermal Conductivity

This section describes the principles of the experimental methods and setups for measuring the ETC of metal hydrides. The obtained results and constraints of the methods are discussed, covering the following main aspects:

→

Time required for the measurements;

→

Possibility to measure in a hydrogen atmosphere;

→

Possible temperature and pressure range;

→

Type of sample;

→

Amount of sample;

→

Commercial availability.

Figure 3 proposes a classification of the existing methods to determine the ETC experimentally. These methods are generally employed not only in metal hydrides but also in polymers, cement, and clays, among other materials. Such methods can be divided into two main categories: steady-state and transient [33,41,51,52,53,54,55]. The main difference here is how the heat source operates during the measurement. In the steady-state methods, the heat is provided from one side of the sample (heat source) and transferred towards the other side of the sample (heat sink). The measurement continues until a constant temperature gradient in either the axial or radial direction, depending on the cell design, is reached. The steady-state methods can be divided into two subgroups, depending on the direction of the heat, in axial and radial methods (Figure 3). Applying the steady-state method, the ETC can be calculated from the temperature values at different points of the sample, the constant power of the heat source, and the distances from the heat source to the temperature measurement points. The steady-state method can be applied directly, just with the sample, or in a comparative way, including a reference material. In this section, both subgroups of steady-state methods will be further analyzed.

In contrast, in the transient methods, the heat for the measurement is given as a single or periodic pulse. The ETC is calculated from the sample’s temperature change over time. In the case of transient methods, it is possible to distinguish two main techniques: the hot wire and the laser flash. While the two main categories of the ETC measurement methods are well defined, their distinction in sub-categories is less strict and depends on the research field and the author.

2.1. Steady-State Methods

2.1.1. Axial Heat Flow

For the most straightforward axial heat flow method, a sample of known dimensions (usually a cylinder or cuboid) is heated until steady-state conditions are reached. At this point, a one-dimensional heat flow through the sample can be assumed (with a constant axial temperature in every sample layer). The Fourier law, Equation (2), describes this case and has the advantage of easy data evaluation [56], as follows:

$$q=-kA{\displaystyle \frac{dT}{dx}}$$

(2)

where q is the heat flow (W), A represents the cross-section area of the sample (m²) through which the heat flows, dT/dx (K m⁻¹) is the temperature gradient, and k is the thermal conductivity of the sample (W m⁻¹ K⁻¹).

Even this straightforward setup has been used for hydride compounds. Eaton et al. measured the ETC of different LaNi₅-lead composites, reporting up to 2 W m⁻¹ K⁻¹ [57]. Gradually, several modifications to the basic setup have been made to reduce the error and to make the measurement more robust. For instance, the system is isolated to prevent heat loss in the radial direction and to ensure a one-dimensional heat flow through the sample. Another possible modification is adopting a comparative setup instead of an absolute one. In this case, the sample is attached to one or two materials of known thermal conductivity, which act as references. The approach is the so-called comparative cut-bar method, and the concept is shown in Figure 4 [58], where the sample is embedded between two references with known thermal conductivities, k_Ref1 and k_Ref2. Assuming a constant heat flow through the whole stack, the ETC can be calculated through Equation (3) by measuring the temperature T (K) at given positions Z (m) instead of evaluating the heat flow (e.g., by the heat source’s energy consumption), as follows:

$$k_{eff}={\displaystyle \frac{Z_4-Z_3}{T_4-T_3}}·\left({\displaystyle \frac{k_{Ref1}}{2}}·{\displaystyle \frac{T_2-T_1}{Z_2-Z_1}}+{\displaystyle \frac{k_{Ref2}}{2}}·{\displaystyle \frac{T_6-T_5}{Z_6-Z_5}}\right)$$

(3)

Another possible application of the axial steady-state method is the guarded hot plate method, where the measurement is isolated and equipped with a so-called guard heater to compensate for possible heat losses. The guard is heated at the same temperature as the measurement’s heat source, so there is no radial driving force for the heat. The heat can only travel through the sample in the axial direction. Another modification is to calculate the heat flow through the setup directly with a heat flux sensor instead of calculating it from the energy consumption, leading to the guarded heat flow meter method. Additionally, there might be some slight modifications, like a setup where the heater is placed between two identical samples, which would not change the overall working principle [59].

Due to certain experimental limitations (e.g., the need for good contact between the sample and the reference material) [58], the axial methods are more suited for samples with a flat surface, like metal hydride powders compacted into pellets. For instance, Sanchez et al. [60] have used the comparative cut-bar method to measure the ETC of LaNi_4.85Sn_0.15 containing expanded natural graphite (ENG) perpendicular to the pressing direction, improving the thermal conductivity up to ≈19 W m⁻¹ K⁻¹. Similarly, Park et al. [61] have used the same comparative cut-bar method to investigate the ETC of La_0.9Ce_0.1Ni₅ pellets with ENG, reaching a thermal conductivity of 8 W m⁻¹ K⁻¹ (not specifying the measurement direction). Yasuda et al. [62] have used the heat flow meter method to measure the heat conductivity of two different metal hydrides (TiFe_0.9Ni_0.1 and La_0.6Y_0.4Ni_4.9Al_0.1) pressed with aramid and carbon fiber into sheets. They measured an ETC up to 3 W m⁻¹ K⁻¹ in the planar direction and nearly no heat conductivity through the plane. Kumar et al. [63] applied the guarded hot-plate method to measure the ETC of MmNi_4.5Al_0.5. They placed the sample in the measurement chamber and slightly compressed it, forming a disk. Measurements have been carried out in a hydrogen atmosphere at 50 bar and about 100 °C, with an ETC value of about 1.2 W m⁻¹ K⁻¹.

The axial steady-state method’s general advantages are the measurements’ high accuracy and the straightforwardness of data acquisition. The guarded hot plate method is one of the most precise methods for measuring heat conductivity [41]. However, to reach this high precision, the cell must be well designed and built, and the measurement devices, such as thermocouples and heat sources, must have a high accuracy. This precision comes at the price of long measurement times, which can be hours or days until steady-state conditions are reached [52]. While there are commercial setups whose measurement principle is based on the axial steady-state method [64], their application is aimed towards insulating materials, and often large sample sizes are required. Hence, the reported results in the scope of metal hydrides are made with custom-made devices, and consequently, the operation conditions (sample size, hydrogen atmosphere, and operative temperature) depend on how the device was constructed.

2.1.2. Radial Heat Flow

In radial methods, the sample is placed inside a circular sample holder with a certain number of thermocouples at varying distances from the heat source. In most setups, the heat flow is directed from the center of the sample holder towards the shell, which acts as a heat sink. Similarly to axial methods, radial methods can operate in absolute and comparative modes. Figure 5A shows a schematic design of the absolute radial method: r₁ and r₂ (m) are the radial positions of two thermocouples, T₁ and T₂ (K) are the corresponding measured temperatures, and h_sh (m) is the sample holder radius. In general, radial heat flow setups utilize more than two thermocouples to improve the precision of the measurement. There are reports of setups using three [65], four [66], or even five thermocouples [67], as seen in Figure 5C. Using the schematic setup of Figure 5A, quantifying the heat generated q (W), the ETC can be calculated using Equation (4):

$$k_{eff}={\displaystyle \frac{1}{T_2-T_1}}·{\displaystyle \frac{q}{2\pi h_{sh}}}\ln \left({\displaystyle \frac{r_2}{r_1}}\right)$$

(4)

In the comparative mode, the sample holder contains the sample and a reference material. The reference should have a low thermal conductivity to establish a high-temperature gradient and reduce the measurement error; typical reference materials are PTFE [47] or Nylon [68]. In Figure 5B, the setup contains four thermocouples in total (two in the sample and two in the reference material) with distances from the heat source of r₁ to r₄ (m) and measured temperatures of T₁ to T₄ (K).

In contrast to the absolute method, in the comparative one, it is not necessary to know the heat generated by the heat source, and with the knowledge of the heat conductivity for the reference material, k_Ref (W m⁻¹ K⁻¹), the ETC for the sample can be calculated according to Equation (5) [33].

$$k_{eff}=k_{Ref}{\displaystyle \frac{T_4-T_3}{T_2-T_1}}·{\displaystyle \frac{\ln \left(\frac{r_2}{r_1}\right)}{\ln \left(\frac{r_4}{r_3}\right)}}$$

(5)

Compared to the axial steady-state method, the radial method is more robust (nearly no heat loss due to the placement of the heat source in the middle of the setup), and less or no isolation is needed to prevent undesired heat losses. The outer shell often acts as a heat sink; in this case, isolation is even undesirable. This makes the radial steady-state method suited not only for samples with flat surfaces, like pellets, but also for powder samples. The comparative method might have a slight advantage because the heat source is responsible for the most significant error in the calculation [60], and a lower amount of sample might be needed compared to the absolute method. Nevertheless, both methods have been used to determine the ETC (k_eff) of hydride compounds with no distinct preference. Due to the flexibility of the radial method mentioned above, it has been used for very different samples. Shim et al. [66] used the radial flow method to measure the ETC of MgH₂ pellets with added NbF₅ and ENG in the absolute mode, using a vessel with four inserted thermocouples, as shown in Figure 5C. They have reached a value up to 13 W m⁻¹ K⁻¹ for samples containing 20 wt% ENG with a particle size of 200 µm and measured up to 70 bar at RT.

In contrast, Wang et al. [68] have measured the ETC of a LaNi₅ bed in the comparative mode, reaching a value of 2.1 W m⁻¹ K⁻¹ at 16 bar of hydrogen pressure. Recently, X. Mou et al. [69] have reported ETC measured for LaNi₅ with different initial particle sizes and porosities under a He atmosphere and an average temperature of 20 °C in an in-house design cell. It was found that a decrease in porosity and an increase in particle size could raise the ETC. Additionally, measurements of the LaNi₅ powder beds were performed under 1 to 40 bar of hydrogen pressure and an average temperature between 20 to 60 °C, providing values from 0.1 to 1.7 W m⁻¹ K⁻¹.

A general disadvantage of the radial steady-state method is that the cells need to be large enough to have the proper space to allocate several thermocouples. Consequently, a large sample is needed. Madeira et al. [70] used 310 g of La_0.8Ce_0.2Ni₅, while Yang et al. [65] needed 290 g of LaNi₅. The radial steady-state method shares the same disadvantage of long measurement times as the axial steady-state method until the steady-state is reached. Furthermore, to our knowledge, there are no available commercial setups for the axial steady-state method, and all reported results have been measured using custom-made devices. The setup design constrains the experimental temperature range, hydrogen pressure, and sample amount among the most relevant parameters.

2.2. Transient Methods

2.2.1. Hot-Wire Method

In contrast to steady-state methods, the heat for ETC measurement is provided as a pulse for transient methods. Hence, no waiting time to reach a steady state is needed, significantly reducing the measurement time from hours to minutes. Figure 6 shows the general working principle of the hot-wire method. It is based on using a wire placed inside the sample material. The wire should be infinitely long and thin and act both as a point heat source and sensor. The electric impulse heats the wire by Joule heating, creating a heat flow towards the surrounding media. Samples with a higher thermal conductivity can dissipate this heat faster, and the wire remains colder than a material with a low thermal conductivity. At a distance from the hot wire, zones with similar temperatures (isothermal lines) exist. Four different sources/sensors can be utilized for such transient measurements: the transient line source (TLS), transient plane source (TPS), modified transient plane source (MTPS), and the eponymous hot wire. They all share a similar working principle, highlighted by the fact that, recently, devices that can operate with all four source types have been developed [71].

The ETC of the sample can then be calculated according to Equation (6), where k_eff is the ETC of the sample W m⁻¹ K⁻¹, q the heating power (W), t₁ and t₂ (s) the times at which the measurement started and ended, respectively, and ΔT (K) the resulting temperature difference [72].

$$k_{eff}={\displaystyle \frac{q}{4\pi }}·{\displaystyle \frac{\ln t_2-\ln t_1}{∆T}}$$

(6)

The temperature difference is mainly measured by the changes in resistance in the wire. The wire should have a linear dependency between the resistance and temperature in an extensive temperature range; therefore, platinum wires are mainly used. However, some devices use two wires with different lengths and orientations (cross-wire and parallel-wire techniques) [72].

Sundqvist et al. [73] used the hot-wire methods to investigate potential complex hydride materials (alanates and borohydrides) and performed measurements between 100 K–300 K and up to a 2 GPa hydrogen pressure. For instance, they measured NaAlH₄ at RT and 2 GPa for a thermal conductivity of about 3.7 W m⁻¹ K⁻¹, while NaBH₄ showed a lower thermal conductivity of about 1.5 W m⁻¹ K⁻¹ at similar conditions. One of the main problems of the hot-wire method is that the wire must be as thin as possible, leading to very delicate systems. The TLS (transient line source or needle probe) overcomes this problem. In the TLS, the wire is encased into a metal needle. This needle contains two thermocouples: one at the middle of the needle (where the highest temperature is expected to be) and one at the tip of the needle as a reference point. The TLS method is more robust than the hot-wire method, sacrificing some accuracy. The TLS is not as thin as the hot wire, which leads to additional axial heat losses [51]. Dedrick et al. [74] used the TLS method to measure the heat conductivity of sodium alanate between RT and 130 °C and up to 100 bar, reaching about 0.5 W m⁻¹ K⁻¹. While the hot wire and TLS have the general advantage of fast measurement times and can be used at temperatures up to 200 °C (which might not be enough for high-temperature hydrides like MgH₂ [75]), they have to be embedded into the sample and can therefore be only used for powder beds. Furthermore, there are no commercially available systems for measurements in a hydrogen atmosphere, and the reported results have been obtained from custom-made pressure cells.

In contrast to the hot-wire and TLS methods, the TPS method utilizes a metal spiral (usually made of nickel) embedded into a temperature-stable polymer, like Kapton or Mica, which also acts as an insulator. For the measurement, different parameters (e.g., the measurement time or the electric power) have to be adjusted to only slightly heat the sample (about 1–2 K), so that the correlation between the resistance and the temperature increase can be assumed as linear [76]. The surrounding material’s thermal conductivity can be calculated depending on the sensor’s temperature change. The standard TPS sensor is planar, and the heat pulse is emitted in the axial direction. Therefore, two similar samples of sufficient thickness and diameter are needed (Figure 7A) so that the heat pulse of the sensor remains in the sample and there are no influences from the surrounding media. However, the two samples do not need to be identical. They only have to be large enough (twice the radius of the TPS and at least as thick as the radius of the TPS) [33,76], and some deviation between the samples can be tolerated. Hence, the requirements for the samples are not as strict as those for the steady-state methods. In this regard, the MPTS method has the advantage that the working principle and application are similar to the standard TPS method but only needs one sample (Figure 7B).

The flexibility of the TPS method is further highlighted by the possibility of measuring both powder [77,78,79] and pellet [80,81,82] samples. This flexibility might contribute to the observation that the use of the TPS method has increased recently, while there are nearly no recent publications on the hot wire with TLS methods for hydride compounds.

Similar to the hot-wire and TLS methods, the TPS sensors can work in an extensive temperature range (up to 300 °C [71]), mainly limited by the stability of the polymer embedding. Another advantage is the fast measurement time. If the starting parameters (electric energy, measurement time, etc.) are set correctly, a measurement can be conducted in minutes. Again, a challenge in metal hydride research is to measure thermal conductivity by TPS in a hydrogen atmosphere. For instance, Jepsen et al. [77] measured the thermal conductivity of Li-Mg reactive hydride composites up to 180 °C in an argon atmosphere inside a glove box, potentially even underestimating the effective thermal conductivity as hydrogen is more conductive than argon [83]. Other authors have even measured the thermal conductivity of metal hydrides without a protective atmosphere: Atalmis et al. [82] have measured the thermal conductivity of LaNi₅ pellets coated with Cu and ENG, reaching up to 13.82 W m⁻¹ K⁻¹ for samples containing 20 wt% ENG at RT. However, efforts have been made to measure the thermal conductivity of metal hydrides directly in hydrogen atmospheres with a TPS sensor system. For instance, Albert et al. [78] constructed a pressure vessel with an embedded TPS sensor so that the thermal conductivity of metal hydrides can be measured directly in a hydrogen atmosphere and investigated a Ni-MgH₂ blend up to 25 bar of hydrogen pressure at 300 °C, reaching up to 0.9 W m⁻¹ K⁻¹.

2.2.2. Laser-Flash Method

In the laser-flash method (LFM), the sample is heated by a laser pulse. The heated sample radiates IR light, which is caught by an IR detector (Figure 8). In order to improve the absorption of the laser radiation and to level surface roughness to a certain degree, the sample should be coated with a layer of carbon before the measurement is performed [84]. The thickness of the coating must be negligible in comparison to the sample thickness, so it does not influence the calculation of the thermal conductivity.

In contrast to other thermal conductivity methods, the LFM does not measure the thermal conductivity directly; instead, it measures the sample’s thermal diffusivity, α_d (m²/s) (Figure 9). However, with knowledge of the density ρ (kg/m³) and the heat capacity c_p (J kg⁻¹ K⁻¹] of the sample, the ETC k_eff (W m⁻¹ K⁻¹) can be calculated according to Equation (7):

$$k_{eff}={\alpha }_d·\rho ·c_p$$

(7)

Due to the measurement principle of the LFM, it is challenging to apply it to metal hydride beds in powder form, so, in the scope of metal hydride research, this method is mainly used for pellets. For samples with low thermal conductivity, the gas trapped in the pores, the convection mechanism, and the different path lengths for the radiation can lead to biased results [86]. While the LFM can be used in an extensive temperature range (up to 3000 °C [59]), the measurement in a reactive atmosphere, like hydrogen, is problematic and might even harm the LFM system [86]. The measurement must be conducted in an inert gas flow (nitrogen or helium) or at a low hydrogen pressure. Pohlmann et al. [87] used the LFM to measure the heat conductivity of Hydralloy C5 (a commercially available AB₂-alloy) pellets with varying amounts of ENG in a custom-made protective cell for air-sensitive samples [88], reaching up to 63 W m⁻¹ K⁻¹ for a sample containing 12.5 wt% ENG orthogonal to the pressing direction. They could further show that the thermal conductivity depends mainly on the ENG content, while the compaction pressure for the pellets has virtually no impact. Popilevsky et al. [89] investigated MgH₂ pellets with embedded 2 wt% multiwall carbon nanotubes and measured the ETC between RT and 250 °C, reaching up to 35 W m⁻¹ K⁻¹.

The LFM has gained popularity for hydride compounds mainly due to fast measurement times, which are in the order of seconds and are even faster than the hot-wire methods. With known measurement parameters, the samples can be placed in the commercially available system, and the measurement can be started. This might be further accelerated by the possible addition of automatic sample changers so that several samples can be measured in a row. In the hot-wire methods, the sensor has to be removed from the sample and placed in a new one before the measurement can be repeated, adding downtime. However, finding the correct set of measurement parameters for the LFM, e.g., the proper voltage for the laser source, is not trivial. Furthermore, optimizing the initial data is not fully automated, and the user must input the measurement’s starting values [90]. Therefore, the measurement’s quality can significantly decrease if wrong initial values are chosen, leading to considerably longer measurement times.

3. Development of Models for Effective Thermal Conductivity

This section comprehensively describes the models used to represent the k_eff in metal hydride beds. A total of 13 models developed from 1873 to 2023 are discussed. Moreover, the theoretical background of the thermal conductivity of the metal–metal hydride as a solid phase is exposed, and the effect of hydrogen in the metal lattice on the solid thermal conductivity is analyzed.

3.1. Description of the Models and Parameters to Calculate the ETC (k_eff)

3.1.1. Maxwell Model (1873)

To the authors′ knowledge, the solution proposed by Maxwell [42] can be considered the oldest model for evaluating the effective thermal conductivity of a composite material made by solid particles dispersed in a fluid phase. The model applies to a three-dimensional heat flow for low concentrations of the particulate [48]. Equation (8) describes the model:

$$k_{eff}=k_g{\displaystyle \frac{1+2\frac{\frac{k_s}{k_g}-1}{\frac{k_s}{k_g}+2}\left(1-\epsilon \right)}{1-\frac{\frac{k_s}{k_g}-1}{\frac{k_s}{k_g}+2}\left(1-\epsilon \right)}}$$

(8)

where k_eff is the ETC, k_g is the thermal conductivity of the gas (or fluid) phase, k_s is the solid thermal conductivity, and ε is the porosity. This model appears very simple and is characterized by only three parameters: k_s, k_g, and ε.

Additional equations to evaluate the listed parameters are not provided. The Maxwell model was not developed for MHs, and the values of the parameters are not supposed to change. The value of k_s is usually available or could be easily measured if a bulk sample of the solid material is available; regarding k_g, in the case of hydrogen, there are several sources in the literature, so it can be considered easily achievable. Instead, the initial value of ε can be calculated, while other models are needed for the evaluation during the hydrogenation/dehydrogenation processes. No comparisons of the Maxwell model to experimental data are reported in [42,48]. However, in [91], where a different model is presented, which will be analyzed further in this section, the authors compare experimental data with several models, including the Maxwell model. According to [91], the agreement between the Maxwell model and experimental data for LaNi₅-air and Fe-air systems is not very good, and in both cases, the ETC is underestimated.

3.1.2. Yagi and Kunii Model (1957)

In 1957, Yagi and Kunii separated into two terms the effects on the effective thermal conductivity of a packed bed of particles surrounded by a fluid phase: one independent of fluid flow and another dependent on the lateral mixing of the fluid in the packed beds [92]. The field was simplified to a packed bed with motionless fluid. In Figure 10, the model of heat transfer is reported. Equation (9) describes the model:

$$k_{eff}={\displaystyle \frac{l_p\left(1-\epsilon \right)}{\frac{l_s}{k_s}+\frac{l_v}{k_g+l_v{h}_{rs}}}}+\epsilon l_p{h}_{rv}$$

(9)

Here, l_p is the effective length between the centers of two adjacent particles, l_s is the effective length of the solid particles related to the heat conduction, and l_v is the effective thickness of the fluid film adjacent to the contact surface of two solid particles, while h_rv and h_rs are the radiative heat transfer coefficients for solid-to-solid and void-to-void, respectively, and are given by the following Equations (10) and (11):

$$h_{rs}=0.1952{\displaystyle \frac{e}{2-e}}{\left({\displaystyle \frac{T}{100}}\right)}^3$$

(10)

$$h_{rv}=0.1952{\left(1+{\displaystyle \frac{\epsilon }{2\left(1-\epsilon \right)}}{\displaystyle \frac{1-e}{e}}\right)}^{-1}{\left({\displaystyle \frac{T}{100}}\right)}^3$$

(11)

Here, e is the emissivity factor of the solid surface. The Yagi and Kunii model appears quite simple, with a limited number of parameters: k_s, k_g, ε, l_p, l_s, l_v, and e.

In addition to the considerations made for the parameters already present in the Maxwell model, the ones added here (l_p, l_s, l_v, and e) must be commented on. By their definition, assuming that particles are rigid spheres in point-contact, l_p and l_s have the exact dimension of the particle diameter, while l_v can be considered as a fitting parameter, and a value for e is not always available in the literature. However, the contribution of radiation to heat transfer is usually negligible, and emissivity is not required. The authors compared the experimental results reported in previous papers with those calculated: the solid materials used are various and include glass, coal, diphenylamine, alumina, and sand, while the gases used are H₂, He, and air, obtaining a good agreement.

3.1.3. Zehner–Schlünder Model (1970)

The model presented by P. Zehner and E. U. Schlünder in 1970 [44] constitutes a base for many other models developed in the subsequent years, as will be addressed in the following. In Figure 11, an eighth of the unit cell of the model is represented [45]. It is made of a cylinder of gaseous phase (empty area) containing one solid particle (shaded area). Adjacent particles only touch each other at one point. Following the assumptions made by Zehner and Schlünder, there are only two parallel paths for heat conduction, one in the fluid area and one in the biphasic region, while there is no heat flow through solid particles. Concerning Figure 11, it appears clear that the influence of the biphasic region depends on the shape of the solid particle, so on A_sf. To stand for this geometrical aspect, the authors introduced the shape factor parameter, B, defined through Equation (12):

$$r^2+{\displaystyle \frac{z^2}{{\left(B-\left(B-1\right)z\right)}^2}}=1$$

(12)

Following Equation (12), for B = 0, the solid particle disappears, so there is only the gas phase. For B = 1, the solid particle becomes a sphere, while for B → infinite, the solid particle occupies the entire inner cylinder so that the biphasic region becomes a solid region. The porosity and shape factor are geometrically related by Equation (13):

$$\epsilon =1-{\left({\displaystyle \frac{B\left(3-4B+B^2+2\mathit{ln}B\right)}{{\left(B-1\right)}^3}}\right)}^2$$

(13)

Equation (13) can be approximated by applying Equation (14):

$$B=C{\left({\displaystyle \frac{1-\epsilon }{\epsilon }}\right)}^m$$

(14)

where C and m are two fitting coefficients, namely the form factor and the exponential coefficient; for C = 1.25 and m = 10/9, Equation (14) corresponds to Equation (13) [45]. The central equation of the Zehner–Schlünder model is described by Equation (15):

$$k_{eff}=k_g·\left(1-\sqrt{1-\epsilon }+{\displaystyle \frac{2\sqrt{1-\epsilon }}{1-\frac{k_g}{k_s}B}}\left({\displaystyle \frac{\left(1-\frac{k_g}{k_s}\right)B}{{\left(1-\frac{k_g}{k_s}B\right)}^2}}\mathit{ln}{\displaystyle \frac{k_s}{k_{g}B}}-{\displaystyle \frac{B+1}{2}}-{\displaystyle \frac{B-1}{1-\frac{k_g}{k_s}B}}\right)\right)$$

(15)

The Zehner–Schlünder model appears quite simple, with a small number of parameters: k_s, k_g, B or C, m, and ε.

If Equation (13) is used, the only new parameter, as compared to the previous models, is the shape factor, B, which has to be treated as a fitting parameter, as it is not possible to determine a specific shape factor starting from measured data, while the porosity is not unknown. Alternatively, Equation (14) can be employed, and C, m, and ε are the unknown parameters. The porosity should be evaluated with a proper model in this latter case, while C and m are fitting parameters. The Zehner–Schlünder model was not developed for metal hydride systems but for composite materials in general, with solid-phase particles dispersed in a continuous fluid phase. Hsu et al. [45] compared the experimental data from Nozad et al. [93,94] to values calculated with the Zehner–Schlünder model. The two-phase systems were obtained by combining different solids and fluids: the solids were glass, stainless steel, urea-formaldehyde, bronze, and aluminum, while the fluids were water, glycerol, and air. The model agrees with the experimental data for k_s/k_g < 10³, while it underpredicts the effective thermal conductivity for higher values of this ratio.

3.1.4. Zehner–Bauer–Schlünder Model (1978)

In the Zehner–Bauer–Schlünder model, a similar approach to the Zehner–Schlünder model [44] is used, but some novelties are introduced. The most relevant is the area contact between adjacent particles, determined by the flattening factor, φ. This way, a third parallel heat conduction path is introduced, the one through solid particles. The flattening factor is correlated to the particle diameter, d, and the contact area diameter, d_c, through the following Equation (16) [95,96]:

$$\phi ={\displaystyle \frac{23{\left(\frac{d_c}{d}\right)}^2}{1+22{\left(\frac{d_c}{d}\right)}^{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}}$$

(16)

Equation (17) represents the primary expression of the model:

$$k_{eff}=\left(1-\sqrt{1-\epsilon }\right)k_h+\sqrt{1-\epsilon }\left(\left(1-\phi \right)k_{bp}+\phi k_s\right)$$

(17)

In Equation (17), the equivalent thermal conductivity of the biphasic region, k_bp, and that of the hollow region, k_h, appear, which are defined as described in Equations (18) and (19):

$$k_{bp}={\displaystyle \frac{2k_g}{Y-Q}}\left({\displaystyle \frac{\left(Y-\left(1+{Kn}^{*}\right)\frac{k_g}{k_s}\right)B}{{\left(Y-Q\right)}^2}}\mathit{ln}{\displaystyle \frac{Y}{Q}}-{\displaystyle \frac{B-1}{Y-Q}}\left(1+{Kn}^{*}\right)-{\displaystyle \frac{B+1}{2B}}{\displaystyle \frac{k_s}{k_g}}\left(1-\left(Y-Q\right)-\left(B-1\right){Kn}^{*}\right)\right)$$

(18)

$$k_h={\displaystyle \frac{k_g}{1+\frac{{Kn}^{*}}{\epsilon }}}+\epsilon k_g{Nu}_r$$

(19)

The parameters of Equations (18) and (19) can be calculated by the Equations (20)–(24) [95,97]:

$$Q=B\left({\displaystyle \frac{k_g}{k_s}}+{Kn}^{*}\left(1+{Nu}_r^{*}\right)\right)$$

(20)

$$Y=\left(1+{Nu}_r^{*}\right)\left(1+{Kn}^{*}\right)$$

(21)

$${Kn}^{*}=2{\displaystyle \frac{2-a}{a}}{\displaystyle \frac{2\frac{k_g}{{\mu }_g{c}_v}}{\gamma +1}}{\displaystyle \frac{l}{d}}$$

(22)

$${Nu}_r^{*}={Nu}_r{\displaystyle \frac{k_g}{k_s}}$$

(23)

$${Nu}_r={\displaystyle \frac{4\sigma T^3}{\frac{2}{e}-1}}{\displaystyle \frac{d}{k_g}}$$

(24)

Furthermore, Equations (12)–(14) are still valid. Referring to the previous equations, Kn^* and Nu_r^* are the modified Knudsen and Nusselt numbers, respectively, Nu_r is the Nusselt number, a is an accommodation coefficient, σ is the Stefan–Boltzmann constant, γ is the hydrogen heat capacity ratio, c_v is the hydrogen heat capacity at a constant volume, l is the hydrogen mean free path, and μ_g is the hydrogen dynamic viscosity.

The unit cell of the Zehner–Bauer–Schlünder model is represented in Figure 12.

The number of parameters is higher, i.e., 13 parameters, than the ones needed for the Zehner–Schlünder model, making this model more complex; they are k_s, k_g, B or C, m and ε, γ, μ_g, c_v, l, e, d, d_c, and a.

In addition to the parameters already commented for the previous models, here are some other hydrogen parameters (γ, μ_g, c_v, and l), which can be considered readily available, the diameter of the particles, d, which, if not available from the literature, might be experimentally evaluated, and two fitting parameters, d_c and a, as it is difficult to measure or evaluate the contact area diameter correctly.

As for the Zehner–Schlünder model, the Zehner–Bauer–Schlünder model was also made for a packing bed of solid particles surrounded by a fluid species. Kallweit et al. [95] used this model for the comparison to experimental data of two not-activated interstitial metal hydrides, LaNi_4.7Al_0.3 and HWT 5800 (Ti_0.98Zr_0.02V_0.43Fe_0.09Cr_0.05Mn_1.5), with inert void gases, He, Ar, N₂, and H₂; hydrogen is considered an inert gas in this case, as materials are not activated. The correspondence between the experimental and measured values is excellent.

3.1.5. Hayashi Model (1987)

In 1987, Hayashi et al. [43] proposed a model for the thermal conductivity of packed beds following the scheme of Masamune and Smith [98], which is the same as that used by Zehner, Bauer, and Schlünder [96], according to which the heat transfer in a packed bed of particles can be split into three terms:

• Conductive and radiative heat transfer in the gaseous phase;
• Conductive and radiative heat transfer in series through the gas and the solid in the biphasic region;
• Conductive heat transfer through the contact surface of solid particles.

The unit cell of the model is similar to the ones presented for the previous cases, with a cylinder enclosing two hemispheres in contact with each other (Figure 13).

The effective thermal conductivity is described by Equation (25):

$$k_{eff}={\displaystyle \frac{1}{2}}\left(3\epsilon -1\right)\left(k_{g1}+{\beta }^{\prime }{h}_{rv}d\mathit{cos}\theta \right)+{\displaystyle \frac{3{\beta }^{\prime }\left(1-\epsilon \right)\left(1-\delta \right)\mathit{cos}\theta }{2\left(\frac{1}{\frac{k_{g2}}{{\varphi }^{*}}+h_{rs}d}+\frac{\mathit{cos}\theta -{\varphi }^{*}}{k_s}\right)}}+{\displaystyle \frac{3}{2}}\left(1-\epsilon \right)\delta k_s$$

(25)

In the following, the terms appearing in Equation (25) are analyzed. The authors use the same equations introduced by Yagi and Kunii in [92] to calculate h_rv and h_rs, which are reported above in Equations (10) and (11). k_g1 is the thermal conductivity associated with mechanism one and is described by Equation (26):

$$k_{g1}={\displaystyle \frac{k_g}{1+\frac{2\beta l}{d_e}}}$$

(26)

where

$$\beta ={\displaystyle \frac{\left(2-a\right)\left(9\gamma -5\right)}{2a\left(\gamma +1\right)}}$$

(27)

Here, a is an accommodation coefficient, and d_e is the equivalent diameter for the void space, defined in Equation (28):

$$d_e=d\sqrt{{\displaystyle \frac{3\epsilon -1}{3\left(1-\epsilon \right)}}}$$

(28)

k_g2 is the thermal conductivity associated with mechanism two, conductive and radiative heat transfer in series through the gas and the solid in the biphasic region, and is equal to

$$k_{g2}={\displaystyle \frac{k_g}{1+\frac{2\beta l}{{\varphi }^{*}d}}}$$

(29)

The fractional area associated with the conductive heat transfer through the contact surface of solid particles, δ, is evaluated as shown in Equation (30):

$$\delta ={\displaystyle \frac{2k_{e0}}{3\left(1-\epsilon \right)k_s}}$$

(30)

Here, k_e0 is the limit of effective thermal conductivity as pressure tends to zero. θ is the angle corresponding to the solid–solid contact area, as described in Equation (31):

$$\theta ={\left({\mathit{sin}}^{-1}{\left({\displaystyle \frac{\delta }{n}}\right)}^{0.5}\right)}^2$$

(31)

where n is the number of contact points for a hemispherical particle surface. β′ is a factor related to the angle between the actual heat flow direction and that parallel to the axis between particles and varies from 1.0 to 0.9 for loose and close packings, respectively [43]. The authors propose the following estimation, Equation (32):

$${\beta }^{\prime }=0.9+{\displaystyle \frac{\left(\epsilon -0.260\right)\left(1-0.9\right)}{0.476-0.260}}$$

(32)

Here, 0.476 and 0.260 are the values of ε for loose and close packing, respectively. ϕ* is a dimensionless length defined as shown in Equation (33):

$${\varphi }^{*}={\displaystyle \frac{l_v}{d}}$$

(33)

Here, l_v is the effective thickness of fluid adjacent to two solid surfaces, evaluated as in Figure 13 and Equation (34).

$$l_s+l_v=d\mathit{cos}\theta$$

(34)

l_s is the solid phase length, evaluated via Equation (35):

$$l_s=d\left(\mathit{cos}\theta -{\varphi }^{*}\right)$$

(35)

By means of a linear approximation similar to that used for β′, it is possible to calculate the dimensionless length ϕ^* as exhibited in Equation (36):

$${\varphi }^{*}={\varphi }_2^{*}+{\displaystyle \frac{\left(\epsilon -0.260\right)\left({\varphi }_1^{*}-{\varphi }_2^{*}\right)}{0.476-0.260}}$$

(36)

Hayashi et al. give calculated values for the two parameters introduced in Equation (36), ϕ^*₁ and ϕ^*₂, depending on the ratio between the solid and the gas thermal conductivities. In Figure 14, the linear interpolation is reported (for numerical data, refer to [43]).

The Hayashi model appears more complex than the previous ones, with 10 parameters: k_s, k_g, k_e0, ε, γ, d, n, e, l, and a. The complexity lies in the large number of unknown parameters, and some are not quickly evaluated. Among these, in addition to the parameters common to the previous models, for which the same considerations are still valid, there are some new ones. The effective thermal conductivity at zero pressure (vacuum), k_e0, can only be obtained experimentally (unless this value is available in the literature from published experimental data on the same material); the number of contact points for a hemispherical particle surface, n, can be treated as a fitting parameter, as it is not possible to obtain this data in other ways (or it would be challenging). Parameters l_s and l_v also appear in the Yagi and Kunii model, even if they are not given a defining equation.

This model was not developed for metal hydride beds, but for solid particles in a fluid media. It has been tested for several combinations of particles and gases: glass beads (0.462 mm Ø), lead balls (1.1 mm Ø), activated alumina (0.23 mm Ø), and cylindrical PVC resin (2 mm Ø by 2 mm L, 2 mm Ø by 4 mm L), in combination with He, H₂, N₂, Ar, C₃H₆, and a He-N₂ mixture, at different temperature and pressure conditions (288.2–313.2 K and 1.33–101.3 kPa). Both for spherical and cylindrical particles, the model shows good predictions. In addition, the authors also compared it with a porous particle system, and, in this case, the correspondence is good.

In the same paper [43], the authors suggest a simplified model in which the radiative heat transfer is neglected, being very small for such systems. The same considerations made for the base model are valid, and the model’s central equation becomes as described in Equation (37):

$$k_{eff}={\displaystyle \frac{1}{2}}\left(3\epsilon -1\right)k_{g1}+{\displaystyle \frac{3{\beta }^{\prime }\left(1-\epsilon \right)\left(1-\delta \right)\mathit{cos}\theta }{2\left(\frac{{\varphi }^{*}}{k_{g2}}+\frac{\mathit{cos}\theta -{\varphi }^{*}}{k_s}\right)}}+{\displaystyle \frac{3}{2}}\left(1-\epsilon \right)\delta k_s$$

(37)

3.1.6. Sun and Deng Model (1990)

Sun and Deng presented a model for evaluating the effective thermal conductivity of metal hydride beds in 1990 [46,47]. Interestingly, this model was made for metal hydrides and not for general composite materials.

The same three heat flow terms listed for the Hayashi model are considered in this case, too. The system is made of identical quasi-spheres surrounded by a continuous gas phase, as represented in Figure 15.

The central equation of the model is Equation (38):

$$k_{eff}=\left(1-{\displaystyle \frac{4}{\pi }}\left(1-\epsilon \right)\right)\left(k_g^{*}+r·h_{rv}\right)+\left({\displaystyle \frac{\frac{4}{\pi }\left(1-\epsilon \right)\left(1-\mathit{tan}{\theta }_0\right)}{\frac{\frac{\pi }{4}-\mathit{tan}{\theta }_0}{\left(1-\mathit{tan}{\theta }_0\right)k_s^{*}}+\frac{1-\frac{\pi }{4}}{\left(1-\mathit{tan}{\theta }_0\right)k_g^{*}+\left(1-\frac{\pi }{4}\right)r·h_{rs}}}}\right)+{\displaystyle \frac{4}{\pi }}\left(1-\epsilon \right)\mathit{tan}{\theta }_0{k}_s^{*}$$

(38)

In Equation (38), the effective thermal conductivity of the gas and the solid phases, k_g* and k_s*, respectively, are introduced and described in Equations (39) and (40):

$$k_g^{*}={\displaystyle \frac{k_g}{1+2\frac{\left(2-a\right)l_0{P}_0}{a·L·P}}}$$

(39)

$$k_s^{*}=k_s\left(1+b{X}_{at}\right)$$

(40)

In Equations (39) and (40), a and b are accommodation coefficients, L is half of the distance between the centers of two particles in the no-contact direction (Figure 15), P is the current pressure, l₀ is the hydrogen mean free path at the pressure P₀, and X_at is the H-to-M atomic ratio. L can be related to ε through the following Equation (41), obtained by simple geometrical considerations (here not reported):

$$L={\displaystyle \frac{\pi d}{12\left(1-\epsilon \right)}}$$

(41)

The radiative heat transfer coefficients described in Equations (10) and (11) can be used in Equation (38). However, the authors use a slightly different numerical coefficient, 0.2269 instead of 0.1952.

Still, in Equation (38), θ₀ is the contact angle between solid particles at zero pressure. When the radiant heat transfer is omitted, it is possible to calculate this parameter from Equation (42):

$${\theta }_0={\tan }^{-1}{\displaystyle \frac{\frac{k_{e0}}{k_s}}{\frac{4}{\pi }\left(1-\epsilon \right)}}$$

(42)

This model appears to be a good compromise in terms of complexity and completeness, as the number of parameters is not as high. On the other hand, the dependence on the hydrogen concentration in the metal hydride is included. The model has 10 unknown parameters: k_s, k_g, k_e0, ε, e, d, ρ_g, l, a, and b.

Compared to the previous models, no new unknown parameters are added, but b, an accommodation coefficient, and hydrogen density, which is readily available. An interesting feature introduced by this model, reflecting that it was made for hydrides, are the equations for solid and hydrogen thermal conductivity, which are not constant, as they depend on hydrogen concentration and pressure, respectively.

In reference [47], the authors also compare the model and the experimental k_eff measured for a MlNi_4.5Mn_0.5 sample in the range 0.1–4 MPa at 40–60 °C, which shows a good agreement.

In the end, similarly to the Hayashi model [43], the authors suggest a simplified version [46,47], where the radiant heat transfer is omitted so that the main Equation (38) becomes Equation (43):

$$k_{eff}=\left(1-{\displaystyle \frac{4}{\pi }}\left(1-\epsilon \right)\right)k_g^{*}+\left({\displaystyle \frac{\frac{4}{\pi }\left(1-\epsilon \right)\left(1-\mathit{tan}{\theta }_0\right)}{\frac{\frac{\pi }{4}-\mathit{tan}{\theta }_0}{\left(1-\mathit{tan}{\theta }_0\right)k_s^{*}}+\frac{1-\frac{\pi }{4}}{\left(1-\mathit{tan}{\theta }_0\right)k_g^{*}}}}\right)+{\displaystyle \frac{4}{\pi }}\left(1-\epsilon \right)\mathit{tan}{\theta }_0{k}_s^{*}$$

(43)

3.1.7. Extended Zehner–Bauer–Schlünder Model (1994)

In 1994, Kallweit and Hane introduced some improvements to the Zehner–Bauer–Schlünder model, realizing an extended version [95]. These improvements were made to consider the interaction between the solid material and the gas, as it happens for metal hydride beds.

The unit cell and the equations are the same as the Zehner–Bauer–Schlünder model (Figure 12 and Equations (12), (14) and (17)–(24), respectively), but hydrogen concentration-dependent equations for the porosity, solid thermal conductivity, and contact radius are introduced. The central equation of the k_eff is still Equation (17).

As the metal lattice absorbs hydrogen atoms, solid particles expand; the authors use the calculations made by Peisl [99], according to which the volume variation, ΔV, recalls Equation (44):

$$∆V=2.9·{10}^{-30}{\displaystyle \frac{m^3}{H atom}}$$

(44)

According to Equation (44), it is possible to calculate the maximum expansion factor of the particles, f_ve, described in Equation (45):

$$f_{ve}=2·{10}^{-5}{\displaystyle \frac{∆V·N_A·X_{max}·{\rho }_s}{{MW}_{H_2}}}$$

(45)

Here, MW_H2 is the molecular weight of molecular hydrogen, ρ_s is the density of the metal, and X_max is the maximum hydrogen capacity in wt%. It also can be defined as the relative expansion of solid particles by Equation (46):

$${\displaystyle \frac{{∆V}_s}{V_{s,0}}}=f_{ve}{\displaystyle \frac{X}{X_{max}}}$$

(46)

Here, V_s,0 is the initial volume of the metal particle, and X is the hydrogen capacity in wt%. The porosity can be related to f_ve, supposing that the volume expansion is entirely compensated by a reduction in the porosity, being ε₀ the initial porosity, according to Equation (47):

$$\epsilon ={\epsilon }_0-\left(1-{\epsilon }_0\right)f_{ve}{\displaystyle \frac{X}{X_{max}}}$$

(47)

Similarly to the variation of the porosity, the authors propose a variation of the solid thermal conductivity with hydrogen to metal concentration, X, as referred to in Equation (48):

$$k_s=k_{s,0}-(k_{s,0}-k_s^{*}){\displaystyle \frac{X}{X_{max}}}$$

(48)

Here, k_s,0 is the solid thermal conductivity of pure metal particles, while k_s^* is the solid thermal conductivity of the hydride phase, evaluated by the authors as shown in Equation (49):

$$k_s^{*}=\left(0.7~0.8\right)·k_{s,0}$$

(49)

The extended Zehner–Bauer–Schlünder model also includes a formulation for the increase in the contact area between adjacent particles. The base idea is that it equals a value between two extreme cases: constant contact area fraction during the expansion (upper limit) and particles fully restrained between two plates (lower limit). The force between two particles corresponding to the lower limit, F₀, can be calculated with Hertz’s first formula as described in Equation (50):

$$F_0={\displaystyle \frac{4}{3}}{r}_{c,0}^3{\displaystyle \frac{E}{1-{\nu }^2}}{\displaystyle \frac{1}{r_0}}$$

(50)

Here, r_c,0 is the initial contact area radius (when no hydrogen has been absorbed), E is Young’s modulus, ν is the Poisson’s ratio, and r₀ is the initial particle radius.

F_X is the force resulting from the expansion of a particle fully restrained between two plates and can be calculated starting from Hertz’s second formula, as shown in Equation (51):

$$F_X={\displaystyle \frac{4}{3}}·{\displaystyle \frac{E·r_0^2}{1-{\nu }^2}}{\left({\displaystyle \frac{f_{ve}}{3}}·{\displaystyle \frac{X}{X_{max}}}\right)}^{{\displaystyle \frac{3}{2}}}$$

(51)

Then, introducing the force factor, f_F, it is possible to obtain an expression for the force dependent on the hydrogen concentration through F_X, according to Equation (52):

$$F=F_0+f_F{F}_X$$

(52)

The authors stated that f_F lays between 0 and 4% and can be related to ε₀ and to f_ve (with f_ve > 15.5% and 40% < ε₀ < 60%) as described in Equation (53):

$$f_F={\epsilon }_0\left(0.55-0.035f_{ve}\right)+2.285f_{ve}-35.8$$

(53)

Now that F can be calculated, it is possible to use Hertz’s first formula, Equation (50), to compute the contact radius, r_c, in dependence of X (through F), as shown in Equation (54):

$$r_c=\sqrt[3]{{\displaystyle \frac{3}{4}}F{\displaystyle \frac{1-{\nu }^2}{E}}r}$$

(54)

The number of unknown parameters is 17, which represents a higher number than the ones required for the Zehner–Bauer–Schlünder model, making this model quite complex; they are: k_s,0, k_g, ε₀, B or C, m, γ, μ_g, c_v, l, e, r_c,0, r₀, a, E, ν, ρ_s, and X_max.

The extended Zehner–Bauer–Schlünder model introduces some new parameters compared to the parent model. E, ν, ρ_s, and X_max are parameters of the alloy: X_max is always available or experimentally evaluable, E and ν are not always known, but if a bulk sample of the alloy is available, they are experimentally evaluable, and ρ_s is usually known or easily experimentally evaluable. The introduction of some mechanical properties of the solid material among the parameters is a novelty factor. Interestingly, ε and k_s vary with hydrogen concentration, and the models proposed differ from those included in the Sun and Deng model.

The authors compared the model to experimental data for LaNi_4.7Al_0.3 and HWT 5800 (Ti_0.98Zr_0.02V_0.43Fe_0.09Cr_0.05Mn_1.5) at different hydrogen pressures (in the same reference, as reported above, the base Zehner–Bauer–Schlünder model was used for comparison with the same materials but with inert gas). They noted that the typical S-shape of the k_eff curve was modified due to the interaction of the metal matrix with hydrogen. The model was able to follow this behavior quite well.

3.1.8. Modified Zehner–Schlünder: Area-Contact Model (1994)

Still starting from the Zehner–Schlünder model, Hsu et al. [45] formulated two modified versions of the model, which are analyzed in the current and the following sections.

The first is the area-contact model. As suggested by its name, in this case, a finite contact area between particles is introduced (Figure 16) because, as observed by the authors, its presence is necessary, as the effective thermal conductivity of packed beds under vacuum is not zero, so a heat conduction path through solid particles exists.

To include the finite area-contact, the deformed factor, α, is introduced so that Equation (12) is modified into Equation (55):

$$r^2+{\displaystyle \frac{z^2}{{\left(\left(1+\alpha \right)B-\left(B-1\right)z\right)}^2}}=1$$

(55)

For r = r_c and z = 1, from Equation (55), it is possible to obtain the measure of r_c, according to Equation (56):

$$r_c^2=1-{\displaystyle \frac{1}{{\left(1+\alpha B\right)}^2}}$$

(56)

It follows the main Equation (57) of the model (for details about its derivation, refer to [45]):

$$k_{eff}=\left(1-\sqrt{1-\epsilon }\right)k_g+k_s\sqrt{1-\epsilon }\left(1-{\displaystyle \frac{1}{{\left(1+\alpha B\right)}^2}}\right)+{\displaystyle \frac{2k_g\sqrt{1-\epsilon }}{1-\frac{k_g}{k_s}B+\left(1-\frac{k_g}{k_s}\right)\alpha B}}·$$

$$·\left({\displaystyle \frac{\left(1-\frac{k_g}{k_s}\right)\left(1+\alpha \right)B}{{\left(1-\frac{k_g}{k_s}B+\left(1-\frac{k_g}{k_s}\right)\alpha B\right)}^2}}·\ln {\displaystyle \frac{1+\alpha B}{\left(1+\alpha \right)B\frac{k_g}{k_s}}}-{\displaystyle \frac{B+1+2\alpha B}{2{\left(1+\alpha B\right)}^2}}-{\displaystyle \frac{B-1}{\left(1-\frac{k_g}{k_s}B+\left(1-\frac{k_g}{k_s}\right)\alpha B\right)\left(1+\alpha B\right)}}\right)$$

(57)

Similarly to Equation (13), also in this case, it is possible to obtain a correlation between the geometrical factors, α and B, and the porosity as shown in Equation (58):

$$\epsilon =1-{\displaystyle \frac{B^2}{{\left(1-B\right)}^6{\left(1+\alpha B\right)}^2}}·{\left(B^2-4B+3+2\left(1+\alpha \right)\left(1+\alpha B\right)\mathit{ln}{\displaystyle \frac{\left(1+\alpha \right)B}{1+\alpha B}}+\alpha \left(B-1\right)\left(B^2-2B-1\right)\right)}^2$$

(58)

Also in this case, Equation (14) is valid, but C and m are functions of α (the expressions of C(α) and m(α) are not reported by the authors).

The area-contact model is quite simple, even if the main equation appears complex, as the number of parameters is limited, these being the same as the base model, with the addition of α, which should be treated as a fitting parameter, as it is not possible to measure it. The list of parameters is reduced as compared with other models (Zehner–Bauer–Schlünder model [95,96], Hayashi model [43], Sun and Deng model [46,47], and Extended Zehner–Bauer–Schlünder model [95]): k_s, k_g, B or C, m, ε, and α.

The experimental k_eff data reported by Nozad et al. [93] and already compared to the numerical results for the Zehner–Schlünder model (see Section 3.1.3) were compared to the area-contact model curves too [45], taking them to a reasonable agreement also for high values of the ratio k_s/k_g (>10³). Despite the improvement of adding a finite contact area, this model was still not explicitly made for metal hydrides.

3.1.9. Modified Zehner–Schlünder: Phase-Symmetry Model (1994)

Here, the second model proposed by Hsu et al. [45], the phase-symmetry model, is analyzed. It is still based on the Zehner–Schlünder model [44] (Section 3.1.3), but the material is considered a sponge-like porous medium, with each phase continuously connected and in phase symmetry. The unit cell (Figure 17) comprises three parallel layers (fluid, biphasic, and solid), and the 1D heat conduction is the sum of the three contributions, each referring to a specific domain.

The central equation of the model is described in Equation (59):

$$k_{eff}=k_g\left(1-\sqrt{1-\epsilon }\right)+k_s\left(1-\sqrt{\epsilon }\right)+k_g\left(\sqrt{1-\epsilon }+\sqrt{\epsilon }-1\right)\left({\displaystyle \frac{B\left(1-\frac{k_g}{k_s}\right)}{{\left(1-\frac{k_g}{k_s}B\right)}^2}}\mathit{ln}{\displaystyle \frac{k_s}{k_{g}B}}-{\displaystyle \frac{B-1}{1-\frac{k_g}{k_s}B}}\right)$$

(59)

Equation (60) gives the relationship between the shape factor and the porosity:

$${\displaystyle \frac{\sqrt{\epsilon }-\epsilon }{\sqrt{1-\epsilon }+\sqrt{\epsilon }-1}}={\displaystyle \frac{B}{{\left(1-B\right)}^2}}\left(B-1-\mathit{ln}B\right)$$

(60)

Also, for this model, Equation (14) is valid. The phase-symmetry model appears quite simple, as the number of parameters is low, but it was not explicitly developed for metal hydride beds, and the assumptions made may not be valid for such systems. The list of parameters, which is the same as for the Zehner–Schlünder model, is k_s, k_g, B, or ε. No comparison to experimental data is reported for this model by the authors [45].

3.1.10. Raghavan–Martin Model (1995)

In the model proposed by Raghavan and Martin in 1995 [49], the unit cell is made of a cubic particle in a more significant cubic volume (Figure 18).

The central equation of the model is described in Equation (61):

$$k_{eff}=k_g\left(1+{\displaystyle \frac{1-\epsilon }{\frac{\frac{k_s}{k_g}}{\frac{k_s}{k_g}-1}-h_{Maxwell}·Z}}\right)$$

(61)

In Equation (61), two new parameters appear: h_Maxwell and Z. The first is a dimensionless height (the h of the unit cell, Figure 18), which takes to Maxwell’s result when it is equal to Equation (62):

$$h_{Maxwell}=1-{\displaystyle \frac{\epsilon }{3}}$$

(62)

Z, instead, is a function accounting for the geometric distortions equivalent to the distortions in isotherms and heat flow lines [48,49], which can be calculated through Equation (63):

$$Z=1+{\displaystyle \frac{\frac{k_s}{k_g}-1}{\frac{k_s}{k_g}+2}}\epsilon \left(1-\epsilon \right)\left({\displaystyle \frac{3}{2}}\epsilon A_0+A_1\left(1-\epsilon \right)\right)$$

(63)

Here, A₀ and A₁ are fitting parameters. The authors suggest, through a comparison to other models, 1/3 and 1/4 for A₀ and A₁, respectively [49].

The Raghavan–Martin model is straightforward, and only a limited number of five parameters are involved: k_s, k_g, ε, A₀, and A₁. Ghafir et al. [48] used the Raghavan–Martin model to inversely evaluate k_s starting from k_eff experimental data from Hahne et al. [100] on LaNi_4.7Al_0.3 and HWT 5800 and showed a good agreement of the model with experimental data (in terms of k_eff, as there are no experimental data for k_s).

3.1.11. Improved Area-Contact Model (2014)

Starting from the above-reported area-contact model by Hsu et al. [45], Matsushita et al. [50] proposed an improved version of the model, in which the deformed factor was calculated by assuming a simple geometrical deformation was caused by the difference between the particle expansion and the packed bed expansion (Figure 19). The unit cell is the same as reported in Figure 11 for the area-contact model.

In this regard, with α₀ and B₀ as the initial deformed and shape factors, respectively, and α_a and B_a as the deformed and shape factor after the expansion, respectively, Equation (64) can be given:

$${\alpha }_a{B}_a=\left({\left({\displaystyle \frac{1+{\varphi }_{p}R}{1+{\varphi }_s}}\right)}^{{\displaystyle \frac{1}{3}}}·\left(1+{\alpha }_0\right)-1\right)B_0$$

(64)

Here, R is the reacted fraction, while ϕ_p and ϕ_s are the particle and bed expansion ratio. This last parameter is defined in another paper from the same group [101], where they present an experimental study on the porosity of metal hydride beds. The bed expansion ratio differs among the absorption and desorption processes. The experimental formula proposed is described by Equations (65) and (66):

$$ABS: {\varphi }_s={\varphi }_{abs}·R$$

(65)

$$DES: \left\{\begin{array}{l}0<R<{\displaystyle \frac{{\varphi }_{abs}}{{\varphi }_p}}: {\varphi }_s={\varphi }_{abs}·R_{max}-\left(1-{\displaystyle \frac{{\varphi }_p}{{\varphi }_{abs}}}R\right){\varphi }_{des}\\ {\displaystyle \frac{{\varphi }_{abs}}{{\varphi }_p}}<R<1: {\varphi }_s={\varphi }_{abs}·R_{max} \end{array}\right.$$

(66)

In Equations (65) and (66), three new parameters appear: R_max, the maximum reacted fraction, ϕ_abs, the expansion ratio, and ϕ_des, the contraction ratio. Usually, R_max is meant to be 1 (the maximum value of the reacted fraction), but the reference paper defines it as “the maximum hydrogen storage in each cycle”. This definition is more appropriate for X than for R. From this interpretation, it could be calculated as the ratio between the equilibrium concentration at the externally imposed conditions, X_eq, in terms of T and P (not the T and P at every time, as the externally imposed should be the equilibrium ones, taking to the highest X value over the cycle) and the maximum reachable hydrogen concentration for the material, X_max. According to this argumentation, the maximum reacted fraction is defined in Equation (67):

$$R_{max}={\displaystyle \frac{X_{eq}}{X_{max}}}$$

(67)

The Equations (68) and (69) for ϕ_abs and ϕ_des are, instead, directly given by Matsushita et al. [101]:

$${\varphi }_{abs}={\displaystyle \frac{V_1-V_0}{V_0}}$$

(68)

$${\varphi }_{des}={\displaystyle \frac{V_1-V_2}{V_2}}$$

(69)

Here, V₀, V₁, and V₂ are the minimum, maximum, and end-of-the-cycle bed volumes. Differently from what is expected, here, the authors give Equation (70) for the reacted fraction:

$$R=\left\{\begin{array}{cc}{b}_1-\sqrt{{\displaystyle \frac{c_1-P}{a_1}}}& R<R_P\\ {\displaystyle \frac{1}{2}}+{\displaystyle \frac{P-P_{eq}}{\frac{dP}{dR}}}& R_P<R<1-R_P\\ b_2+\sqrt{{\displaystyle \frac{P-c_2}{a_1}}}& 1-R_P<R\end{array}\right.$$

(70)

where the dependence on parameters are described by Equations (71)–(75):

$$a_1={\displaystyle \frac{1}{R_P^2}}\left(P_{eq}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{dP}{dR}}\right)$$

(71)

$$b_1={\displaystyle \frac{dP}{dR}}{\displaystyle \frac{1}{2a_1}}+R_P$$

(72)

$$c_1=a_1{b}_1^2$$

(73)

$$b_2=1-R_P-{\displaystyle \frac{dP}{dR}}{\displaystyle \frac{1}{2a_1}}$$

(74)

$$c_2=P_{eq}+\left({\displaystyle \frac{1}{2}}-R_P\right){\displaystyle \frac{dP}{dR}}-a_1{\left(1-R_P-b_2\right)}^2$$

(75)

In Equations (71)–(75), P_eq is the equilibrium pressure, while R_P has been interpreted as the reacted fraction at the plateau beginning (the definition is not present in the reference paper [50] but is consistent with the equations).

The central equation of the improved area-contact model is the same as the area-contact model, Equation (57), just with α_a and B_a instead of α and B, leading to Equation (76); other parameters are also re-defined, as follows:

$$k_{eff}=\left(1-\sqrt{1-\epsilon }\right)k_g+\sqrt{1-\epsilon }\left(1-{\displaystyle \frac{1}{{\left(1+{\alpha }_a{B}_a\right)}^2}}\right)k_s+{\displaystyle \frac{2k_g\sqrt{1-\epsilon }}{1-\frac{k_g}{k_s}{B}_a+\left(1-\frac{k_g}{k_s}\right){\alpha }_a{B}_a}}·$$

$$\left({\displaystyle \frac{\left(1-\frac{k_g}{k_s}\right)\left(1+{\alpha }_a\right)B_a}{{\left(1-\frac{k_g}{k_s}{B}_a+\left(1-\frac{k_g}{k_s}\right){\alpha }_a{B}_a\right)}^2}}·\mathit{ln}{\displaystyle \frac{1+{\alpha }_a{B}_a}{\left(1+{\alpha }_a\right)B_a\frac{k_g}{k_s}}}-{\displaystyle \frac{B_a+1+2{\alpha }_a{B}_a}{2{\left(1+{\alpha }_a{B}_a\right)}^2}}-{\displaystyle \frac{B-1}{\left(1-\frac{k_g}{k_s}{B}_a+\left(1-\frac{k_g}{k_s}\right){\alpha }_a{B}_a\right)\left(1+{\alpha }_a{B}_a\right)}}\right)$$

(76)

This expression is very long, and several parameters appear. B_a and α_a have already been analyzed. The porosity is related to R, ϕ_p, and ϕ_s via Equation (77):

$$\epsilon =1-\left(1-{\epsilon }_0\right){\displaystyle \frac{1+{\varphi }_{p}R}{1+{\varphi }_s}}$$

(77)

The Eucken equation, Equation (78), is used to calculate the hydrogen thermal conductivity:

$$k_g=k_{g,Eucken}=\mu \left(c_p+{\displaystyle \frac{5R_{gas}}{4{MW}_{H_2}}}\right)$$

(78)

Here, c_p is the specific heat at a constant pressure. At a relatively low pressure, the authors suggest including the Smoluchowski effect:

$$k_g=k_{g,Eucken}{\displaystyle \frac{1}{1+2\frac{\left(2-a\right)\left(9\gamma -5\right)}{2a\left(\gamma +1\right)}·\frac{l}{d}}} at P<10 \mathrm{b}\mathrm{a}\mathrm{r}$$

(79)

Here, the already mentioned particle diameter, d, and hydrogen mean free path, l, are present, while a is an accommodation coefficient; k_g,Eucken is the hydrogen thermal conductivity calculated with Equation (78). An equation to evaluate how the particle diameter changes with the reacted fraction is also given (linear volume expansion):

$$d=d_0{\left(1+{\varphi }_{p}R\right)}^{{\displaystyle \frac{1}{3}}}$$

(80)

The improved area-contact model appears very complex, with 17 parameters, but it is also very detailed, as many equations besides the k_eff equation are given. The list of the 17 unknown parameters is k_s, ε₀, α₀, B₀, α_a or B_a, ϕ_p, V₀, V₁, V₂, P_eq, R_P, d₀, l, γ, a, μ, and c_p.

As already said for the previous models, the parameters related to hydrogen are straightforward to find in the literature. Regarding α₀, B₀, α_a, and B_a, only one equation relates these parameters to each other, so the others should be treated as fitting parameters, as it is challenging to evaluate the initial ones, α₀ and B₀, and almost impossible to evaluate the during-process ones, α_a and B_a. The volumes V₀, V₁, and V₂ could be easy to evaluate, following the procedure described in [50]. The particle expansion ratio ϕ_p is usually not available and should be experimentally evaluated, even if, for some “common” alloy, this value can be found in the literature. P_eq and R_P are known, as they can be deduced from the PCI curves, which are usually available in the literature; for new alloys, the PCIs need to be experimentally realized. One more comment is dedicated to the reacted fraction, R. This parameter is modeled through Equation (70), while in general, it is treated as a variable.

The comparison of the experimental data for the LaNi₅ to the improved area-contact model and a version of the model with no change in the contact area is reported in [50], showing how the introduction of a variation of the contact area influences the results, especially at high values of P, leading to a better predictivity.

3.1.12. Abdin–Webb–Gray Model (2018)

In 2018, Abdin et al. [102] developed a model for a metal-hydride tank for hydrogen storage made up of several sub-models. Among them, an effective thermal conductivity model was included. They started from the result of Gusarov and Kovalev [103], who obtained Equation (81):

$$k_{eff}={\displaystyle \frac{N\left(1-\epsilon \right)}{\pi d·R_c}}$$

(81)

where N is the particle coordination number, and R_c is the thermal contact resistance. The parameters of Equation (81), N, ε, d, and R_c, have been evaluated by Abdin et al. [102] as follows. A linear volume expansion is assumed, so Equation (80) also expresses the relation between d and R in this case. N, instead, is evaluated with an empirical expression found by van de Lagemaat et al. [104], as shown in Equation (82):

$$N={\displaystyle \frac{3.08}{\epsilon }}-1.13$$

(82)

The geometrical model used to evaluate R_c is represented in Figure 20.

The contact between two particles can be divided into two different scale regions: the macro-gap and the micro-gap. The heat flows through the macro-gap, encountering two parallel resistances: the gas (R_G) and the contact between solid particles. This last one is divided into two in-series contributions. One is given by assuming a perfect contact between particles (R_L) and the other, looking at the micro-gap scale, by including the micro-contacts (R_s) in parallel with the resistance of the interstitial gas contained in the micro-gaps (R_g). In this way, R_c can be described by Equation (83):

$${\displaystyle \frac{1}{R_c}}={\displaystyle \frac{1}{\frac{R_s{R}_g}{R_s+R_g}+R_L}}+{\displaystyle \frac{1}{R_G}}$$

(83)

The resistance R_L is described in the Equation (84):

$$R_L={\displaystyle \frac{1}{2k_s{r}_c}}$$

(84)

The model for the thermal resistance of micro-contact is the same one developed by Bahrami et al. [105] and shown in Equation (85):

$$R_s={\displaystyle \frac{0.565H_v{d}_v}{k_s{F}_n}}$$

(85)

In Equation (85), F_n is the normal contact force, H_v is the Vickers microhardness, and d_v is the mean indentation diagonal depth (the depth from the sample’s surface to the indenter tip). Introducing C₁ and C₂, the Vickers microhardness coefficients, H_v can be expressed as in Equation (86) [106]:

$$H_v=C_1{d_v}^{C_2}$$

(86)

The normal contact force, instead, can be calculated as expressed in Equation (87):

$$F_n={\displaystyle \frac{4}{3}}{E}^{\prime }\sqrt{r·{d_v}^3}$$

(87)

Here, E′ is the effective Young’s modulus, and r is the particle radius. The macro-contact radius, r_c, can be calculated as in Equation (88) [107]:

$$r_c=\left\{\begin{array}{c} {\displaystyle \frac{1.605}{\sqrt{\epsilon }}}{\left(0.75{\displaystyle \frac{F_{n}r}{E^{\prime }}}\right)}^{{\displaystyle \frac{1}{3}}} 0.01\le \epsilon \le 0.47\\ \left(3.51-2.51\epsilon \right){\left(0.75{\displaystyle \frac{F_{n}r}{E^{\prime }}}\right)}^{{\displaystyle \frac{1}{3}}} 0.47<\epsilon \le 1 \end{array}\right.$$

(88)

Also, the micro-gap resistance, R_g, has been obtained by Bahrami et al. [105] according to Equation (89):

$$R_g={\displaystyle \frac{\sqrt{2}{\sigma }_R{a}_2}{\pi k_g{r}_c^2\mathit{ln}\left(1+\frac{a_2}{a_3+\frac{M}{\sqrt{2}{\sigma }_R}}\right)}}$$

(89)

Here, the surface roughness, σ_R, two operators, a₂ and a₃, the contact microhardness, H_c, and the gas parameter, M, appear. They can be calculated as described in Equations (90)–(93) [102,106]:

$$a_2={erfc}^{-1}\left(0.03{\displaystyle \frac{P_{load,max}}{H_c}}\right)-a_3$$

(90)

$$a_3={erfc}^{-1}\left(2{\displaystyle \frac{P_{load,max}}{H_c}}\right)$$

(91)

$$H_c=C_1{\left(1.62d_v\right)}^{C_2}$$

(92)

$$M=\left({\displaystyle \frac{2-{\alpha }_{T_1}}{{\alpha }_{T_1}}}+{\displaystyle \frac{2-{\alpha }_{T_2}}{{\alpha }_{T_2}}}\right)\left({\displaystyle \frac{2\gamma }{1+\gamma }}\right){\displaystyle \frac{1}{Pr}}l$$

(93)

In Equations (90) and (91), the maximum contact load pressure, P_load,max, has been introduced, which can be calculated through Equation (94):

$$P_{load,max}={\displaystyle \frac{2}{\pi }}{E}^{\prime }\sqrt{{\displaystyle \frac{d_v}{r}}}$$

(94)

In Equation (93), α_T1 and α_T2 are thermal accommodation coefficients, depending on the gas–solid combination, and Pr is the Prandtl number. The hydrogen thermal conductivity is taken from [108] and calculated with Equation (95):

$$k_g={\displaystyle \frac{k_{g,ref}}{1+2b_{g}Kn}}$$

(95)

Here, k_g,ref is the hydrogen thermal conductivity at the gas reference pressure (1 bar), b_g is a constant depending on the gas and is equal to 9.87 for hydrogen, and Kn is the Knudsen number, which can be calculated as shown in Equation (96):

$$Kn={\displaystyle \frac{l}{d}}$$

(96)

In [108], it is also given an equation to calculate the mean free path as shown in Equation (97):

$$l={\displaystyle \frac{6\left(\gamma -1\right)}{9\gamma -5}}·{\displaystyle \frac{k_g}{P}}\sqrt{{\displaystyle \frac{{MW}_{H_2}·T}{2k_B}}}$$

(97)

In the end, R_G can be calculated as in Equation (98) [103]:

$$R_G={\displaystyle \frac{1}{2\pi k_{g}d\left(\frac{1}{2}\mathit{ln}\left(1+L\right)+\mathit{ln}\left(1+\sqrt{L}\right)+\frac{1}{1+\sqrt{L}}-1\right)}}$$

(98)

with L (thickness) described by Equation (99):

$$L={\displaystyle \frac{\gamma +1}{9\gamma -5}}·{\displaystyle \frac{3}{4Kn\sqrt{\pi }}}$$

(99)

Substituting Equations (84), (85), (89), and (98) in Equation (83) and then Equation (83) in (81), it is possible to calculate k_eff.

Concerning the porosity, Abdin et al. [102] propose to use the same equations given by Matsushita et al. [50], i.e., Equations (65), (66) and (77); even if they do not mention Equations (67)–(75), it is possible to assume that such equations are still needed, or R might be treated as a variable.

This model is complex and characterized by many parameters and equations, many of which are taken from different references. It also appears very interesting because of the approach, which is entirely different from that used in the above-presented models, with the distinctions of the contributions from the macro-gaps and the micro-gaps, which are strongly related to the material’s mechanical properties. The unknown parameters are 20 as listed: k_s, k_g,ref, ε₀, d₀, d_v, C₁, C₂, E′, σ_R, Pr, α_T1,α_T2, γ, b_g, P_eq, R_P, ϕ_p, V₀, V₁, and V₂. Apart from the parameters that have already been commented on, here are some new parameters. C₁, C₂, E′, σ_R, and d_v are directly related to the mechanical properties of the solid material. Often, these properties are unknown for the alloys used in MH systems, commonly sold as powders, while a bulk sample would be needed for the experimental characterization. Pr is easy to evaluate (in [102], a value of 0.6 is suggested), as well as k_g,ref.

Two comparisons with experimental data are reported in [102]. The first is with LaNi₅ data from Pons et al. [109], demonstrating a good agreement, also in comparison to a model published by the same authors [110]; the second is with a LaNi_4.7Al_0.3 alloy, and in this case, the model from Abdin et al. showed a better approximation to experimental data than other three models (the above presented Extended Zehner–Bauer–Schlünder model and Sun and Deng model, and the model from Asakuma et al. [111]).

3.1.13. Heat Transfer Concentrating Model (2023)

Bai et al. [91] proposed a heat transfer concentrating (HTC) model based on the observation that small gas volumes are responsible for a high quantity of transferred heat. According to the authors, only a minimal amount of heat transfer happens in gas channels, as all the channels in the direction of the heat flow are interrupted by solid surfaces, as shown in Figure 21. In this way, the pure gas channel turns into the pathway of particle–gas film–particle, and the contribution of pure gas channels can be neglected. Furthermore, as the measured effective thermal conductivity in a vacuum is usually extremely low, the contribution of the contact between particles can also be neglected.

The central expression of the heat concentrating model is described in Equation (100):

$$k_{eff}={\displaystyle \frac{k_s{k}_g\left(\frac{\epsilon }{G}+1-\epsilon \right)}{\frac{\epsilon }{G}·k_s+k_g\left(1-\epsilon \right)}}$$

(100)

Compared to the equation reported in the reference [91], Equation (100) appears quite different to the reader because, apart from the different symbols used, some substitutions have been applied to obtain a compacted form.

In Equation (100), the effective gas coefficient, G, appears. It is defined as the reciprocal of the effective gas film volume ratio, η₃, which is the ratio between the gas volume of the effective gas film region, V_gas,eff, and the total gas volume in the characteristic cubic unit, V_gas,all (refer to Figure 21), as shown in Equation (101):

$$G={\displaystyle \frac{1}{{\eta }_3}}$$

(101)

$${\eta }_3={\displaystyle \frac{V_{gas,eff}}{V_{gas,all}}}={\displaystyle \frac{V_{gas,cylinder}\left|\begin{array}{c}\tau =\omega :\left({\eta }_2=0.9\right)\\ \tau =0 \end{array}\right.}{V_{gas,all}}}$$

(102)

In Equation (102), some new parameters are introduced: τ, the angle from the center of the spherical particle to the annular cylinder, η₂, the ratio between the heat quantity transferred through the cylinder region (τ from 0 to ω) and the total heat quantity transferred through the characteristic cubic unit, and ω, the angle of the cylinder region. Imposing a value of ω so that η₂ = 0.9 corresponds to finding the volume of the effective gas film region, V_gas,eff, which is the gas volume in the cylinder, V_{gas, cylinder}, contributing to 90% of the heat quantity transfer. The expression for η₂ is described by Equation (103):

$${\eta }_2={\displaystyle \frac{\mathit{ln}\frac{k_s}{k_g}+\mathit{ln}\left(1-\mathit{cos}\omega \right)-\left(1-\mathit{cos}\omega \right)}{\mathit{ln}\frac{k_s}{k_g}+\frac{4-3\pi }{2\pi }}}$$

(103)

Taking into consideration the geometry of the system (Figure 21), the successive Equations (104)–(106) describe such characteristics:

$$V_{gas,cylinder}=\pi {\left(r·\mathit{sin}\omega \right)}^2·d-\left(2\pi h_{sc}^2\left(r-{\displaystyle \frac{h_{sc}}{3}}\right)+2\pi {\left(r·\mathit{sin}\omega \right)}^2\left(r-h_{sc}\right)\right)$$

(104)

$$V_{gas,all}=d^3-{\displaystyle \frac{4}{3}}\pi r^3$$

(105)

$$h_{sc}=r\left(1-\mathit{cos}\omega \right)$$

(106)

The HTC model is simple from a mathematical point of view, and the number of parameters is limited, but it has been designed without including the particle expansion typical of hydrides. The four unknown parameters are k_s, k_g, ε, and d. No new unknown parameter is introduced in the HTC model. Here, ω is not listed as an unknown parameter as it can be calculated by Equation (103) imposing η₂ = 0.9. The comparison of the HTC model to experimental data reported in [93,100,112] and to results from several other models (Maxwell, Yagi and Kunii, Kunii and Smith, Zehner–Schlünder, Abrahamsen, and Geldart, and improved area-contact) gives interesting results. The HTC model best agrees with a comprehensive class of materials regarding solid particles and gases. The tested solid materials are Al₂O₃, steel, Si, SiC, Cr, Al, Fe, LaNi₅, LaNi_4.7Al_0.3, and HWT 5800 with air, Ar, He, N₂, and H₂. Even if the results are excellent, it needs to be underlined that the comparison was made for not MH systems or, in the case of LaNi_4.7Al_0.3 and HWT 5800 with H₂, for not-activated hydrides (in [100] experimental data for both activated and not-activated systems are reported, but it appears that data used in [91] are those of not-activated systems).

3.2. On the Thermal Conductivity of the Solid

As described in the previous sections, the hydride-forming metal/alloy falls into very small pieces during cycling due to the swelling upon hydrogen absorption and the subsequent shrinkage during desorption. The apparent thermal conductivity of both the hydrogen gas and solid phase particles together constitutes the ETC. The thermal conductivity of the gas phase is related to its actual pressure. The solid particles may have a two-phase microstructure, consisting of a hydride phase and an unhydrided phase or partial solid solutions of hydrogen in the alloy phase. Three models contemplate the effect of the concentration of hydrogen in the solid phase and the behavior of what is called the solid thermal conductivity of the hydride. Sun and Deng’s model (Section 3.1.6) [46,47] and the extended Zehner–Bauer–Schlünder’s model (Section 3.1.7) [95] introduced the k_s* in Equations (40) and (48), respectively, to account for the change in the thermal conductivity of the solid phase during the phase transition from metal/alloy to hydride. In the literature, there are two works [48,63] in which the thermal conductivity of the solid (alloy to hydride) is calculated through an inverse methodology from the ETC and applying models. Ghafir et al. [48] calculated the solid thermal conductivity from the ETC by applying the Yagi and Kunii model (Section 3.1.10) [92] to LaNi_4.7Al_0.3H_x and HWT 5800—Ti_0.98Zr_0.02V_0.43Fe_0.09Cr_0.05Mn_1.5H_x [49]. Kumar et al. [63] used the same strategy but with the Yagi and Kunii model (Section 3.1.2) [92] to MmNi_4.5Al_0.5. An increase in the solid thermal conductivity was observed during the transition from the alloy to the hydride compound. Both models (Section 3.1.2 and Section 3.1.10) were not developed for hydrides, and this is an aspect to be considered when analyzing the results. There is practically no review work in the field of hydrides that discusses this important aspect. Therefore, in this section, the main question to address is whether the thermal conductivity of the hydride-forming metal/alloy is lower or higher than that of the solid metal hydride (at the micrometric/nanometric level). The authors propose an analysis of an isolated small particle without considering the gas phase by an analysis of the fundamental physical principles governing the conduction of heat in solids.

Initially, the heat conduction mechanism for the perfect crystalline structures of solids is presented. In the second step, imperfections in the crystalline structure of solids are introduced, and their effects on the thermal conductivity of solids are discussed. Lastly, the thermal conductivity of metal hydrides is considered in this context. Referring to the Fourier law and as already described by Equation (2) in the axial direction (Section 2.1.1), the heat flux transported by conduction in all directions can be defined as in Equation (107):

$$Q=k\nabla T$$

(107)

Here, Q is the heat flux perpendicular to an area in the solid with a temperature gradient of ∇T. In general, thermal conductivity k is a tensor; however, for the sake of simplicity, it will be treated in the following as a scalar (assuming isotropic cubic crystal symmetry).

Various energy carriers can transfer thermal energy in solids. These can be phonons (lattice vibrations), electric charges (electrons and holes, among the most relevant), excitons (electron–hole pairs), photons (electromagnetic waves), plasmons (collective electron oscillations), spin waves, and other excitation processes. Therefore, the thermal conductivity of a solid is the sum of all these heat conduction channels, as described by Equation (108):

$$k={\sum }_i{k}_i$$

(108)

where i indicates a specific heat conduction channel. It must be noted that Equation (108) does not consider the cross-interaction between different heat conduction channels (no cross-terms). In reality, the electrons are scattered by phonons or photons, and phonons scatter photons. These cross-interactions are inharmonic and require a third or higher order of perturbation terms. The contribution of these processes to the overall thermal conductivity can be accounted for by perturbation theory [113,114,115,116,117,118] and the Born–Oppenheimer approximation [119,120,121,122,123,124], assuming an adiabatic perturbation [125,126,127]. A more comprehensive understanding of these processes can only be adequately addressed using quantum statistical methods. Readers who are more interested in a rigorous theoretical treatment of this subject are referred to the literature in references [128,129,130,131,132,133]. Another valid simplification to Equation (108) can be made based on the observation that, in most cases, heat conduction in solids is predominantly facilitated by phonons and electrons. Equation (109) describes this simplification:

$$k=k_p+k_e$$

(109)

In the first step, we will address the mechanism of heat conductivity through phonons (k_p), followed by the mechanism of heat transport via electrons (k_e). The final section will address the comparative implications of these mechanisms for different categories of materials.

In 1911, Arnold Euken [134] experimentally observed very high thermal conductivity for diamond, which was considerably more significant than silver or copper. Furthermore, he reported an inverse proportional relationship between the temperature and its thermal conductivity (k~T⁻¹) at elevated temperatures. Peter Debye [135], Max Born [136,137], and others [138,139] could theoretically derive the T⁻¹ thermal conductivity dependency is based on the assumption that heat in solids is carried by a quantum of lattice vibration field (phonons). In particular, Debye [140] could derive, based on the kinetic theory of monoatomic gases, a relatively simple expression for thermal conductivity that could reproduce many experimental results, described in Equation (110):

$$k_p={\displaystyle \frac{1}{3}}{c}_v^p{v}_p{\mathsf{\Lambda }}_p={\displaystyle \frac{1}{3}}{c}_v^p{v}_p^2{\mathsf{\tau }}_p$$

(110)

Here, c_v^p stands for the total heat capacity of the phonon gas (c_v^p refers to the volume-specific heat capacity of the phonon gas and c_v^p = n_p·c_p, where n_p is the concentration of phonons, and c_p is their heat capacity), v_p is the average group velocity of the phonons and Λ_p represents the particle mean free path between two collision events (Ʌ_p = v_p·τ_p, where τ_p is the relaxation time of a phonon). The heat capacity of phonons can be approximated at lower temperatures with the Debye law, as described in Equations (111) and (112) [141]:

$$c_v^p\approx {\displaystyle \frac{1}{3}}{\alpha }_{Debye}{T}^3$$

(111)

with

$${\alpha }_{Debye}={\displaystyle \frac{234N_{atoms}{k}_B}{{\theta }_D}}$$

(112)

N_atoms is the number of atoms in the sample, k_B is the Boltzmann constant, and θ_D is the Debye temperature. With an increasing temperature, the phononic heat capacity reaches the limit of high temperatures, the Dulong–Petit value of 3Natoms·k_B. It is essential to recognize that Debye’s theory was predicated on elastic scattering among phonons, leaning on the classical notion of the kinetic theory of monoatomic gases (phonon “gas”), where energy and momentum are conserved. Furthermore, he assumed a linear dispersion relationship (for all phonon polarizations) in the solid (ω_f~K, ω_f and K being the angular frequency and the wave-vector of the phonons, respectively). This assumption, however, implies an elastic continuum where the group velocity is independent of the frequency. It is worth mentioning that the exact derivation can also be made for other excitation channels as well, which leads to the same expression as in Equation (110) [142,143,144]. Peierls [145] showed that Debye’s assumptions were rough approximations, where no physical mechanism was involved in the theory that could lead to thermal resistance in a solid body. He demonstrated that viewing a solid as an elastic continuum would be insufficient for an accurate theoretical representation of the thermal conductivity of a crystalline solid body. A crystalline solid should, instead, be described as a discontinuous lattice, respecting the boundary conditions set by its periodicity. By expanding the perturbation terms for phonons to third order, including all the possible three phonon–phonon scatterings, Peierls [145] was able to demonstrate a mechanism (specified in Equations (111) and (112)) that leads to a statistical equilibrium among different oscillatory states in a solid according to Equation (113):

$${\mathit{K}}_{\mathbf{1}}+{\mathit{K}}_{\mathbf{2}}={\mathit{K}}_{\mathbf{3}}+\mathit{G}$$

(113)

Here, Κ_1,2,3 is the wave vector of the ith phonon (I = 1,2,3), and G is the reciprocal lattice vector. One possible process is obtained when the reciprocal lattice vector vanishes (G = 0). This is the so-called “normal” process (N-process), where the momentum and energy of the phonons are preserved. The following process is possible when the reciprocal lattice vector does not vanish (G ≠ 0), which is the central statement of Peierls’ theory [145], stating that if the two phonons are scattered, their total momentum exceeds the first Brillouin zone. Therefore, the lattice vector G has the “obligation” to fold the resulting phonon vector back into the first Brillouin zone. This is the so-called “Umklapp” process or short U-process (Umklappen is a German word that means to fold over). Also, in the U-process, the energy is preserved. However, it does not preserve the momentum, providing the driving force for an equilibrium among different oscillatory states in the solid. Both processes are shown schematically in Figure 22. Therefore, the relaxation time τ_p of the phonons can be constructed from two distinct additive relaxation times, namely the relaxation time of the N-process and the relaxation time of the U-process.

The relaxation time for the U-process, τ_U, was derived by Klemens et al. [143] and others, which is proportional to an exponential decay with temperature, as shown in Equation (114) [142,146,147]:

$${\tau }_U^{-1}~T·e^{{\displaystyle \frac{-{\theta }_D}{\alpha ·T}}}$$

(114)

This process becomes more important at higher temperatures, since the wavelengths of the phonons are small enough to exceed the first Brillouin zone in the scattering process. If we move away from a perfect crystal and consider real crystalline solids, then additional scattering processes need to be taken into account. In the first step, a real solid is spatially limited by its boundaries, and phonons can scatter at these boundaries (grain boundaries, interfaces, and surfaces). The corresponding relaxation time for the phonon-boundary scattering, τ_P-B, is given by Equation (115) [148]:

$${\tau }_{P-B}^{-1}~{\displaystyle \frac{v_s}{D}}$$

(115)

In this equation, vs is the speed of sound inside the solid with which the phonons are propagating, and D is the size of the boundary of the solid (the size of the grain in a polycrystalline solid or the thickness of a thin film, for instance). For a bulk crystal, the phonon-boundary effect becomes only dominant at very low temperatures (below the Debye temperature, θ_D), as the wavelength of the phonons can be compared to the size of the crystal. In contrast, in nanocrystals, the phonon-boundary effect only becomes dominant at very high temperatures since the wavelengths of the phonons become comparable or smaller than the size of nanocrystals. In both cases, the mean free paths (MFP) of phonons are reduced, hence posing a thermal resistance that reduces the overall thermal conductivity of the solid. Furthermore, real crystalline solids can contain impurities, defects, and isotope variations, to which the phonon field can couple. The most significant hampering effect on the phonons MFP, thereby reducing the overall thermal conductivity, is delivered by the randomly distributed isotopes in a crystal solid. The relaxation time for the phonons MFP, τ_P-is, is given by Equation (116) [149]:

$${\tau }_{P-is}^{-1}~{\left({\displaystyle \frac{∆M}{M_{av}}}\right)}^2{\omega }_f^4$$

(116)

Here, ΔM is the difference in mass between the two isotopes, M_av is the mass of the average atomic mass, and ω_f is the angular moment of the phonon. This scattering mechanism becomes more prominent at temperatures below the Debye temperature since their wavelengths become comparable to interatomic distances. There are also other scattering mechanisms originating from the lattice imperfections of crystals, such as dislocations, vacancies, interstitials, stacking faults, phonon-resonance scattering, etc. [150,151,152]. Since the relaxation times of all of them, τ_P-d, are proportional to their concentration, here we summarize them in Equation (117):

$${\tau }_{P-d}^{-1}~{\sigma }_d{n}_d$$

(117)

In Equation (117), n_d is the concentration of defects, and σ_d represents the corresponding scattering cross-section for phonons. In conclusion, the relaxation times of all phononic scattering processes can be summarized according to Matthiessen’s rule, as shown in Equation (118) [153]:

$${\tau }_p^{-1}\approx {\tau }_U^{-1}+{\tau }_{P-B}^{-1}+{\tau }_{P-is}^{-1}+{\tau }_{P-d}^{-1}$$

(118)

In Figure 23, an example of the aforementioned phononic scattering processes belonging to CoSb₃ is shown [154].

Considering the electronic contribution (k_e) to the thermal conductivity in solids, primarily metals, alloys, and semiconductors can make a significant contribution. In these materials, the electronic thermal conductivity, alongside phononic thermal conductivity, plays a role in the total thermal conductivity of the material, which, for simplicity, we consider only pure metals. The conduction electrons in metals are delocalized, which can be treated as an electron “gas” within the metal. Therefore, Equation (110) can also be applied here and expressed as in Equation (119):

$$k_e={\displaystyle \frac{1}{3}}{c}_v^e{v}_F{\Lambda }_e={\displaystyle \frac{1}{3}}{c}_v^e{v}_F^2{\tau }_e$$

(119)

Here, c_v^e refers to the electronic heat capacity per unit volume, Ʌ_e represents the particle mean free path between two collision events of electrons, τ_e describes the relaxation time of an electron, and v_F represents the Fermi velocity of the conducting electrons at the Fermi surface. The heat capacity c_v^e of the electron gas is given by Equations (120) and (121) [155]:

$$c_v^e={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\pi }^2{n}_e{k}_B^{2}T}{{\in }_F}}$$

(120)

and

$${\in }_F={\displaystyle \frac{1}{2}}{m}_e{v}_F^2$$

(121)

with n_e as the electron concentration, m_e as the electron mass, and ∈_F as the Fermi energy. Inserting Equations (120) and (121) in Equation (119), the electrical thermal conductivity of a solid is obtained, as expressed by Equations (122) and (123):

$$k_e={\displaystyle \frac{{\pi }^2{n}_e{k}_B^2{\tau }_e}{3m_e}}T$$

(122)

$$k_e={\gamma }_{Sommerfeld}T$$

(123)

with the definition of the so-called Sommerfeld constant via Equation (124):

$${\gamma }_{Sommerfeld}={\displaystyle \frac{{\pi }^2{k}_B^2{n}_e{\tau }_e}{3m_e}}$$

(124)

It is worth noting the difference in the temperature dependency between the phonon heat capacity (c_v^p~T³) in Equations (111) and (112) and the heat capacity of the electron gas (c_v^e~T). The ratio between them reflects their temperature behavior: c_v^e/c_v^p~T⁻². At very low temperatures, the electronic heat capacity is predominant, whereas at higher temperatures, the phonon heat capacity becomes more prominent. Electrical conductivity, θ_electrical, is defined as in Equation (125):

$${\theta }_{electrical}={\displaystyle \frac{n_e·{\underset{\_}{e}}^2·{\tau }_e}{m_e}}$$

(125)

Here, e is the elementary charge, equal to 1.602176634 × 10⁻¹⁹ C [156]. Combining Equations (122) and (125) results in Equation (126):

$$k_e={\displaystyle \frac{{\pi }^2}{3}}{\left({\displaystyle \frac{k_B}{\underset{\_}{e}}}\right)}^2{\theta }_{electrical}T$$

(126)

Equation (126) shows a proportionality relation between thermal and electronic conductivity, which is an essential relationship since the thermal conductivity can be indirectly determined by measuring its electronic conductivity. Such an electronic conductivity measurement has to be carried out, of course, at very low temperatures where the electronic heat capacity is predominant, and all phonon modes are mostly “frozen”, as discussed above. Equation (126) is the first approximation of the electronic thermal conductivity of a metal. Landau’s Fermi theory of liquids gives a more realistic treatment of metals. In this theory, the overall electronic relaxation time, τ_e, is split into many processes, such as electron–electron scattering processes, and is given by Equation (127):

$${\tau }_{e-e}^{-1}~T^{-2}$$

(127)

It is worth mentioning that this relaxation time is only valid for bulk metals. In the case of thin films or nanotubes, the proportionality with respect to temperature changes significantly. As a further significant interaction, the electron-impurity scattering can be considered, which has a relaxation time described by Equation (128):

$${\tau }_{e-i}^{-1}~{\displaystyle \frac{N_i·{\sigma }_{e-i}}{V}}$$

(128)

Here, V is the volume of the metal, σ_e-i is the cross-section for the scattering of the impurity, and N_i is the concentration of the impurity in the metal. Lastly, the electrons can be scattered by phonons, and their relaxation time at low temperatures is given by Equation (129) [142]:

$${\tau }_{e-p}^{-1}~T^5$$

(129)

This is only valid at temperatures below the Debye temperature (T < θ_D). For higher temperatures, the expression for τ_e-p⁻¹ becomes proportional to the temperature (T) due to the increased phonon density.

Other electronic scattering processes are not further relevant to this work. In conclusion, the overall relaxation time for the electronic thermal conductivity in metals can be summarized as follows in Equation (130):

$${\tau }_e^{-1}\approx {\tau }_{e-e}^{-1}+{\tau }_{e-i}^{-1}+{\tau }_{e-p}^{-1}$$

(130)

Combining Equation (130) with Equation (118), it is possible to obtain the total relaxation time for a crystalline metal, as expressed in Equation (131):

$${\tau }_{tot}^{-1}\approx {\tau }_p^{-1}+{\tau }_e^{-1}\approx {\tau }_U^{-1}+{\tau }_{P-B}^{-1}+{\tau }_{P-is}^{-1}+{\tau }_{P-d}^{-1}+{\tau }_{e-e}^{-1}+{\tau }_{e-i}^{-1}+{\tau }_{e-p}^{-1}$$

(131)

Considering the phononic and electronic thermal conductivities, phonons are quasi-quantum particles of the solid lattice; hence, they are present in all solids, whereas electronic thermal conductivity is primarily present in metals, alloys, and semiconductors. However, to predict the thermal conductivity of a solid, a quantitative expression of thermal conductivity is desirable. Leibfried and Schlömann [157] derived such a quantitative description explicitly, focusing on the U-processes, which was further improved by Slack [158] by excluding the contribution of optical phonon branches to the group velocity of phonons, as shown in Equation (132):

$$k_p=A_{const}{\displaystyle \frac{M_{av} {\delta }_v {\theta }_D^3}{{\gamma }_{Gr\ddot{u}}^2 N_{atoms-cell}^{3/2} T}}$$

(132)

Here, A_const is a proportionality constant, M_av is the average atomic mass, δ_v is the average atomic volume, N_atoms-cell is the number of atoms in the unit cell, and γ_Grü is the Grüneisen parameter. The Debye temperature is also proportional to the group velocity (v_p) according to Equation (133):

$$v_p=\sqrt{B_{bulk}{\delta }_v^3{M}_{av}^{-1}}$$

(133)

where B_bulk is the bulk modulus. Hence, the thermal conductivity is inversely proportional to M_av^1/2 and proportional to B_bulk^1/2. This formula is beneficial for predicting the phononic thermal conductivity of materials. It indicates that materials possessing a high phononic thermal conductivity are made of light atoms, with a crystal structure containing one atom in the unit cell, strong interatomic (stiff) bonds, and a small Grüneisen parameter. Diamond, Ge, and Si, for instance, fulfill all these requisites, and indeed, they show an extremely high thermal conductivity [159,160,161,162]. This phenomenological formula can be applied as a directive measure for the phononic thermal conductivity value of a material. However, it does not consider the details of the phonon band structure of the material [163].

When comparing the thermal conductivity of metals with their corresponding alloys, it can be stated that, in general, the thermal conductivity of metals is higher than that of their corresponding alloys. This is due to their simple crystal structure, as most of them are solidified in FCC or BCC lattice structures, except for a few (Mg, Ti), which have an HCP lattice structure, which allows a uniform propagation of phonons across the metal lattice. Moreover, the contribution of the electronic heat conductivity of metals is higher, as well, in comparison to their corresponding alloys. This is due to their delocalized conducting band electrons, which facilitate heat transfer in metals. Alloys, in contrast, are mixtures of two or more elements exhibiting more complex crystal structures with more than one atom in the unit cell. Furthermore, the introduction of different atoms in the metal lattice can create significantly higher defects (vacancies, dislocations, stacking faults, precipitates, etc.), varying the electron density, and phase variation within the crystal lattice. All these imperfections provide significantly higher scattering sites for phonons and electrons in alloys relative to metals [164]. Similarly, introducing hydrogen into a metal lattice can create a more complex crystal structure with more than one atom in the unit cell. It also introduces more imperfections through embrittlement, blister formation, flacking, and dislocations, among other imperfections [165], which again provide additional scattering sites for phonons and electrons. Therefore, in general, it can be stated that metals are better thermal conductors in comparison to their corresponding alloys and their counterparts in the hydride state. Furthermore, in intermediate stages, two phases (hydride/dehydrided) with very different heat conduction mechanisms can be present within one particle of the hydride bed. While the fraction of each phase can be correlated with the phase transformation kinetics and resembled by respective models, the distribution of phases within the microstructure of a particle is challenging to account for, but may have a significant influence on the apparent thermal conductivity of the solid phase. E.g., upon hydriding, nuclei of the hydride phase may form as islands in an otherwise continuous metallic alloy, or a continuous hydride phase develops all along the crystallite boundaries, where hydrogen diffusion is fastest, then growing into the crystallites from there. In the first case, heat conductivity will be governed by the metallic phase, in the latter, the hydride phase will control it. A similar effect applies to the desorption case. Depending on the kinetic rate-limiting step, the appropriate phase transformation model may be applied to approximate the heat conductivity in intermediate stages as a function of the transformed fraction a.

4. Effect of the Pressure, Temperature, and Composition on the ETC

This section analyses how the effective thermal conductivity of a metal hydride bed varies with the main hydride’s operation variables, which are the pressure (P), hydrogen concentration (X), and temperature (T), for the models described in Section 3.1. Moreover, the experimental results obtained in some papers [63,100,166] are analyzed and compared with the presented models [42,43,44,45,46,47,49,50,91,92,95,96,97,102]. Furthermore, the dependence of the ETC on the hydride composition is not contemplated in models developed for materials other than hydrides. Models developed for MHs are, in particular, the Sun and Deng model (Section 3.1.6) [46,47], the extended Zehner–Bauer–Schlünder model (Section 3.1.7) [95], the improved area-contact model (Section 3.1.11) [50], and the Abdin–Webb–Gray model (Section 3.1.12) [102].

As k_s, k_g, and ε are standard parameters for all the models, some brief considerations are needed. Even if k_eff is not directly dependent on P, T, or X in several cases, these dependencies appear in k_s, k_g, and ε. In a metal hydride system, it can be stated that there is no dependence of k_s on P, as well as of k_g on X. Porosity mainly depends on X, even if a weak dependence on P and T is present, due to particle compaction and thermal expansion, but these parameters are usually not included in the models for ε. Similarly, the hydrogen mean free path is another parameter appearing in several models. When not differently indicated, it is calculated through the well-known relation described in Equation (134) [167]:

$$l={\displaystyle \frac{k_{B}T}{\sqrt{2}\pi P{d}_{kin}^2}}$$

(134)

In this equation, the kinetic diameter, d_kin, appears. In addition, a dependence on P and T is expressed, which is introduced in the k_eff models where this expression for l is included. In the following analysis, these assumptions are made:

• The k_s, k_g, ε, and l expressions given by the authors of every k_eff model are used;
• If the authors did not include any model for k_s, a constant value is used for this parameter;
• If no information or model for k_g is given, as it can be considered readily available and generally acceptable, it is calculated through data from RefProp v.10.0 (which is available from the free Mini-RefProp v.10.0 too) [168], where a dependence on P and T is present;
• If the authors did not include any model for ε, a constant value is used for this parameter;
• If no information on l is given, it is calculated through Equation (134);
• To analyze the influence of the composition on k_eff in the improved area-contact model, X, through R, is treated as a variable and not as a parameter.

The comparison was made using the parameters referring to the LaNi₅ alloy, as most of the solid material parameters appearing in the models for this alloy are available in the literature (

Table 1. Standard parameters for models’ comparison based on LaNi₅.

Parameter	Value
ε	0.52
μg	8.8 × 10−6 Pa·s
cp,g	14,300 J kg−1·K−1
γ	1.4
ks	8 W m−1·K−1
d	15 × 10−6 m
dkin	289 × 10−12 m
lv	0.034·d
e	0.5
Xmax	1.45 wt %
E	155 × 109 Pa

). The remaining unknown parameters were evaluated by fitting the models’ equations to the experimental data from Yoshida et al. [166] using the Curve Fitter tool from Matlab v.R2023B. These data were chosen because experimental data of k_eff vs. P, T, or X for LaNi₅ were needed for fitting, and while some examples are available for k_eff vs. P or X, there are no experimental k_eff vs. T graphs reported in the literature for LaNi₅, to the best of our knowledge. In addition, experimental conditions used to obtain these data are comparable to those at which material parameters are available in the literature. However, measurements of k_eff vs. P, T, or X for other materials have been reported, such as MmNi_x-yAl_y [63], HWT 5800 (Ti_0.98Zr_0.02V_0.43Fe_0.09Cr_0.05Mn_1.5), and LaNi_x-yAl_y [48,100].

The dependence on P, T, and X are analyzed from 0.001 to 100 bar, from −20 to 100 °C, and from 0 to 1.45 wt%, respectively.

In the following, the graphs of the simulated curves are reported. To express the dependence of k_eff on one variable, the other two were fixed: the fixed values for P, T, and X are 10 bar, 20 °C, and 1.0 wt%, respectively. The Abdin–Webb–Gray model (Section 3.1.11) [102] is not present because, using the standard parameters listed in Table 1, the argument of the erfc⁻¹ function, present in the expressions to evaluate the parameters a₂ and a₃, assumes values out of the domain for this function (that is, 0–2). This makes it impossible to evaluate these parameters, so k_eff. The simplified Hayashi [43] and the simplified Sun and Deng [46,47] models are not reported either, as the results are almost identical to those of the correspondent base models in the analyzed ranges (as it was said, the contribution of radiation is relevant only at high temperatures).

To obtain a more straightforward and contracted way to refer to every model, the following abbreviations obtained by the initials of the authors’ name or the acronym of the models’ name are used:

• 3.1.1 Maxwell model [42]: M;
• 3.1.2 Yagi and Kunii model [92]: YK;
• 3.1.3 Zehner–Schlünder model [44]: ZS;
• 3.1.4 Zehner–Bauer–Schlünder model [95,96,97]: ZBS;
• 3.1.5 Hayashi model [43]: H;
• 3.1.6 Sun and Deng model [46,47]: SD;
• 3.1.7 Extended Zehner–Bauer–Schlünder model [95]: EZBS’
• 3.1.8 Area-contact model [45]: AC;
• 3.1.9 Phase-symmetry model [45]: PS;
• 3.1.10 Raghavan–Martin model [49]: RM;
• 3.1.11 Improved area-contact model [50]: IAC;
• 3.1.12 Abdin–Webb–Gray [102]: AWG;
• 3.1.13 Heat transfer concentrating model [91]: HTC.

4.1. Effect of the Pressure on the ETC

The experimental outcomes from LaNi₅ obtained from [166] and other hydride-forming alloys such as MmNi_x-yAl_y [63], HWT 5800 (Ti_0.98Zr_0.02V_0.43Fe_0.09Cr_0.05Mn_1.5), and LaNi_x-yAl_y [48,100] have shown S-shape curves. The gas pressure causes a marked effect on the ETC, mainly at intermediate pressures. At extremely low pressures, lower than 10⁻⁴ bar, the Kn number ≥ 10 (Kn = mean free path (l) over characteristic length (Ʌ): Kn = l/Ʌ), the ETC is independent of pressure, molecules fly without collisions throughout the pores, the heat is transported by heat conduction in the particle bed, and the effect of radiation is neglectable. Between 10⁻³ and 10 bar, i.e., 0.01 < Kn < 10, the collisions between molecules increase with pressure, and the exchange of energy among the particles is dependent on the filling gas and proportional to the number of gas molecules. Hence, the ETC becomes pressure-dependent. At higher pressures, over 10 bar and 0.01 ≥ Kn, the pores are saturated with gas, and the gas phase can be considered a continuum, with the ETC independent of the pressure [100].

Figure 24 compares the dependence of ETC on P calculated by applying the models and the experimental data for LaNi₅ obtained from [166]. Reporting this graph using a logarithmic scale for the pressure axis improves the visual clarity. It can be seen that only a few models show the typical S-shaped curve, which are the ZBS (Section 3.1.4), H (Section 3.1.5), SD (Section 3.1.6), EZBS (Section 3.1.7), and IAC (Section 3.1.11). Those models showing S-shape curves include the effect of the gas in the biphasic region and the hollow region, as shown in Figure 12. The other models do not consider such features. Thus, it is possible to note that the other models give almost constant values of ETC with P, slightly increasing at high pressures. The PS model (Section 3.1.9) shows values for k_eff around 2.8 W m⁻¹ K⁻¹, which is much higher than all the other models, while YK (Section 3.1.2) and AC (Section 3.1.8) are around 1.5 W m⁻¹ K⁻¹, barely higher than the other models.

4.2. Effect of the Concentration on the ETC

The experimental evidence with hydride-forming materials such as MmNi_x-yAl_y [63], HWT 5800 (Ti_0.98Zr_0.02V_0.43Fe_0.09Cr_0.05Mn_1.5), LaNi_x-yAl_y [48,100], and LaNi₅ [166] have shown a steep increase in the ETC (in the solid solution region), followed by slight increases as the amount of hydride phase is considerable. This experimental behavior of the ETC has been attributed to the changes from the free hydrogen molecule flow (high Kn numbers) to the pressure-dependent flow (transition regimen: 0.01 < Kn < 10) with minor changes in the concentration (solid–solution). Then, as the concentration increases, a region of continuous gas (continuous regimen: 0.01 ≥ Kn) independent of the pressure is formed, and the ETC slightly increases or remains practically constant. It is also important to mention that there are further effects, such as the gas–solid reaction, which overlaps with the pressure effect [100].

In Figure 25, the dependence of the ETC on the composition is compared. Only three models are reported since they express a dependence on X: the SD model (Section 3.1.6), the EZBS model (Section 3.1.7), and the IAC model (Section 3.1.11). In addition, the experimental data for LaNi₅ obtained from [166] are also reported. The variation is almost linear for the three models, but while it increases for SD and IAC, it decreases for EZBS, as seen from the experimental results [48,63,100,166], the pressure always increases as the concentration rises, and therefore, there is a marked effect of the pressure. In the case of the calculation with the SD, IAC, and EZBS models, the pressure remains constant, which might be the reason why it is not possible to see the steepness at low values of X. However, with the increase in X in the case of the SD and IAC models, the ETC continuously increases. Meanwhile, for the EZBS model, the ETC decreases following the opposite trend. The main differences between the SD and IAC models and the EZBS model are the treatment of the contact area between particles and the model for k_s. The EZBS model might simplify the treatment of the contact area and use an equation for the calculation of the solid thermal conductivity with a decreasing trend with X (k_s is increasing with X for the SD model and constant for the IAC model).

4.3. Effect of the Temperature on the ETC

To the best of our knowledge, there are no experimental data on the dependence of ETC on the temperature for LaNi₅. However, an experimental example of the dependence of the ETC on the temperature was presented in ref. [63] for MmNi_4.5Al_0.5. The authors claimed that it is difficult to maintain a constant pressure and concentration as the temperature changes. Hence, the effect of the pressure and composition cannot be decoupled completely at the time to measure the dependence of the ETC on the temperature. For MmNi_4.5Al_0.5 with a porosity of 0.43, it was observed that, at low pressures, the effect of temperature is negligible, and the concentration also remains constant (the radiation effects are negligible). As explained in Section 4.1, this effect can be related to the molecules flying without collisions throughout the pores (high Kn numbers). Over a certain pressure (in the ref. [63]: 0.015 bar), as the temperature increases, the concentration slightly decreases, causing decreases in the ETC. It is not possible to speculate about this behavior since it is in the transition range defined by the values of the Kn number (0.01 < Kn < 10), and, according to the effect of the concentration on ETC, the ETC should show a slightly increasing trend. It was concluded that the effect of the temperature on ETC is minimal, increasing slightly as the temperature rises [63].

Figure 26 reports the dependence of ETC on T for the models. The variation of the effective thermal conductivity is almost linear and increasing for all the models, as expected from the experimental outcomes. Nonetheless, SD (Section 3.1.6) and IAC (Section 3.1.11) are slightly decreasing. Regarding the values, similar to the comparison on the pressure dependence, the PS (Section 3.1.9) curve assumes higher values than all the others, while the YK (Section 3.1.2) and AC (Section 3.1.8) curves are just higher than the others. The calculations were conducted at 10 bar, which is the lower limit of the continuous gas phases where the ETC does not depend on the pressure (0.01 ≥ Kn). The concentration used for the calculations was set at 1.0 wt%, which is relatively high, and the ETC remains constant. It is possible to propose that the slight increase in the ETC with T for all the models, except SD (Section 3.1.6) and IAC (Section 3.1.11) models, is directly related to the increase in the k_g of hydrogen [168]. In the case of the SD (Section 3.1.6) and IAC (Section 3.1.11), instead, an equation for the thermal conductivity of hydrogen showing a decreasing trend with temperature was used, leading to an ETC that decreases with temperature.

5. Summary and Conclusions

This work provides a comprehensive overview of the most important methods to measure the ETC, thirteen models for the calculation of the ETC employed for MHBs in chronological order of publication, and the effect of the hydride formation on solid thermal conductivity. The main equations, needed data, and parameters for each model were analyzed, and their advantages and disadvantages were investigated. Moreover, the models were applied to a test hydride-forming alloy, LaNi₅, and the dependence of the ETC on the temperature, pressure, and concentration of hydrogen was evaluated.

The experimental methods used to measure thermal conductivity have been classified as steady-state and transient methods. In the case of the steady-state methods, the axial and radial methods need longer measurement times since temperature equilibrium conditions should be reached. To the best of our knowledge, steady-state methods are carried out using in-house developed devices. The setup design constrains the experimental temperature range, hydrogen pressure, and sample amount among the most relevant parameters.

The transient methods, on the contrary, are commercially available: the transient line source (TLS), transient plane source (TPS), modified transient plane source (MTPS), and the eponymous hot-wire. Laser-flash and transient plane methods are the most frequently used ones. The TPS method shows a better accuracy than the other ones, though the determination of the ETC at elevated temperatures and pressures and its applicability to the industry remains a challenge [41].

The development of models to calculate the ETC of hydride compounds has been a matter of extensive research and has been in constant development. Figure S1 shows the temporal evolution of the models and their main features. Moreover, in the ESI, Table S1 provides a summary of the models with the central equation, brief description, parameters, and advantages and disadvantages, and Table S2 shows experimental and model-calculated ETC values and parameters for different hydride-forming and non-forming materials under different atmospheres, temperatures, and pressure conditions. Figure S1 shows a model classification into those for regular spherical particles (gas-film modification models), those for regular spherical and non-spherical particles (particle modification models), and the models developed for MHBs, which involve gas film and particle modification models. As also seen, there are common parameters for all the models: k_g, k_s, and ε, specific parameters for the gas film modification models: e, γ, l, and d, and specific parameters for the particle modification models: B, C, α which are related to the form and contacting of the particles. The models developed for MHBs need additional parameters, such as V, ϕ_p, P_eq, R_P, X_max, E′ (AWG [102]), and σ_R (AWG [102]). The parameters used in the gas film modification models are also utilized in the particle modification models and the MHBs models. Therefore, the number of parameters needed, especially for MHBs, notably increases, thus also increasing the complexity. In this regard, the simple additive equation k_eff = k_s·(1 − ε) + k_g·ε is widely used to model the ETC of MHBs [169]. Despite its simplicity, it does not represent the complexity of the effective thermal conductivity of MHBs, limiting the evaluation of the effects of several parameters on the ETC and leading to overestimating the value of k_eff (quick calculations would give values around 5 W m⁻¹·K⁻¹). The comparison of data reported in Table S2 also gives information on how accurately the models predict the experimental values of the ETC, making it possible to analyze how the complexity of the models relates to their accuracy. Nevertheless, they must take into consideration the type of systems, in terms of materials, and the operative conditions, in terms of temperature and pressure, at which these evaluations have been performed.

Table S1 exhibits a clear differentiation between the parameters that can be measured or acquired from the literature and those that are hard to measure or used as fitting parameters of the models. In most models, intrinsic parameters related to the changes of the material upon hydrogenation/dehydrogenation, such as void volumes, the form and contact of the particles, and mechanical properties, are usually employed as fitting parameters since it is challenging to quantify them. In Table S2, it is observed that in several cases, such intrinsic parameters are not even reported. From the presented and employed models for the calculation of the ETC, there are just four models developed for MHBs: SD/SSD (Section 3.1.6), EZBS (Section 3.1.7), IAC (Section 3.1.11), and AWG (Section 3.1.12) [45,46,47,50,95,102]. The SD/SSD was developed as a gas film model, while the EZBS, IAC, and AWG were created as particle modification models. Comparing the ETC experimental results and calculated with the ETC model, some models such as the M [42], YK [92], ZBS [95,96,97], H [43], AC [45], PS [45], RM [49], and HTC [91] present deviations from 2% to 80%. The model with the best agreement, RM, was evaluated in quite a narrow range of ETC so it might contribute to the slight deviation of less than 2%. In the case of the models developed for MHBs, the deviation between the experimental and model values of ETC stays between 5% and 15%. The gas film SD model and particle modification EZBS and IAC models include the change in the concentration of hydrogen. However, the AWG model does not. The SD, EZBS, and IAC models reproduce the dependence of ETC well on pressure (S-shape curve), temperature (slight increase), and concentration (increase/decrease). From the analyzed models, the EZBS model is considered the most suitable one to calculate the ETC, considering the complexity of the hydrogenation/dehydrogenation processes since it considers different gas regions (far from the particles and near the surface of the particles), change in hydrogen concentration, morphological, and mechanical properties.

It is a common practice [170] to add heat conductivity enhancers, like expanded natural graphite (ENG), copper, etc., as a mixture with the powder active material. Commonly, the ETC models consider the solid phase to be made only of the hydride-forming alloy. In the future, it would be of great interest to develop novel models in order to account for a multiphase solid material for the representation of the ETC through models that are as realistic as possible.