Modelling a kinetic deviation of the magnesium hydrogenation reaction at conditions close to equilibrium

A model has been derived for the magnesium hydrogenation reaction at conditions close to equilibrium. The reaction mechanism involves an adsorption element, where the model is an extension of the Langmuir adsorption model. The concept of site availability ( s s ) is introduced, whereby it has the capability to reduce the reaction rate. To improve representation of s s , an adaptable semi-empirical equation has been developed. Supple-ment to the surface reaction, a rate equation has been derived considering resistance effects. It was found that close to equilibrium, surface resistance dominated the reaction.


Introduction
Research efforts to store surplus energy for concentrated solar power (CSP) plants have looked beyond latent heat technologies, to a potentially cheaper and more effective technology, thermochemical energy storage (TCES). [1,2], A potential TCES system is based on utilising the enthalpy of reaction of a reversible metal hydride reaction, where metal hydrides are selected based on several governing factors, such as operating temperature, enthalpy of reaction, plateau pressure and cost. Many hydrides have been identified for a TCES system, in part due to their versatility, for example combining magnesium with a transition metal such as nickel or iron enables the tuning of the operating temperatures at a reasonable cost [3]. In addition, adding a destabilisation agent, such as In or Ge, enables thermodynamic modification for further tuning within a TCES system if needed [4,5]. However, this paper focuses on magnesium hydride, due to the simple case of a binary system and the suitability for a thermal oil circuit [6]. A key challenge with the design of this system are the reaction kinetics, where sufficient understanding is needed to facilitate an efficient energy store [2]. As such, this paper focuses on the kinetics, specifically the hydrogenation reaction.
Analysis of the literature reveals the intricate nature of the Mg hydrogenation reaction, whereby the reaction is sensitive to temperature, particle size, gas pressure, irregularities in particle structure, concentration, and diffusion effects [7e9] A popular choice are the nucleation-growth-impingement models known as the JohnsoneMehleAvramieKolomogorov (JMAK) [10]. This model is based on the concept of an extended volume. However, there are drawbacks with this concept, as in practice the reactive surface area would decrease as more reactant depletes.
Rate equations developed for metal hydrogenation that include temperature dependency are commonly governed by an Arrhenius term, and this is widely accepted, however the pressure dependency is still being debated. These rate equations can be written in many forms, such as dimensionless [11], mass concentration [12] or mole concentration [13], of which they commonly resemble a first order equation based on the original work of Mayer and Groll [14]. Chaise et al. also attempted to fit their data using an Arrhenius term and other reaction mechanisms [15].
A pressure term typically contains the gas pressure P and equilibrium pressure P eq . These range from a straight pressure driving force (conventional Langmuir) [16,17], or via a pressure ratio. It has been observed that modelling P eq near to the equilibrium isotherm yields a closer result to experimental data [13]. To model the equilibrium pressure, approaches have varied from entirely empirical equations, such as a partition of two separate exponential functions [13], to semi-empirical functions which model the overall shape with mathematical functions and use the Van't Hoff equation to determine the temperature dependence [12].
Switching attention to the overall format of the pressure ratio, Mayer et al. modelled a LaNi 4$7 Al 0.3 bed using the term ln (P/P eq ) [14]. Herbrig et al. modelled a hydride-graphite composite bed using (P-P eq )/P eq [13]. This term originated from Ron (1999), who experimented with different pressure terms [18]. A pressing question is how the pressure dependency is interpreted, whether that be by conventional Langmuir assumptions or perhaps via a new mechanistic concept, or, not at all.
As such, this paper presents the findings of an alternative derived model that encompasses an expansion of Langmuir, temperature & pressure dependency, and diffusion & surface effects. In addition, this model is compared to experimental data.

Site availability model (SAM)
The SAM contains elements of the Langmuir's adsorption model (LAM) [19], and a shrinking core approach [20], whereby the overall reaction progresses through steps, shown below. Fig. 1 illustrates the mechanism.
The simple case whereby step 3 is assumed the ratedetermining step is shown in Adsorption reaction. Inclusion of Step 1 is discussed in section resistances. The SAM is compared with LAM below:

Site de-activation
For H 2 sorption to a metal to occur, the gas pressure (P) must be greater than the equilibrium pressure (P eq ). As the exothermic reaction invokes a rise in temperature, which causes the equilibrium pressure P eq to rise and approach P; the thermodynamic driving force is reduced and so must the reaction rate. In effect, P eq can be considered the thermodynamic boundary.
Expanding on site theory, if the reaction rate must fall, the rate of a successful collision to an active site must decrease. Thus, it is assumed that the rise in temperature and subsequently the reduction in pressure driving force, must temporarily de-activate a site, even if that site is unoccupied/ unavailable, to reduce the amount of successful collisions, and therefore reduce the reaction rate.
Within SAM, this theory is represented by the ratio of available sites to unavailable sites. As temperature is linked to the equilibrium pressure, then across each time step, P eq indicates the unavailable sites. At the next time step, the change in (P-P eq ) will be an indication of sites that have become inactive across that time step. Therefore, the ratio that sites 1 Transport of H 2 to Mg surface a. Initial permeation through the porous bed (external diffusion) b. Diffusion through product layer (internal diffusion) 2 Dissociation of H 2 at Mg surface: 3 Reaction via adsorption and depletion of core: The rate of attachment to an active site (S) is proportional to rate of collisions that H make with active sites Same An active site is an unoccupied site. An active site is a site that is unoccupied and available. The collision rate to an active site is proportional to the pressure over the surface (P H ).
The collision rate to an active site is proportional to the site availability (s S ).
As H only adsorbs on active sites, it is proportional to concentration of active/ vacant sites. (C v ) Same r H,S ¼ k"P H C v r H,S ¼ k"s s C v i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 9 1 2 3 e2 9 1 3 1 that are available to ones that are unavailable can be estimated by the site availability, s s . s s z S av S Unav z P À P eq P eq (1)

Isotherm representation
A semi-empirical isotherm equation has been developed to represent metal hydrogenation systems. The equilibrium pressure is calculated by Eq. (2).
P N is the normalisation pressure. For a typical PCT, this point would be where q ¼ 0.5, often this position is the plateau pressure. Note, Eq. (2) does not work at q ¼ 0 or at q ¼ 1. The slope parameter y 1 can be constant or a variable, such as f(q), where increasing y 1 increases the slope. Eq. (2) can be modified for a system that exhibits two separate growth phases: where 'y 2 ' is a constant that can shift the location of the normalisation point.

Adsorption reaction
If step 3 is considered the rate-determining step, the rate of forming H,S is, Assuming a site balance where the total sites is the sum of vacant and filled sites (C T ¼ C v þ C H$S ), and incorporating the normalised hydride fraction q ¼ C H$S =C T , gives a first order rate equation, This is an ideal scenario, where the sample can reach full capacity. If C T is, where C ref is the reference concentration (taken as 7.66 wt%) and x m is the maximum capacity fraction at operating conditions. Including the hydride fraction (a helpful performance measurement), x ¼ qx m the rate equation can be modified,

Resistances
The model now includes surface and diffusion resistance. The SAM neglects the bulk movement of hydrogen through the bed; assumes a constant bed porosity (ε b ) and the expansion of the metal hydride forms surface cracks, grain boundaries and defects in the structure allowing the permeation of molecular hydrogen direct to the inner metal surface. This is represented by Knudsen diffusion [21]. If the reaction of H þ S/H,S is considered, where the rate equation is based on gaseous atomic hydrogen which takes part in the reaction, the surface resistance balance for a spherical pellet is, Where C Hc is the concentration of H at the core interface. Converted into effective volume terms and rearranging: Fig. 1 e Proposed general mechanism of hydrogenation kinetics for Magnesium. This is a 2D cutout section of a sphere. R 0 is the radius of the sphere, R c is the radius of shrinking metal reactant. i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 9 1 2 3 e2 9 1 3 1 where k e ¼ k '' RSA and RSA is the reactive surface area. If RSA is based on a sphere, The radius of the core (R c ) is related to the constant particle size (R 0 ) by: For the diffusion resistance balance: Applying , the pseudo steady-state assumption and integrating [20].
Combining the resistances, (Eqs. 10 & 14) gives, The surface resistance is represented by R 2 0 RSA=R 2 c k e s s and the diffusion resistance is R 0 =D e R c ðR 0 À R c Þ. Eq. (8) is based on the concentration of atomic hydrogen in the solid state, whereas Eq. (15) is based on the concentration of atomic hydrogen in the gaseous phase, which takes part in the reaction. E.g. if the initial gas pressure was 20 bar and the final pressure was 18 bar, Eq. (15) is based on the moles corresponding to 2 bar. Therefore, a mole balance can be used where if the hydrogen is not bound, it must be in the gas phase. This enables calculation of the reacted fraction.

Experiment
Experiments were performed using a Sievert's apparatus, where high purity hydrogen (99.9999%) was used. The sample was atomised Mg powder (26 mm) by SFM -FluorsidGroup Company. The bulk density was 1000 kg/m 3 , thus giving an assumed porosity of (1e1000/1740) z 0.4. Kinetic runs were at several constant initial temperatures and when the initial overpressure ([P-P Pl ]/P Pl ) ¼ 1 (P Pl is the plateau pressure). We assume that when ([P-P Pl ]/P Pl ) ¼ 0, it is at equilibrium, when the overpressure is 1, it is close to equilibrium, and when > 1, we are far from equilibrium for the temperature range 330-360 C.
The equivalent volume method was used [22], with a manifold volume of 62.9 cm 3 and an effective sample cell volume of 11e13 cm 3 depending on the temperature. The sample temperature was measured with a thermocouple in the sample.
Each sample was activated at 360 C and 40 bara over 4 cycles. These cycles were typically 2e3 days each. For the first sample, the capacity successfully recovered by the fourth cycle and underwent experiments in triplicates at 330 C and 360 C, where the total number of cycles was kept low at approximately 10 cycles, in order to minimise the possibility of sintering. To further minimise sintering, both samples were always kept in the hydrogenated state overnight, and once fully dehydrogenated (all dehydrogenation was at 360 C and took approximately 1 h), hydrogenated again within an hour. The second sample was activated under identical conditions as the first sample, whereby capacity was recovered by the fourth cycle. Experiments of sample two at 345 C were also performed in triplicates.

Geometry
The sample sat at the base of the sample holder hole where a thermocouple was placed 15 mm from the edge of the sample. It was observed to remain relatively steady at the set point. Thus, a constant temperature boundary condition was assumed. Fig. 2 illustrates this.

Sample dimensions
L 2 was assumed to be L 2 ¼ 0.2L 1 . The sample mass was 0.2 g, and assuming the magnesium's maximum capacity, wt m ¼ 0:0766 the reference mass of hydrogen at reaction completion (m ref ) is given by To account for the changing density, an assumed average porous matrix density (r p ) was used. A simple average density of Mg and MgH 2 gave a r p ¼ 1590 kg/m 3 . This created a pseudo value for L 1 and L 2 , equating to an effective sample domain volume (V de ). V de is calculated by considering the average maximum uptake of hydrogen fraction (constant) at operating conditions x av m (i.e. the experimental end capacity).

Balance equations
To model the experimental data, the finite element method was used and achieved by using COMSOL Multiphysics 5.3a.

Mole balance
External diffusion effects are neglected, and the balance is based on the gas phase H concentration taking part in the reaction considering the bed void volume.
The Arrhenius equation calculated the effective specific rate constant k e , with activation energy E e and frequency factor A e .
As the maximum capacity (x m ) varies with temperature, an empirical relation was used, with a 1 and a 2 being constants. The equation for x m changes for each operating condition. i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 9 1 2 3 e2 9 1 3 1

Energy balance
Conductive heat transfer was expected to dominate, thus advective heat transfer was neglected and there was thermal equilibrium between the porous matrix and hydrogen [12].
With DH being the enthalpy of reaction. Considering a porous matrix and gas domain, the volumetric heat capacity ðrc p Þ e and thermal conductivity l e were expressed in effective terms.
where c pH 2 and l H 2 were governed by an empirical equation which was a f(T,P), using data from NIST database at the system pressure ranges [23].
The specific heat capacity (c pp ) and porous matrix thermal conductivity ðl p Þ were assumed to be constant. l p was selected to reach the desired value of l e . At 300 K and 1 bar, the measured l e of MgH 2 under air and argon was 0.116 and 0.088 W m À1 K À1 respectively [24]. Together with the gas thermal conductivities, linear interpolation at the same conditions (but in H 2 ) resulted in an $ l e ¼ 0.7 W m À1 K À1 . Even though l H2 is higher at operating conditions, this conservative estimate of l e was used.

Isotherm modelling
The calculated polynomial for y 1 from experimental magnesium sorption PCTs between 300 and 390 C [25] and used within COMSOL simulations is shown in Eq. (26). The paper calculated a DH ¼ À76.07 kJ (mol H 2 ) À1 ± 1.21 and a DS ¼ À137.89 JK À1 (mol H2) À1 ± 1.97. The values of DH & DS calculated by these authors are comparable to Refs. [26,27]. This equation sufficiently represents the plateau and b phase but does not include the a phase. This is because the total a-phase reaction time was less than 10 s, and thus could be omitted. The PCT is shown in Fig. 3, where y 1 was regressed using data at 300 C and assumed not to be a function of temperature. y 1 ¼ À245:09q 4 þ 310:95q 3 À 123:74q 2 À 12:268q þ 73:943 (26) Kinetics Using the equations outlined in previous sections, enabled the calculation of the hydride fraction as the reaction progressed, shown in Fig. 3. The calculated effective activation energy i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 9 1 2 3 e2 9 1 3 1 (E e ¼ 172 kJ mol À1 ; A ¼ 4:5 Â 10 11 s À1 ) was within the error of one standard deviation to the values determined by DSC dehydrogenation under argon via Kissinger plots of the same magnesium batch (169±9 kJ (mol H2) À1 & A ¼ 4:5 Â 10 11 s À1 (lower limit À1:8 Â 10 11 s À1 and upper limit þ 6:8 Â 10 11 ) [28].
To calculate E e , we started with the initial estimate of the experimental data. Then A e was kept constant with E e as a regression parameter within the COMSOL simulation. The reaction at 360 C has progressed further compared to 330 C and 345 C (of which no reaction reached completion) and results indicate that lower temperatures result in higher capacities at conditions close to equilibrium. Altering the yintercept of Eq. (20) is required with a change in initial temperature. The gas pressure experimental data was directly used within the rate equations and shown in Fig. 3 for reproducible purposes. However, for a typical Sievert's apparatus kinetic experiment, if one plots P versus q, it is a straight line. Therefore, it is possible to calculate the gas pressure alongside q as the reaction progresses with a suitable expression.

Diffusion & surface resistance
Eq. (15) enables the resistances to be plotted independently. In the first case, the diffusion resistance is assumed to have a product layer porosity of 0.3 (expansion by 30%) with an assumed pore radius of 1 nm. In case two, the diffusion resistance is maximised by minimising the porosity and pore radius. Here, it is assumed a practical minimum of porosity ¼ 0.1 and the pore radius ¼ 0.5 nm (z2x the diameter of 1 hydrogen molecule) to give a maximum Knudsen diffusion of 10 À8 m 2 /s. From Fig. 4, it is apparent that the diffusion resistance does not exceed 8000 s m À1 for when the effective diffusion coefficient is maximised.
In comparison, the surface reaction resistance is in the order of 10 7 . In effect, a small diffusion resistance implies that hydrogen diffusing through the hydride has little effect on limiting the kinetics, if based on Knudsen diffusion. However, one cannot eliminate the possibility that diffusion resistance becomes influential at far from equilibrium operating conditions, where the additional overpressure could change the process in how the hydride forms, minimising cracks and forming an encasing shell, i.e. no longer entirely Knudsen diffusion.

Surface resistance
Consequently, if it is assumed that the diffusion resistance is negligible, Eq. (15) simplifies to Eq. (27): Analysing Eq. (27) shows that the effective rate constant ðk e Þ dominates in comparison to ðR c =R 0 Þ 2 and s S . As k e ¼ k '' RSA where RSA ¼ 3R 2 c =R 3 o ; an increase in RSA would increase k e and thus raise the reaction rate, which is achieved by a reduction in particle size (R 0 ). This mathematically represents that reducing the particle size, is an effective method of increasing the reaction rate. If Eq. (27) is further simplified by substituting (12), this gives the result in terms of concentration only ((28)), and if expressed in terms of the normalised hydride fraction gives (29). This result indicates that if there is a sufficient reduction in the RSA, feasibly a high wt% material, where the rate determining step is a reaction occurring at a surface of a sphere at conditions close to equilibrium, the reaction order is 5/3. If in contrast, the material is low wt%, where there is not a significant reduction in R c , thus R c /R 0 z 1; Eq. (27) simplifies to a first order rate Eq. (6) e similar to kinetics currently used for low wt% materials based on the work of Mayer and Groll [14] and popularised by Jemni & Nasrallah [12]. In effect, a low wt% material potentially contains a low total amount of available sites for hydrogen relative to the total surface, so the metal reactant surface does not shrink. In contrast, a high wt% material contains a high total amount of available sites, resulting in a shrinking metal reactant and therefore, the analysis indicates inclusion of the reacting surface is important when deriving a suitable rate equation for a high wt% material.
In addition, if s s is assumed a constant, Eq. (28) is suitable for integration: However, there are issues with the integral method. Firstly, reaction conditions should strictly be isothermal, to reliably assign a single 'k' value at a given temperature. This can be problematic with metal hydride samples, due to low effective thermal conductivities. Further, if the inclusion of a pressure term is a valid assumption, then the site availability is not constant due to the variation in P and P eq . Thus, the accuracy of the calculated effective activation energy is mainly dependent on the estimated value of s s . Even if partial differentiation was applied to Eq. (28), regarding s s , an estimation over time would still be required. In effect, the analysis indicates the integral method is not recommended to completely model H 2 hydrogenation of metals. Comparably, non-linear regression involves solving the reaction rate differential equation. If coupled to an energy balance, this would better represent the key changing parameters within the system, significantly s s ; thus, mitigating the issues of the integral method.

Site availability
Analysis of the site availability results in several key observations. If the equilibrium pressure cannot be greater than the gas pressure, then site de-activation would limit the initial temperature spike. This indicates that with a system of exclusively hydrogen and metal, thermal runaway cannot occur. As the rise in temperature increases site de-activation, an improvement in heat transfer would improve site-reactivation and thus increase the reaction rate. Further, as temperature varies with space and time, the wall region would experience relatively no site deactivation, whereas in the reactor centre, the opposite would occur. This important observation implies that regions experiencing a higher rate of heat transfer will transition to slower phases sooner, such as the transition from aþb to b phase. Thus, modelling localised phase transitions sufficiently is required.
However, once the reaction reaches isothermal conditions (latter reaction stage), site de-activation is no longer influential. It is proposed that the site availability is still necessary, but s s is now representing phase behaviour effects. As the metal reacts, an increase in H-H interactions occur. If these interactions cause sites to be blocked, rendering these sites unavailable e then the reaction rate reduces & capacity drops. Therefore, the site availability could be a generalised measure of the phase behaviour (and temporary site de-activation) hence; to give a good approximation of the site availability throughout the reaction, it is speculated that an accurate representation of the PCT or "phase transition curve", which is related to H-H interactions, at operating conditions is required.
We note that although s s can recognise the change in phase behaviour across each time step, it can only describe the reaction endpoint (full capacity) based on the equation used to model the equilibrium pressure. Within the equations in this paper, the capacity is only determined through a simple empirical relation (Eq. (20)) and thus this is a pitfall. Given that Eq. (20) estimates the maximum capacity when the reaction has completed, this can be influenced by the presence of a surface oxide. Although most of the magnesium activated, there is a possibility that a certain fraction did not, blocking potential sites and reducing capacity. However, we are confident that the oxygen content does not significantly influence the reaction rate as the empirical relation (4.5) mainly governs the reaction endpoint, and has little effect on the kinetic curve trajectory, which is primarily governed by the rate-determining step.

Conclusion
Non-linear regression of a derived rate law considering both diffusion and surface effects resulted in a similar effective activation energy determined through DSC desorption experiments under Argon. At conditions close to equilibrium, the analysis suggests that the reaction is surface resistance dominating with a reaction order of 5/3. There is confirmation within the rate equation that the reactive surface area increases the reaction rate. In addition, it is known that increasing the heat transfer rate improves the kinetics, where this term is explained with the concept of site deactivation. When the reaction enters isothermal conditions, it is proposed that the site availability represents the potential influence of solid phase behaviour, which can affect the reaction rate and helps to explain the reason for variation in capacity.