Optimized Modeling and Design of a PCM-Enhanced H2 Storage

Thermal and mechanical energy storage is pivotal for the effective exploitation of renewable energy sources, thus fostering the transition to a sustainable economy. Hydrogen-based systems are among the most promising solutions for electrical energy storage. However, several technical and economic barriers (e.g., high costs, low energy and power density, advanced material requirements) still hinder the diffusion of such solutions. Similarly, the realization of latent heat storages through phase change materials is particularly attractive because it provides high energy density in addition to allowing for the storage of the heat of fusion at a (nearly) constant temperature. In this paper, we posit the challenge to couple a metal hydride H2 canister with a latent heat storage, in order to improve the overall power density and realize a passive control of the system temperature. A highly flexible numerical solver based on a hybrid Lattice Boltzmann Phase-Field (LB-PF) algorithm is developed to assist the design of the hybrid PCM-MH tank by studying the melting and solidification processes of paraffin-like materials. The present approach is used to model the storage of the heat released by the hydride during the H2 loading process in a phase change material (PCM). The results in terms of Nusselt numbers are used to design an enhanced metal-hydride storage for H2-based energy systems, relevant for a reliable and cost-effective “Hydrogen Economy”. The application of the developed numerical model to the case study demonstrates the feasibility of the posited design. Specifically, the phase change material application significantly increases the heat flux at the metal hydride surface, thus improving the overall system power density.


Introduction
The exploitation of renewable energies has attracted a growing number of efforts in recent years [1], ranging from the theoretical physical aspects to the engineering implementation of novel energy systems. However, efficient and economically convenient implementations require the presence of energy buffers, due to the aleatory nature of renewable energies, especially from solar and wind sources [2,3].
Among the most promising technologies, both for power generation and energy storage, hydrogen-based energy systems have shown interesting potential in recent years, due to their innumerable advantages for the environment, economy and energy availability and security for final users [3][4][5][6][7][8][9][10]. One of the key issues for the actual development of H 2 economy is the distribution and storage of such an element [11]. In this respect, metal hydride (MH) H 2 storage systems provide a safe, compact and reversible storage solution, even for large H 2 amounts: this is the reason why such a technology has increasingly attracted research efforts in recent years [12,13]. More specifically, many reversible H 2 storage devices exploit light intermetallic hydrides, such as AB5 and AB2, reaching hydrogen storage capacities of about 2 wt% [14]. Nonetheless, the hydride formation/hydrogen absorption is an exothermic reaction that generates a quantity of heat that is approximately equal to the absolute value of the reaction enthalpy [15]: the generated heat decreases H 2 absorption rate, reducing the kinetics of the hydride formation. Thus, system performances are largely affected by thermal management of the whole storage [16].
At present, technical and economic feasibility of hydrogen storage systems are still far from being achieved: in fact, none of the current metal hydride technologies fulfill all the essential criteria for a practical Hydrogen Economy, due to low storage capacity, slow kinetics and difficulty in managing hydrogen absorption/desorption temperatures [11].
Thus, the optimal thermal management of Metal Hydride storage systems provides the framework of the present work: since hydrogen absorption and desorption are exothermic and endothermic processes, respectively, an optimized thermal control and management during H 2 storage and release dramatically increases the storage performance and power density [17].
In this work, we address the problem of enhancing the performance of a metal hydride tank for hydrogen storage, both for transportation applications [13] and micro-CHP (Combined Heat and Power) implementations. More specifically, we first show that by coupling the PCM-based storage to a metal hydride tank, the power density can be dramatically increased by storing/releasing the heat connected to the absoprion/desorption of H 2 from the hydride. Then, we develop a fast and reliable numerical method to evaluate the performance of the PCM storage by means of the Lattice Boltzmann Method [18]. LBM has been applied to a vast panorama of problems of scientific and technological interest in recent years [19][20][21][22][23][24], providing a reliable and accurate alternative to traditional solvers based on the continuum assumption, such as the Navier-Stokes approach. The results provide a promising framework for a new generation of safe and sustainable energy storage systems to promote the diffusion of an effective hydrogen economy.
However, controlling the heat transfer from/to the canisters is still a relevant challenge for such systems [16]. In fact, H 2 absorption is an exothermic reaction in which the enthalpy variation ∆h mh is in the range [10 kJ/(kmol)-30 kJ/(kmol)] and the thermal power that can be removed/supplied to the canisters limits the overall hydrogen mass flow [15,16]: such a thermal flux is limited both by the heat conduction inside the hydride and convection from the canister to environment [15,16].
In this paper, we tackle the heat convection limitation. More specifically, we posit the challenge to enhance the performance of a metal hydride-based H 2 storage through its hybridization with a thermal energy storage based on a phase change material (PCM-TES), in order to increase the thermal power exchanged by the canister surface and to enhance the temperature and equilibrium pressure control [33]. We leverage the design of an existing canister [27] to evaluate the feasibility of the hybrid PCM-hydride storage and to predict the gain in system performance in terms of H 2 flux increment.

Hybrid Metal Hydride PCM-TES Design
We consider the off-the-shelf canister MyH2-2000 produced by H2Planet [27], whose relevant properties are reported in Table 1, as the starting point to design the hybrid MH-PCM storage system (see Figure 1a).  To model the heat generated by the H 2 absorption process, or equivalently required by its desorption, we select from [15] an AB 2 type alloy that has an equilibrium pressure in the range reported in Table 1.
Specifically, such an alloy is characterized by ∆h mh = −14 × 10 3 J/mol and by an entropy variation ∆s mh = −64 J/(mol K), resulting in an equilibrium pressure of 12 bar at a temperature of 323 K (see Equation (1) [15] where p 0 = 1.0 bar is the reference pressure, T is the temperature, and R u = 8314 J/(kmol K) is the universal gas constant).
The H 2 gravimetric density of the storage is w %,H 2 = M H 2 /M mh = 1.17 % where M H 2 and M mh are the masses of the H 2 and of the canister, respectively. A complete absorption cycle releases a quantity Q = M H 2 ∆h mh /M = 1155 kJ , of heat, where M = 2 g/mol is the molecular weight of H 2 . Equivalently, 58 W of thermal power must be drained from the canister surface to allow an H 2 flux of 1 kW while keeping its temperature constant. The same power should be provided during the desorption phase.
The heat exchange process at the canister surface is summarized by the Nusselt Number, Nu = λl/k [34], where λ, is the convection coefficient, l is the reference length, and k is the thermal conductivity of the surrounding fluid.
For the baseline design, we proxy the canister with a vertical cylinder immersed in quiescent air. Thus, k = 0.0256 W/(m K), l = 0.54 m is the height of the canister and, for natural convection [34]: Nu = 0.0210Ra 2/5 for Ra > 10 9 (2) where Ra is the Rayleigh number, g = 9.81 m/s 2 the magnitude of the gravitational acceleration, β, is the volumetric expansion coefficient, T env = 293 K is the temperature of the surrounding fluid, ν = 1.50 × 10 −5 m 2 /s its kinematic viscosity, and Pr = 0.71 the Prandtl Number. We schematized the canister as a lumped parameter system with constant temperature at the outer wall, T = 323 K. Table 2a summarizes the resulting mail physical parameters relative to the heat exchange at the canister surface. Specifically, the thermal power exchanged at the surface of the baseline canister is limited toQ surf = 22 W. In turn, the H 2 mass flow is limited to 3.17 × 10 −6 kg/s (see Equation (3)) that are equivalent to a storage power ofẆ =ṁ H 2 LHV H 2 = 308 W, where LHV H 2 is the Lower Heating Value of hydrogen.ṁ In such conditions, about 14 h are required to completely fill or unload the reservoir. We finally evaluate that the Biot number, Bi = λk mh /l * , where l * is a reference length equal to the ratio between the canister volume and the heat-exchanging surface. According to [34], the conductive resistance can be discarded, and the solid can be assumed at constant temperature for Bi < 0.1: for the baseline design, Bi < 0.1 if k mh > 1 W/(m K). Enhanced configurations have higher conductivity [15,[30][31][32] but k mh 1 W/(m K) is also feasible for pure MH canisters [28,29]. In the following, the lumped parameter hypothesis holds, and λ is effectively limiting the input/output power of the storage.
To increase the heat exchanged at the surface, we improve the design of the baseline storage by flooding the canister within a PCM thermostatic bath (see Figure 1b). As PCM, we consider the paraffin octadecane, whose relevant thermo-chemical properties are reported in Table 2b.  In the following, we focus on the melting process. Heat transfer mechanisms at the solid wall follow three regimes, namely pure conduction, mixed conduction/convection, and pure convection. In two dimensions, we proxy the whole system as a rectangular cavity with one hot wall (at the inner) and a cold wall (at the outer, see Figure 1b). For this geometry, the different heat transfer regimes can be summarized through Equaton (4) [36]: where θ = Fo St is the non-dimensional time, St = C (T − T env )/L is the Stefan and Fo = αt/l 2 the Fourier numbers, respectively, and c 1 , c 2 , and n are empirical parameters to be evaluated according to the melting regime. Before the upper boundary of the domain is completely melted c 1 = 0.35, c 2 = 0.0175, and n = −2 [36]. For the convective regime, Equation (4) asymptotically approaches Nu ∞ calculated through Equation (5) (see Figure 2).
For the MH-PCM hybrid system Ra = 8.90 × 10 10 (see Table 2b) and Nu > Nu ∞ for θ > 10 −4 (see Figure 2). Such a threshold corresponds to t > 0.2 s, which is much lower compared to the expected time to load the canister. Thus, we can safely assume that Nu = Nu ∞ for the whole process. The results summarized in Table 2b highlight that the PCM significantly enhances Nu, λ and, in turn,Q surf . The thermal power exchanged at the surface is about 10 times larger compared to the baseline design. This increment would allow an H 2 mass flow oḟ m H 2 = 2.97 × 10 −5 kg/s, corresponding toẆ = 3565 W. With such a mass flow, the canister is completely charged/discharged in about 1.5 h.
However, we wish to emphasize that, in this case, Bi < 0.1 only for k mh > 15 W/m K, which is feasible only for advanced canister design with enhanced inner thermal conductivity [15,[30][31][32]. For pure MH alloys Bi > 0.1, the lumped parameter hypothesis is not strictly valid, and thermal conduction inside the canister might also limit the storage power.
The total heat released during the charge process is Q = m H 2 ∆h mh = 1155 kJ. The amount of PCM required to store such thermal energy is m PCM = Q/L = 4.73 kg, corresponding to a volume of V PMC = 5.92 × 10 −3 m 3 that can be obtained through a jacket with external diameter D = 0.16 m (see Figure 1).
The results reported in Table 2 provide evidence of the feasibility of the hybrid energy storage concept. However, more thorough analysis of the heat exchange processes at the canister surface is required for the design of a working prototype, especially to evaluate the system performance during transient operation.
To this aim, in Section 4, we propose a novel, accurate, and fast method to predict the local values of Ra and Nu, in order to convey more thorough information for the design of heat-recovery systems for metal-hydride-based H 2 storage solutions. Specifically, leveraging the unsteady nature of the Lattice Boltzmann Method (LBM) and its mesoscopic foundations, we propose a multi-dimensional method to track the evolution of PCM melting and solidification, accounting for the energy storage and release, as a function of Nu in the different regimes of the phase change [36].

Multidimensional Modeling of PCM Melting & Solidification
We adopt a standard 2-population LB implementation [18,37,38] within the Lattice Boltzmann Method (LBM) to solve the fluid dynamics and the thermal evolution: In the above, we have τ F the (standard) fluid relaxation time, given by τ F = ν c 2 s + 1 2 (with ν the kinematic viscosity); τ T the thermal relaxation time, provided by τ T = α c 2 s + 1 2 , with α representing the thermal diffusivity.
Equation (6) accounts for the fluid dynamic evolution of the fluid, while Equation (7) is responsible for the thermal evolution. With this approach, the temperature is treated as a passive scalar, whose effects are felt by the fluid according to ad-hoc tuned forces (i.e., buoyancy effects).
The three forcing terms at the rhs of Equations (6) and (7) represent, respectively: F B the buoyancy effect; F Drag the fluid-solid interaction, given by an empirical mesoscopic drag, which enforces no-flow inside the solid region [39]; finally, F L stands for the latent heat release/absorption during melting/solidification phenomena, respectively.
To compute these forces, one has to evaluate the evolution in time of the phase field φ that governs the melting/solidification front. More specifically, we have the following rate equation, in the discretized approach: in which, R is the reaction term expressed as in [39], in which f + and f − are frequency factors (basically, the inverse time scale for melting and solidification, respectively) and K + and K − are switch functions controlling the onset of melting and solidification, respectively, around a critical temperature T cr . Here, we directly adopt the modified switches proposed in [39], given by In the above, T w controls the energy range of the transition and T s provides two different activation energies for melting and solidification. In our basic test cases, T s = 0.
From the evolution in time of φ, we find the three force contributions in Equations (6) and (7): In the above, we have • is a relaxation parameter, chosen to be in the order of 0.5, to ensure that the relaxation of F Drag is faster than the other dynamics; • St is the non-dimensional Stefan number, given by St = c , with c the specific heat and L the latent heat. • ∆T is the characteristic temperature difference: we have fixed ∆T = T liq − T sol ;

Numerical Results & Discussion
In the following, we first validate the proposed numerical method against the classicl 1D Stefan Problem, in which the buoyancy effects are neglected. Then, we focus on a 2D layout, enabling the effects of gravity and buoyancy and validate our code against results acquired from the literature.

1D Stefan Problem
We take as benchmark the results in [40], considering three values for thermal diffusivity, α = 0.033, 0.36 and 0.70 in lattice units. To reproduce these results, we have considered the very same Stefan number as in [40], St = 1: fixing T sol = 273.0 K, we have T liq = 352.5 K.
The domain is a 100 × 5, 2D box, as in [40]. Fist, we investigate the type of boundary conditions that must be implemented to accurately capture the evolution of the 1D melting problem. In our double-population approach, we have to dissect the nature of the boundary conditions according to the population we are considering: from the hydrodynamic point of view (i.e., for the f population in Equation (6)), our computational domain is a closed box with wall-no slip conditions on all boundaries. For what concerns population g in Equation (7), that is, for the thermal evolution of our system, we have a different picture: the left boundary (West) is the heated side, thus it has a constant temperature, higher than that of the solid PCM; the right boundary (East) has a fixed temperature, equal to that of the solid PCM; the top and bottom sides of the domain (the North and South) are schematized as adiabatic walls.
To implement such adiabatic conditions, at least two methods can be retrieved from the literature [37,41]: a (thermal) free-slip wall and an equilibrium condition. While the first one can be immediately recognized as a strightforward adiabatic condition, the second one can provide accurate results, provided that the thermal gradients (both in time and space) are small compared to the characteristic space and time scales of the problem under investigation [41].
The results of the 1D simulation for the two described boundary conditions are reported in Figure 3, which shows the comparison with the results reported in [40]: Figure 3a reports the results obtained with the free-slip boundary condition (bc) is set at the North and South walls, while in Figure Figure 3 highlights that no appreciable differences emerge when implementing equilibrium or free-slip bc's in the 1D cases, due to the absence of gravity and buoyancy.

2D Melting: The Effects of Gravity and Buoyancy
In this Section, we consider a square cavity filled with solid PCM. As in the previous simulations, the left wall is heated (i.e., it is at constant temperature T wall = T liq ), the solid is kept at T sol = T melt . Our 2D computational domain reproduces the portion of the PCM jacket represented by the dotted rectangle in Figure 1b: the left wall is the heated side, due to the exothermic reaction of H 2 absorption within the hydride; the right wall is schematized at constant temperature and the top and bottom sides are assumed as adiabatic. Considering a small portion of the entire canister allows us to spare considerable CPU time, but the local Ra and Nu numbers will be lower than those of the case analyzed in the previous sections. However, such values will be compliant to those of the benchmark cases reported in [40,42]. Figure 4 reports the evolution in time of the Nu number computed at the heated wall, as a function of the Ra number. The figure highlights that our model is capable of following the transient evolution in Nu, due to the different regimes associated with the melting of the soldi PCM [36,42]: first, a conduction-dominated phase, with a sudden drop in the Nu number; then, Nu re-increases and reaches a plateau in the conduction-dominated phase; finally, a progressive reduction in Nu, that starts when the melted front reaches the right (cold) wall. Moreover, compared to the numerical benchmarks proposed in [40], our model provides values of the Nu numbers in the plateau region that are closer to the predictions provided by the equations developed in [42].   Figure 5, the evolution of the temperature field in time in our computational domain, compared to the contours provided in [40]: the panel highlights the good matching between our results and those in [40] for the same Pr, Ra and St numbers. Visto lo Statuto dell'Università degli Studi di Roma "Tor Vergata" Visto l'art. 4 della legge 3/7/1998, n. 210, che prevede che le Università con proprio regolamento disciplinino l'istituzione dei corsi di dottorato, le modalità di accesso e di conseguimento del titolo, gli obiettivi formativi ed il relativo programma di studi, la durata, il contributo per l'accesso e la frequenza, le modalità di conferimento e l'importo delle borse di studio, nonché le convenzioni con soggetti pubblici e privati, in conformità ai criteri generali e ai requisiti di idoneità delle sedi determinati con decreto del Ministro; Visto l'art. 19 Figure 5 highlights the very good matching between our results and those provided in [40]. The time evolution of the interface between solid and liquid paraffin, in fact, follows the very same steps in the two models, at Ra = 6.8 × 10 6 : • θ ∼ 0.01: the end of the conduction-dominated phase, with the formation of a rounded sack on top of the computational domain, due to the inception of convective motions; • θ > 0.05: advance of the interface towards the cold wall, with the formation of a rounded profile in the upper part of the computational domain, due to the melting front pointing downwards; • θ > 0.15: the melting front reaches the cold wall: this corresponds to the end of the plateau region of Nu in Figure 4; • θ > 0.2: the progressive shrinkage of the solid phase.

Conclusions
Hydrogen storage solutions will play a pivotal role, in the coming years, in implementing a reliable and cost-effective "Hydrogen Economy". In this framework, metal hydrides have proven to have a remarkable potential in delivering safe and reliable H 2 storage systems. In this paper, we have provided novel, promising results for the implementation of metal hydride H 2 storages, by coupling them with a PCM-based thermal energy storage. More specifically, we have designed an engineering solution to store the heat released by the hydrides during the H 2 loading process by means of a Phase Change Material (PCM). Specifically, preliminary design results demonstrate the high potential of the proposed solution in terms of increased power density. In fact, the PCM allows for better control of the MH temperature, and, in turn equilibrium pressure. Moreover, we provide evidence that such a solution would have a negligible impact on the system cost, mass, and volume. However, in order to reach a more accurate system design capable to address the local (transient) heat transfer phenomena during the melting/solidification phases, we have developed an innovative LBM implementation to assist the engineering design of optimized PCM-MH storage systems.
More specifically, to account for the local evolution in time of the Nu number, we have developed a fast and reliable 2D algorithm to track the melting of the paraffine-like material chosen as PCM. Results reported in Section 5 demonstrate the the proposed methodology The model provides an accurate and reliable description of the melting phenomenon, compared to the results provided in the literature.
Further developments hinge on the promising results evidenced both by the innovative numerical LBM implementation and by its application to the MH-PCH hybrid system. Specifically, further work should be conducted to account for density variations within the proposed numerical framework. Similarly, the proposed methodology will be used to numerically model MH-PCM system.

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

1D
One