A critical analysis of transport models for refueling of MOF-5 based hydrogen adsorption system

Journal Pre-proof A critical analysis of transport models for refueling of MOF-5 based hydrogen adsorption system P. Sridhar, Niket S. Kaisare

PII:

S1226-086X(20)30058-7

DOI:

https://doi.org/10.1016/j.jiec.2020.01.038

Reference:

JIEC 4958

To appear in:

Journal of Industrial and Engineering Chemistry

Received Date:

26 July 2019

Revised Date:

22 December 2019

Accepted Date:

30 January 2020

Please cite this article as: Sridhar P, Kaisare NS, A critical analysis of transport models for refueling of MOF-5 based hydrogen adsorption system, Journal of Industrial and Engineering Chemistry (2020), doi: https://doi.org/10.1016/j.jiec.2020.01.038

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

A critical analysis of transport models for refueling of MOF-5 based hydrogen adsorption system

of

P. Sridhar, Niket S. Kaisare*

Department of Chemical Engineering, Indian Institute of Technology Madras, Chennai

Author: [email protected] Phone: [+91] (44) 22574176

re

-p

*Corresponding

ro

600036, India

ur na

lP

Graphical Abstract

An analysis of models for refueling of hydrogen in MOF-5 packed-bed revealed tht the choice of

Jo

adsorption isotherm has significant effect on system dynamics for hydrogen uptake as well as temperature contours within the bed.

Highlights 

The role of various transport equations for H2 adsorption in MOF-5 bed was investigated



The choice of adsorption isotherm model, though fitted to the same data, had significant impact on prediction of adsorbent bed dynamics



Temperature-dependent thermo-physical properties have moderate effect on predictions



Other parameters and assumptions had comparatively minor effect on model predictions

of

Abstract This work presents a critical analysis of 2D model of refueling a hydrogen storage tank containing

ro

MOF-5 adsorbent under cryogenic and room temperature conditions. Three isotherm models, viz. Unilan, modified Dubinin-Astakhov DA and Tóth, were fitted with experimental data of hydrogen

-p

adsorption from the literature, and their effect on refueling of adsorbent bed is analyzed. The choice of isotherm has the greatest impact on predicting refueling behavior, for both flow-through and end-

re

flow (i.e., batch) systems. A comparison between the ideal gas assumption and virial equation of state indicates that the former reasonably captures the overall behavior at lower computational cost.

lP

Analysis of the role of momentum balance, solved as full Navier-Stokes model and with Darcy’s approximation, indicates same hydrogen uptake in the adsorbent bed and only a minor difference in

ur na

the predicted velocity profiles. Finally, effect of temperature-dependent physical properties of gas, adsorbent and walls is analyzed.

Keywords

Jo

Hydrogen storage; MOF-5; Adsorption; Finite element method; Momentum transport

INTRODUCTION Energy and transport systems predominantly rely on non-renewable fossil fuel sources. There is a growing realization that emission of pollutants and greenhouse gases needs to be reduced.

Hydrogen is considered a potential energy carrier owing to higher energy efficiency and lower emissions of hydrogen-powered vehicles compared to the conventional onees that use fossil fuels [1]. An efficient hydrogen storage method, which has minimum impact on design and performance of the vehicle, is required [2]. Onboard storage of hydrogen by adsorption at cryogenic conditions (liquid nitrogen temperature) and moderately high pressures (less than 60 bar) is considered competitive compared to other storage technologies, such as liquid hydrogen, compressed gas, or

of

hydrides [3][4]. Adsorptive storage is simpler and cleaner than onboard conversion of

ro

hydrocarbon, and is safer than compressed or liquefied hydrogen [3]. Adsorptive storage can offer good gravimetric capacities, and rapid and reversible kinetics at moderate pressures [5]. In recent

-p

years, metal-organic frameworks (MOFs) have received attention for adsorptive hydrogen storage

re

because of their well-defined structures, high pore volumes [6], and high hydrogen uptake capacities, especially at cryogenic conditions [7].

lP

Adsorption on an adsorbent is characterized over a wide range of temperature and pressure operating conditions using an adsorption isotherm model [3]. These models differ in their ability

ur na

to fit the amount of gas adsorbed at equilibrium, interpretation of the observed phenomena, and prediction of the adsorption behavior beyond data [8]. The isosteric heat of adsorption, which is an important thermodynamic property, can be calculated from the isotherm models using ClasiusClapeyron equation [9]. Various adsorption isotherm models have been used for adsorption of

Jo

hydrogen at cryogenic conditions and high pressures. Dundar et al. [10] compared the Unilan, modified Dubinin-Astakhov (D-A) and MPTA-MDP for hydrogen adsorption in MOF-5 powders. Purewal et al. [11] studied the effect of densification and addition on graphite on hydrogen uptake in MOF-5 pellets of varying density. They used Unilan, modified D-A and Tóth isotherm models in the range of 77 K to 300 K for compressed pellets of MOF-5. Although the adsorption isotherm

models showed reasonable fits with the experimental data, both authors concluded that the Unilan isotherm model provided the best prediction of the data in the entire temperature range [10][11]. However, the effect of the choice of adsorption isotherm models in predicting the system-level dynamics of a packed bed containing MOF-5 adsorbent has not yet been investigated. In this paper, we apply two mono-layer adsorption models (Unilan and Tóth) and one porefilling adsorption model (modified D-A) to study the system performance during charging of

of

hydrogen into a fixed-bed with MOF-5 adsorbent. The effect of using temperature-dependent heat

ro

of adsorption is also investigated. Specifically, analytical expressions for heat of adsorption can be derived for the corresponding models from Clausius-Clapeyron equation [12]. This heat of

-p

adsorption shares common parameters with the isotherm model, so that the model parameters

re

computed by fitting equilibrium adsorption data are sufficient to compute heat of adsorption as well. Kloutse et al. [12] compared the predictions of heat of adsorption computed from Unilan,

lP

modified D-A and Tóth models, but its role on fixed-bed dynamics was not investigated. Cryo-adsorption at low temperatures and high pressures may be needed for meeting the

ur na

storage capacity targets. At these conditions, the ideal gas assumption may not be valid. Hence, even though modeling real gas behavior is computationally involved, an appropriate equation of state is used. When the operating conditions span a range of temperatures and pressures, the physical properties may vary, and temperature-dependent expressions are used. For example,

Jo

Ortmann and Kaisare [13] used Unilan isotherm, virial equation of state and modeled thermophysical properties of hydrogen and solids as polynomial functions of temperature and pressure; Hardy et al. [14] applied D-A isotherm, real gas equation with detailed treatment of heat of adsorption; whereas Kumar [15] used Langmuir adsorption isotherm with varying thermo-physical properties, but assumed a constant heat of adsorption of –4 kJ/mol. Furthermore, various

assumptions have been used for momentum balance equations. Hardy et al. [14] used a detailed treatment of momentum balance for porous bed accounting for inertial, viscous and body forces. Xiao et al. [16] applied computational fluid dynamics (CFD) approach to study bed dynamics, with Navier-Stokes equation used for momentum balance. In contrast, Momen et al., [17] used Darcy’s law (with experimentally computed constant permeability) to simplify the pressure–flow coupling during charging of hydrogen in cryo-adsorption. Simplified treatment of momentum

of

balance is often used in works that employ finite-element-based solution approach, where the

ro

model equations are solved at each mesh element. Poly-tree based adaptive mesh refinement techniques [18][19] have been recently studied for accurate solutions at lower computational cost.

-p

Such approaches have been used in solid mechanics for cracked structures [19], multi-material

re

topology [20], and fluid-submerged break water interactions [21]. In fixed bed adsorber, complex geometry of bed internals may require mesh refinement to yield accurate results. We apply free

lP

trigonal elements with extremely finer mesh for the domain, without mesh refinement. In this work, a comparative analysis of various model choices, including adsorption isotherm

ur na

model, heat of adsorption, ideal gas assumption, momentum balance assumptions, and temperature invariant physical properties is presented. The paper is organized as follows: the next section describes the adsorbent bed geometries and the model formulations, focusing on both end-flow (semi-batch) and flow-through geometries. Various modeling assumptions in the literature are also

Jo

summarized in this section. Thereafter, the three adsorption isotherms models are introduced and their role in the refueling is discussed. In the subsequent section, various assumptions used in the literature are evaluated specifically. The paper ends with the conclusion of the work.

MODEL FORMULATION Adsorbent Bed Geometry Figure 1 shows a stainless-steel vessel containing a fixed-bed of MOF-5 adsorbent modeled in this work. This hydrogen storage tank is adopted from the one considered by Chahine et al., [22] and Ubaid et al. [23], for hydrogen adsorption on activated carbon. Figure 1(a) and 1(b) show the

of

end-flow and flow-through systems, respectively. An end-flow system is similar to a (semi-)batch system, where hydrogen gas is introduced in the fixed-bed and gets either adsorbed or stored in

ro

inter-particle pores; whereas a flow-through system is a continuous system with a hydrogen inlet

-p

as well as the outlet. The geometric parameters are indicated in Figure 1, whereas the material

re

properties of MOF-5, stainless steel, and hydrogen are given in Table 1.

Governing equations

lP

The refueling behavior of the fixed-bed adsorber containing MOF-5 adsorbent is simulated using finite element package, COMSOL Multiphysics 5.2a. A two-dimensional axisymmetric

ur na

geometry with “extremely finer mesh” is solved. The model consists of partial differential equations written for mass, momentum and energy conservation. Local mass and thermal equilibrium are assumed, and intra-particle mass transfer resistance is ignored. To elucidate the role of individual terms in predicting the system refueling dynamics, a generic model [24] will be

Jo

discussed first, followed by a discussion on various simplifying assumptions made in the literature [25]. The governing equations are same for end-flow and flow-through systems, but they differ in the boundary conditions.

Mass conservation equation The mass balance is derived noting that hydrogen is present as an adsorbed phase as well as gaseous phase in the voids. The mass balance equation is written as: 𝜕(𝜌𝑔 𝜀𝑏 +𝑀H2 𝜌𝑏 𝑛𝑎 ) + ∇. [𝜌𝑔 𝑣] = 0 𝜕𝑡

(1)

where, 𝑛𝑎 is the absolute amount of hydrogen adsorbed in moles per unit mass of adsorbent and

𝜕𝑛𝑎 𝜕𝑡

(2)

is the source term in a standard mass balance equation. Various mass balance

-p

where, 𝑆𝑚 = 𝜌𝑏

𝜕𝜌𝑔 + ∇. [𝜌𝑔 𝑣] = −𝑀𝐻2 𝑆𝑚 𝜕𝑡

ro

𝜀𝑏

of

𝜌𝑔 is the density of gaseous hydrogen. The above equation is rearranged into the following form

equations used in the literature can be converted to the above form. Axial dispersion has also been

re

considered by some authors, such as Delahaye et al. [25]. For most cases of interest, Péclet number

lP

is high enough that this term may be neglected. The overall density is calculated from the pellet density and bed void fraction as:

(3)

ur na

𝜌𝑏 = (1 − 𝜀𝑏 )𝜌𝑝

Note that the absolute amount of adsorbed hydrogen may be expressed as kg hydrogen per unit mass of adsorbent, i.e., 𝑞 = 𝑀H2 𝑛𝑎 . Various formulations used in the literature are summarized in Table 2. With the exception of Momen et al. [17], all the other models can be converted into similar

Jo

form. Specifically, due to very low flowrates compared to other works, Momen et al. [17] include axial and radial diffusion as well. Momentum balance The velocity of hydrogen in the porous bed is obtained using some form of momentum balance equation. Various types of model have been used in the literature, with Darcy’s law being

the most common choice. CFD-based approaches tend to favor the use of a modified Navier-Stokes equation. Modified Navier-Stokes (N-S) equation have viscous and inertial terms added to Darcy’s law. The Brinkmann model is an intermediate between Darcy’s law and Navier-Stokes equation since it includes viscous term added to Darcy’s equation or inertial term dropped from the N-S equation [26]. The nonlinear unsteady N-S momentum balance equation, with correction for flow

𝜕𝑣 + 𝜌𝑔 (𝑣. ∇𝑣) = −∇𝑃 + ∇. [𝜇(∇𝑣 + ∇𝑣 𝑇 )] 𝜕𝑡 −∇. [(

𝜇 1 F = − 𝑣 − 𝐶2 𝜌𝑔 𝑣 2 𝜅 2

-p

where,

(4)

2𝜇 1 − ɳ𝑑 ) ( ) (∇. 𝑣)𝐼] + F 3 𝜀𝑏

ro

𝜌𝑔

of

through fixed bed is defined as:

(5)

3.5 (1 − 𝜀𝑏 ) , 𝐷𝑝 𝜀𝑏3

1 150(1 − 𝜀𝑏 )2 = 𝜅 𝐷2𝑝 𝜀3

lP

𝐶2 =

re

is the momentum source term, C2 is inertial resistance coefficient and 𝜅 is the permeabiity:

𝑏

Darcy’s law results in significant simplification, and is often used in the literature that

ur na

employs finite element based solvers [22] and in-house codes [27]. Darcy’s law describes the relationship between superficial bulk gas velocity and pressure drop[15] as: 𝜅 𝑣 = 𝜌𝑔 (− 𝛻𝑃) µ

(6)

Jo

The above equation is substituted in Eq. (2) to yield the following “Darcy’s law approximation” for the porous bed: 𝜕(𝜀𝑏 𝜌𝑔 ) 𝜅 𝜕𝑛𝑎 + ∇. [𝜌𝑔 (− ∇𝑃)] = −𝑀𝐻2 𝜌𝑏 𝜕𝑡 µ 𝜕𝑡

(7)

In summary, either a full N-S equation (5) is solved accounting for all major resistances in the system, or a modified species balance based on Darcy’s law approximation equation (7) is used in the literature. Both these approaches will be compared in this work. Energy conservation in Adsorbent Bed Myers [28] describes the use of solution thermodynamics to derive the energy balance. The

of

total enthalpy and internal energy in the system comprise of those for the bulk gas phase hydrogen, adsorbed hydrogen and the adsorbent (MOF-5). The derivation of the overall energy balance,

𝜕𝑇 + 𝜌𝑔 𝑐𝑝𝑔 𝑣. ∇𝑇 𝜕𝑡

(8)

-p

(𝜀𝑏 𝜌𝑔 𝑐𝑝𝑔 + 𝜌𝑏 (𝑞𝑐𝑝𝑔 + 𝑐𝑝𝑠 ))

ro

assuming local equilibrium, is given in Kumar and Kumar [24] as:

𝜕𝑝 − (1 − 𝛾𝑇)𝑣. 𝛻𝑝] 𝜕𝑡

re

= ∇. (𝑘eff ∇𝑇) + 𝑆𝑚 (−∆𝐻) + [𝛾𝑇𝜀𝑏

where 𝑐pg and 𝑐𝑝𝑠 are the specific heat capacities of hydrogen and adsorbent, respectively,

lP

𝑆𝑚 (−∆𝐻) represents heat of adsorption whereas the last term in Eq. (8) is the compression heat source. The use of ideal gas equation also simplifies the compression heat source. The effective

ur na

thermal conductivity, 𝑘eff is given by:

k eff = 𝑘𝑏 (1 − 𝜀𝑏 ) + 𝑘𝑔 𝜀𝑏

(9)

Delahaye et al. [25] have their model general enough to include asymmetric thermal

Jo

conductivity. Various assumptions used in the literature are summarized in Table 2. Energy Conservation in Solid Walls Finally, the energy balance equation is also written for the solid stainless steel wall as: (𝜌𝑤 𝐶𝑝𝑤 )

𝜕𝑇𝑤 = 𝛻. (𝑘𝑤 ∇𝑇) 𝜕𝑡

(10)

All the physical properties are summarized in Table 1, whereas temperature-dependent thermophysical properties of the steel wall are discussed in Appendix C. Table 2 summarizes the most common modeling options for mass, momentum and energy balance equations for adsorbent storage bed in the literature. Authors have used in-house finite difference codes [17][30], finite element [14][22] and finite volume (CFD) [29] solvers for modeling hydrogen adsorption in fixed-bed. The major difference is the use of different adsorption

of

isotherm, such as Langmuir [30], Radke-Prausnitz [25] and D-A [31][17][14][22]. Although

ro

majority of the studies applied D-A isotherm, Dundar et al. [10] reported Unilan isotherm provides best fit for hydrogen adsorption on MOF-5. Another key observation is that the Darcy’s law and

-p

real gas equation is often preferred by authors using finite-element-based commercial solvers,

lP

Operating and Boundary Conditions

re

whereas Navier-Stokes models and ideal gas assumption is more common in CFD.

Based on the survey of literature, we have chosen refueling of adsorbent bed under cryogenic and room-temperature conditions. In this work, both end-flow (semi-batch) and flow-through

ur na

(continuous) systems are modeled.

The initial condition in cryogenic refueling represent the fixed-bed adsorber under depleted condition: 1.1 bar pressure, 140 K temperature and the adsorbed hydrogen in equilibrium as

Jo

determined by the isotherm model. At the inlet, hydrogen is charged into the system at 80 K and 8.04926 kg/m2-s mass flux (which yields 0.4 g/s flowrate). The walls of storage tank are insulated. For the end-flow system, there is no outlet and hydrogen is charged for 40 seconds. The flow-through system is operated assuming downstream pressure of 20 bar. Flow-through refueling system operates in two steps. First, hydrogen is fed into the fixed-bed by keeping it in

batch mode until the pressure near the outlet reaches 20 bar. Thereafter, the outlet of the bed is “opened”, which is modeled using pressure-outlet boundary condition with 20 bar outlet pressure. The operating conditions for room-temperature refueling are adopted from Chahine et al. [22], where the tank is initally depleted at 0.3208 bar and 281 K. Hydrogen is introduced in the system at 297 K temperature and mass flux of 0.402463 kg/m2-s (i.e., 0.02 g/s net flowrate). The

of

system is assumed to be placed in an ice bath at 282 K.

ro

ROLE OF ADSORPTION ISOTHERM Isotherm Models

-p

Adsorption isotherms, which quantify the amount of hydrogen adsorbed across a wide range

re

of temperature and pressure, are required for modeling the adsorption process. The Unilan and modified Dubinin-Astakhov (DA) isotherm models have been used for hydrogen adsorption on

lP

MOF-5 powders over a broad range of pressures and temperatures[10]. Purewal et al. [11] used Unilan, modified D-A and Tóth models for densified pellets of MOF-5. Isosteric heat of adsorption

ur na

can be calculated using Clausius-Clapeyron relation from these models. This section discusses the effect of the choice of isotherm models on predicting refueling in an adsorbent bed under cryogenic and room-temperature conditions. For sake of completeness, adsorption isotherm models were refitted to the experimental data of Zhou et al.[32] and the fits are discussed in Appendix A.

Jo

Unilan Isotherm

Unilan, which stands for uniform energy distribution and Langmuir local isotherm, is an

empirical model that uses monolayer Langmuir equation to describe the local isotherms, but averages the Langmuir equation over a continuous energy interval Emin to Emax. Dundar et al.[10] report that Unilan model gives a better fit to adsorption data for hydrogen on MOF-5 in 77K and

298K temperature range. Likewise, Purewal et al.[11] concluded that the Unilan isotherm model is more appropriate for MOF-5 due to its larger pore-size and lower surface heterogeneity. Absolute adsorption expression for Unilan model given by[11][33], 𝐸

(11)

max 𝑎 + 𝑝 exp ( 𝑅𝑇 ) 𝑛𝑚𝑎𝑥 𝑅𝑇 𝑛𝑎 = 𝑙𝑛 [ ] 𝐸 𝐸𝑚𝑎𝑥 − 𝐸𝑚𝑖𝑛 𝑎 + 𝑝 exp ( min )

𝑅𝑇

of

The four parameters required to fit the experimental adsorption data include: The maximum amount adsorbed, 𝑛max (mol/kg), the minimum and maximum energies, Emin and Emax (J/mol), and

ro

the entropy difference ∆S (J/mol.K). The heat of adsorption for the Unilan model, using ClausiusClapeyron equation is given by[12][13]

(1−𝑥)𝑠

1 − exp (

where, 𝑥 = 𝑛⁄𝑛max and 𝑠 = 𝐸𝑚𝑎𝑥 − 𝐸𝑚𝑖𝑛

)

−

𝑥𝑠

1 − exp (𝑅𝑇)

(12)

lP

Modified Dubinin-Astakhov Isotherm

𝑅𝑇

𝑥𝑠

-p

(1 − 𝑥)𝑠

re

Δ𝐻 = 𝐸𝑚𝑎𝑥 −

The modified Dubinin-Astakhov (D-A) isotherm is based on Polanyi adsorption potential

ur na

theory[34]. The adsorbate attraction with the surface is assumed to be due to the Vander Waals interactions[34]. This model is widely applied for micropore materials like activated carbon[35].

Jo

The modified Dubinin-Astakhov isotherm model is given by[36][37], 𝑚

𝑅𝑇 𝑝𝑜 𝑛𝑎 = 𝑛𝑚𝑎𝑥 exp [− ( 𝑙𝑛 ( )) ] 𝛼 + 𝛽𝑇 𝑝

(13)

where, 𝛼 is enthalpic factor and 𝛽 is entropic factor, 𝑝𝑜 is the saturation pressure of the bulk gas, and 𝑝 is the operating pressure. The analytic expression for isosteric heat of adsorption (∆𝐻) for modified DA isotherm is given by[38]:

(14)

𝑛𝑎 ∆𝐻 = 𝛼 √− ln ( ) 𝑛𝑚𝑎𝑥 𝑚

The five parameters need to be fit with the experimental adsorption data include 𝑛max , 𝛼, 𝛽, 𝑝𝑜 and 𝑚 . Along with Unilan, modified D-A is the most commonly-used isotherm to describe adsorption of hydrogen in microporous adsorbents over cryogenic temperature range.

of

Tóth Isotherm model Tóth is an empirical model developed from potential theory for heterogeneous adsorption, to

ro

model low and high-pressure adsorption. It assumes Quasi-Gaussian energy distribution with the

𝑏0 exp(𝑞𝑇 /𝑇) 𝑝 (1 + (𝑏0 exp(𝑞𝑇 /𝑇) 𝑝)𝑛 )1/𝑛

(15)

re

𝑛𝑎 = 𝑛max

-p

adsorption sites[39]. The absolute amount of hydrogen adsorbed for the Tóth model is given by:

where 𝑝 is the equilibruim pressure (Pa) and the heterogeneous parameter, 𝑛 = 𝑛𝑜 + 𝛼𝑇, depends

lP

on temperature. The parameters 𝑏0 , 𝑞𝑇 , 𝛼, 𝑛0 are fitted to the experimental data of hydrogen adsorption on MOF-5[32]. The heat of adsorption is given by[12], 𝛼 𝑛𝑎 𝑅𝑇 2 [ln ( ) (1 + 𝑏𝑝𝑛 ) − ln(𝑏𝑝)] 𝑛 𝑛𝑚𝑎𝑥

(16)

ur na

𝛥𝐻 = 𝑞𝑇 +

In Appendix A, we have refitted the three isotherms to the same data of Zhou et al. [32].

Jo

Table 3 shows the parameters obtained following the procedure in Appendix A. The model parameters for Unilan and modified D-A models are similar to Dundar et al.[10]; we also fitted the parameters for the Tóth isotherm model using the same procedure. The root means square error (RMSE) is also reported in Table 3 and model fits are shown in Figure A-1. These indicate a slightly better fit for the Unilan model in spite of having one less fitting parameter, although all three models show reasonable fits to the isotherm data.

Refueling of End-Flow System Refueling under Cryogenic Conditions We first compare the effect of adsorption isotherm model on the refueling behavior of an end-flow system under cryogenic conditions. All nonidealities are considered in the baseline model: Material properties are given in Table 1, hydrogen is modeled as real gas using virial EOS

of

(see Appendix B) and all the physical properties are assumed to be temperature-dependent (see

ro

Appendix C). The system, initially at 140 K temperature and 1.1 bar pressure, is refueled using hydrogen at 80 K in an insulated stainless-steel tank. Figure 2 compares the temperature contours

-p

from the three adsorption isotherms models. In all cases, the adsorbent bed gets cooled due to the incoming hydrogen at 80 K temperature. However, the heat released on adsorption and rapid

re

pressurization due to refueling raises the temperature, especially in the middle of the tank.

lP

The comparison between the amount of hydrogen adsorbed, predicted by the three models is shown in Figure 3. This figure shows the average amount of hydrogen adsorbed, which is obtained

ur na

by computing the total weight of adsorbed H2 divided by the total weight of adsorbate in the tank. At lower pressures (until about 15 seconds of refueling) the temperature contours, as well as an amount of hydrogen adsorbed, show reasonable match for all the three isotherms. However, towards the end of 40 seconds, a noticeable difference is observed in temperature contours.

Jo

Specifically, the low-temperature region near the entrance is smaller for the Unilan model. Since local mass equilibrium is established, the amount of hydrogen adsorbed is correspondingly lower when using Unilan isotherm. Recall that the isotherm models were fitted to the same data, and the overall fits seemed reasonable, both from the RMSE values (Table 3) and eyeballing the model fits (Figure A-1). The different rates of adsorption affect the overall temperature profile.

Figure 4 shows the heat of adsorption profiles as a function of time to further analyze the observed behavior. The heat of adsorption decreases with time for all the isotherm models because the saturation of high energetic sites leads the molecules to get adsorbed on low energetic sites. Modified DA model shows the highest ΔH initially, which drops rapidly as sites get filled. While the heat of adsorption is higher for modified DA model for the first ten seconds, Unilan model shows highest heat of adsorption thereafter, and is least for Toth model. Some differences in the

of

temperature profiles are noticeable in temperature contours of Figure 2, but these differences are

ro

not significant because the heat of adsorption is similar for all the three models.

-p

Refueling under Room-Temperature Conditions

Recall that the isotherm models were fitted for 77 K to room temperature data and can be

re

used for room-temperature refueling as well. Since the model performance in capturing the adsorption isotherm varies over the range of temperatures, we now compare refueling of end-flow

lP

system under room-temperature refueling conditions. For sake of brevity, we only show the absolute amount of hydrogen adsorbed in the adsorbent bed. The results of Figure 5 indicate that

ur na

the adsorption isotherm influences refueling behavior in end-flow for room-temperature case as well. It should be noted that the inlet hydrogen flowrate is an order of magnitude lower than the cryogenic case, resulting in a higher refueling time of 300 seconds. The difference in prediction

Jo

from the Unilan model may be attributed to the lower adsorption volume (𝑣𝑎 ) and the fact that it provides a better fit to the data at higher temperatures. These results are qualitatively similar to those in cryogenic refueling, indicating that the choice of isotherm model affects the overall prediction of hydrogen uptake.

Refueling of Flow-Through System under Cryogenic Conditions Temperature contours for flow through refueling with 20 bar downstream pressure and a spent adsorbent bed (initial Temperature 140 K and pressure 1.1 bar) are shown in Figure 6. Initially, as hydrogen enters the system at 80 K, the bed starts cooling down and the pressure increases. Until the pressure reaches 20 bar, there is no hydrogen outflow and the bed resembles

of

the end flow system. The pressure reaches 20 bar for Unilan in 17 seconds, modified DA and Toth models take 22 seconds and 20 seconds, respectively. Thereafter, the pressure outlet condition is

-p

from the system and pushes the hot region out of the system.

ro

used and unadsorbed hydrogen exits the bed. The incoming cold hydrogen reduces the temperature

Figure 7 shows that the absolute amount of hydrogen adsorbed increases with time until the

re

system saturates at around 200 seconds for all three models. The initial shaded region in the figure

lP

represents the time taken for the bed to reach pressure of 20 bar (i.e., 17 seconds with Unilan model), until which time the bed operates similar to an end-flow case. When the predictions from the isotherm models are compared, we observe that the Toth and modified DA models over-predict

ur na

the amount adsorbed compared to the Unilan model. The three models show differences in the absolute amount of hydrogen adsorbed, which could be attributed to the different adsorption volume, 𝑣𝑎 , as well as to different iso-steric heat of adsorption (ΔH) and error in isotherm fitting.

Jo

Finally, Figure 8 shows the heat of adsorption comparison profile of the three models for flow through refueling case. The heat of adsorption decreases with the time and Unilan model predicted 3.5 kJ/mol, whereas modified DA and Toth models estimated 2.5 kJ/mol by the end of refueling. In summary, although the three isotherm models predict the adsorption isotherm data reasonably, choice of isotherm has a significant effect on predicting refueling dynamics, especially across a wide range of temperature.

RESULTS AND DISCUSSIONS We now investigate the effect of individual terms in flow-through refueling of hydrogen under cryogenic conditions in the remainder of this paper. The simulations for the end-flow case, under room-temperature conditions, gave qualitatively similar results and are skipped for brevity. Effect of Heat of Adsorption

of

The isosteric heat of adsorption is a measure of the interaction between adsorbate molecules

ro

and the surface atoms of the adsorbent. Figure 8 shows that the isosteric heat of adsorption varies with time. Due to energetic heterogeneity of the solid surface, the most energetic sites get filled

-p

first, resulting in a high value of (– 𝛥𝐻) initially, which drops as adsorption proceeds[40]. To

re

further analyze the effect of heat of adsorption on refueling behavior, we compared the case of constant heat of adsorption (𝛥𝐻 = −4 kJ/mol) and the expression given in Equation (12). Figure

lP

9 shows the adsorbate profile in the bed for flow-through refueling for the two cases. Operating conditions used is similar to the previous cryogenic flow through cases. The adsorbate amount

ur na

predicted is the same for both “regimes”, the initial 17 seconds when the bed pressurizes and the remainder (17 to 300 s) when the hydrogen breaks through. The difference in the heats of adsorption is small (𝛥𝐻 varies from –3.9 to –3.45 kJ/mol during the refueling) and therefore has a negligible effect on predictions of the absolute amount of hydrogen adsorbed in the bed. We

Jo

verified that using constant 𝛥𝐻 gives the exact same pressure profiles and negligible variation of upto 2 K in the temperature profiles (results not shown). Hence, with Unilan isotherm model, it is feasible to use a constant heat of adsorption since it does not significantly affect the overall refueling dynamics.

Momentum balance As seen in Table 2, various assumptions have been used in the literature to model the flow in the packed bed of adsorbent. In this section, two of the popular approaches are examined. Solving the full Navier Stokes model of Equation (4)-(5) is computationally more demanding, since inertial, body-force and viscous terms are included in the momentum balance. In contrast, with the

of

simplified Darcy’s formulation (Equation 6), the momentum balance equation is not solved independently, but is used to yield the following mass balance for the porous bed:

ro

𝜕(𝜀𝑏 𝜌𝑔 ) 𝜅 𝜕𝑛𝑎 + ∇. [𝜌𝑔 (− ∇𝑃)] = −𝑀𝐻2 𝜌𝑏 𝜕𝑡 µ 𝜕𝑡

Eq. (7)

-p

The use of the above equation is quite popular in modeling system dynamics of hydrogen storage,

re

possibly due to the popularity of COMSOL Multiphysics in this field, where this model is available as one of the default options for porous bed. Figure 10 shows that for flow-through refueling under

lP

cryogenic operating conditions, both the methods show a good agreement. The same observations were applicable for higher flowrates, as well as when we compared the refueling dynamics of end-

ur na

flow system, where no difference in predicted hydrogen profiles was observed (results skipped for brevity). Figure 11 shows the velocity contours calculated using the original Navier-Stokes equation (right panels) and Darcy’s law (left panels) at various times. The velocity contours for the inlet and outlet sections, which are of the order of 1 m/s, are not shown and only the contours

Jo

for the fixed bed are shown. This indicates that the inertial and viscous effects are negligible under conditions of interest, resulting in similar predictions from the two models. Effect of Ideal Gas Assumption Critical temperature and pressure of hydrogen are about 33 K and 13 bar, respectively. Thus, hydrogen is a supercritical fluid at the operating conditions for moderate- and high-pressure

adsorptive storage. Hence, it is typically modeled as a real gas, where the density is calculated from compressibility factor computed using the virial equation of state (Appendix B). This also affects the compression heat source term: [𝛾𝑇𝜀𝑏

𝜕𝑝 − (1 − 𝛾𝑇)𝑣. ∇𝑝] 𝜕𝑡

(18)

Heat generated by compression of the hydrogen gas has a major contribution in the total heat

−1 𝜕𝜌𝑔 ( ) 𝜌𝑔 𝜕𝑇 𝑝

ro

𝛾=

of

source. The isobaric thermal expansion coefficient is

-p

For an ideal gas, 𝛾𝑇 = 1 and the above term in Eq. (18) simplifies significantly. We therefore investigate the effect of real gas law on prediction of system dynamics. Figure 12 shows that the

re

using ideal gas approximation results in nearly the same prediction of the amount of hydrogen adsorbed as the real gas equation of state. The difference between pressures predicted at any time

lP

from real gas and ideal gas assumption is less than 0.2 bar and that in temperature is less than 1 K for the entire duration of refueling. While this work considered moderate refueling pressures,

ur na

Lamari et al. [41] argued that ideal gas is a reasonable assumption even for high pressure refueling, since the compressibility factor varies less than 8% for pressures in the range of 1 to 150 bar. In summary, ideal gas approximation, which is computationally beneficial, gives similar

Jo

results as the real gas model in predicting the system dynamics of hydrogen adsorption. Effect of Temperature-Dependent Properties The polynomial fits for themo-physical properties of hydrogen, MOF-5 and steel wall are

given in Appendix C, whereas constant values at representative conditions are summarized in Table 1. Thermal conductivity and specific heat capacity are sometimes assumed constant during simulation studies. During flow-through refueling, the system operates over a wide range of

temperature and pressure, which affects the values of these properties. Figure 13 shows the effect of assuming constant thermo-physical properties. Solid line represents the base-case for Unilan model, with all properties computed using temperature-dependent relations. These results are compared with dotted and dashed lines, which represent the cases when only solid-phase properties and hydrogen properties were held constant, respectively. All three cases show relatively similar eventual behavior (i.e., towards the end of refueling). However, when the properties were assumed

of

constant, the model over-predicted the refueling rate, especially in the middle of the process. For

ro

example, with variable thermophysical properties, the net amount of hydrogen stored in the bed was predicted as 3.35% (by weight) at 100 seconds. At the same time, using constant adsorbent

-p

properties predicted 3.67% by weight and constant gas properties predicted 3.68% by weight. We

re

repeated the same study for room-temperature refueling in an end-flow system and observed that

lP

the effect of thermo-physical properties was relatively insignificant.

In summary, unlike equation of state and heat of adsorption, using temperature-dependent

ur na

thermophysical properties show low-to-moderate effect on refueling dynamics. CONCLUSIONS

This work analyzed the role of various transport model equations for hydrogen refueling in

Jo

MOF-5 system and presented a comparison between several options considered in the literature, such as adsorption isotherm models, ideal gas assumption, momentum balance and the temperature dependency of the thermophysical properties. The major findings of this study are as follows: 

The choice of adsorption isotherm model has a significant impact on the prediction of hydrogen uptake in adsorbent bed during refueling. Although Unilan, modified D-A, and Tóth models showed reasonable fit with the excess adsorption isotherm, they differed in

predictions of the amount of hydrogen adsorbed for cryogenic as well as room-temperature refueling. The effect was maximum for flow-through refueling under cryogenic conditions. 

Further analysis revealed that difference in the adsorption volume (𝑉𝑎 ), error of the model fits and heat of adsorption are likely to be responsible for the observed difference between predictions from the isotherm models. With Unilan isotherm model, constant heat of adsoprtion resulted in similar hydrogen uptake

of



behavior for all refueling cases. Likewise, ideal gas assumption provided similar behavior as

ro

virial equation of state. Consequently, we recommend using ideal gas assumption and constant heat of adsorption (–4 kJ/mol) in order to simplify the computational requirements. Since the system operates under laminar conditions, there was no prominent effect of using

-p





re

Darcy’s law approximation, instead of solving a full momentum balance for porous bed. Choosing constant values of adsorbent and gas specific heat capacities slightly overpredicted

through the case.

lP

the hydrogen adsorption amount than using temperature-dependent properties for the flow-

ur na

In conclusion, the choice of the adsorption isotherm model is crucial, whereas temperaturedependent thermo-physical properties have some effect, in system-level simulation of adsorbent bed. On the other hand, the use of ideal gas assumption, constant heat of adsorption and Darcy’s approximation are recommended to simplify the calculations under typical conditions relevant to

Jo

refueling of hydrogen in an adsorbent bed.

APPENDIX A: EFFECT OF ISOTHERM ON HYDROGEN STORAGE Purewal et al. [11] compared Unilan, modified Dubinin Astakhov and Toth isotherm models for pellets of MOF-5 with and without expanded natural graphite (ENG). They concluded that for

open pore structure and homogenous nature of storage sites, Unilan is the most appropriate model. Dundar et al.[10] compared Unilan and modified DA isotherm models, which were fitted to data of Zhou et al. [32] in the range 77-300 K. In this appendix, for sake of consistency, we refit the Unilan and DA model for the same data (Zhou et al. [32]), as well as obtain parameters for Tóth isotherm model. The adsorption volume, 𝑣𝑎 , is simultaneously fitted along with all other model parameters. In case of the Unilan isotherm, a total of five parameters (including 𝑣𝑎 ) are fitted,

of

whereas six parameters (including 𝑣𝑎 ) are fitted for the modified DA and Tóth isotherms.

ro

The experimental excess adsorption data was obtained by Zhou et al. [32] at various pressures and temperatures. A nonlinear least squares approach is used to obtain the optimum values of the

-p

parameters that minimize the root mean square error (RMSE) objective:[10]

re

1 (𝑖) (𝑖) 2 𝐹=√ ∑(𝑛𝑒𝑥 − 𝑛𝑚 ) 𝑛−𝑝 (𝑖)

(18)

(𝑖)

lP

𝑛𝑒𝑥 is the excess adsorption value obtained from experiments and 𝑛𝑚 is the model prediction at the same conditions and 𝑛 is number oa f data points and 𝑝 is the number of fitted parameters. The

ur na

model parameters obtained using this procedure are summarized in Table 3, where the RMSE value is also noted. As observed by Dundar et al. [10] and Purewal et al. [11], the Unilan isotherm model shows a slightly better fit to the data, in terms of the lower RMSE values. The model predictions for the three isotherm models are plotted as lines in Figure A.1,

Jo

whereas the experimental excess adsorption data of Zhou et al. [32] is shown as symbols, for five different temperatures. These results match those reported by Dundar et al. [10], with Unilan model showing a better fit than the other two models when the entire temperature range is considered.

APPENDIX B: EQUATION OF STATE FOR HYDROGEN GAS Ortmann and Kaisare[13] fitted a third order polynomial to the compressibility data of hydrogen in the temperature range of 65 to 305 K and pressure range of 1-65 bar. The compressibility factor was given by 𝑍(𝑃, 𝑇) = 1 + 𝐵𝑧 (𝑇)𝑃 + 𝐶𝑧 (𝑇)𝑃2 + 𝐷𝑧 (𝑇)𝑃3 𝑋𝑧 = 𝑥𝑧1 +

𝑥𝑧2 𝑥𝑧3 𝑥𝑧4 + 2+ 3 𝑇𝑟 𝑇𝑟 𝑇𝑟

(A.2)

of

where,

( A.1)

ro

Here, Xz represents Bz, Cz, and Dz, and correspondingly𝑥𝑐𝑖 ≡ {𝑏𝑧𝑖 , 𝑐𝑧𝑖 , 𝑑𝑧𝑖 }. 𝑇𝑟 = 𝑇/𝑇𝑐 is the reduced temperature and 𝑇𝑐 = 33.135 K. The compressibility factor is used to compute density of

𝑃𝑀H2 𝑍𝑅𝑇

re

𝜌𝑔 =

-p

hydrogen as:

(A.3)

and the volume expansion coefficient is computed using the definition: 1 𝜕𝜌𝑔 𝜌𝑔 𝜕𝑇

lP

𝛾=−

(A.4)

ur na

Table B.1 summarizes the model parameters used in this work.

Appendix C: Thermo-Physical Properties

Jo

The effect of temperature dependent thermo-physical properties is also analyzed in this work. Like the previous section, the properties of hydrogen depend on temperature and pressure. The specific heat of hydrogen is expressed as: 𝐶𝑝 = 𝐴𝑐 (𝑝) + where,

𝐵𝑐 (𝑝) 𝐶𝑐 (𝑝) + + 𝐷𝑐 (𝑝)𝑇 𝑇 𝑇2

(A.5)

𝑋𝑐 (P) = 𝑥𝑐1 + x𝑐2 P + 𝑥𝑐3 P 2 + 𝑥𝑐4 P 3

(A.6)

Here, Xc represents Ac, Bc, Cc and Dc, and correspondingly𝑥𝑐𝑖 ≡ {𝑎𝑐𝑖 , 𝑏𝑐𝑖 ,𝑐𝑐𝑖 }. The specific heat is given in J/kg.K, temperature the in Kelvins and pressure in bar. In a similar manner, viscosity of the hydrogen (in Pa-s) is expressed as 1

𝜇𝑔 = 𝐴𝑣 (𝑃) + 𝐵𝑣 (𝑃)𝑇𝑟3 +

𝐶𝑣 (𝑃) 𝑇𝑟

of

where,

(A.7)

𝑋𝑣 (𝑃) = 𝑥𝑣1 + 𝑥𝑣2 𝑃 + 𝑥𝑣3 𝑃2

ro

(A.8)

and 𝑥𝑣𝑖 ≡ {𝑎𝑣𝑖 , 𝑏𝑣𝑖 ,𝑐𝑣𝑖 }. The parameters for specific heat and viscosity are given in Table C.1 and

-p

C.2, respectively.

re

The solid properties are also expressed as a polynomial function of temperature: 6

(A.9)

lP

𝑐𝑝𝑠 = ∑ 𝑎𝑖 𝑇 𝑖 i=0

ur na

where 𝐶𝑃 is specific heat of MOF-5 and steel wall in J/kg.K. The parameters for computing specific heat are given in Table C.3. Finally, the wall thermal conductivity is given by: log(𝑘𝑤 ) = −1.4087 + 1.3982[log 𝑇] + 0.2543[log 𝑇]2 − 0.626[log 𝑇]3

(A.10)

+ 0.2334[log 𝑇]4 + 0.4256[log 𝑇]5 − 0.4658[log 𝑇]6

Jo

+ 0.165[log 𝑇]7 − 0.0199[log 𝑇]8

Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

of

ro

-p

re

lP

ur na

Jo

REFERENCES M. Ball, M. Weeda, Int. J. Hydrogen Energy 40 (2015) 7903–19.

[2]

D.J. Durbin, Int. J. Hydrogen Energy 38 (2013) 14595–617.

[3]

R.G. Paggiaro, PhD Thesis, Technical University of Munich, Germany (2008).

[4]

L. Zhou, Renew. Sustain. Energy Rev. 9 (2005) 395–408.

[5]

V.S. Kumar, K. Raghunathan, S. Kumar, Int. J. Hydrogen Energy 34 (2009) 5466–75.

[6]

R.E. Morris, P.S. Wheatley, Gas Storage Mater. (2008) 4966–81.

[7]

M.P. Suh, H.J. Park, T.K. Prasad, D.W. Lim, Chem. Rev. 112 (2012) 782–835.

[8]

X. Tang, N. Ripepi, K. Luxbacher, E. Pitcher, Energy & Fuels 31 (2017) 10787–801.

[9]

Y. Tian, J. Wu, Langmuir 33 (2017) 996–1003.

-p

ro

of

[1]

re

[10] E. Dundar, R. Zacharia, R. Chahine, P. Bénard, Fluid Phase Equilib. 363 (2014) 74–85. [11] J. Purewal, D. Liu, A. Sudik, M. Veenstra, J. Yang, S. Maurer, U. Muller, D.J. Siegel, J.

lP

Phys. Chem. C 116 (2012) 20199–212.

[12] A.F. Kloutse, R. Zacharia, D. Cossement, R. Chahine, R. Balderas-Xicohténcatl, H. Oh, B.

ur na

Streppel, M. Schlichtenmayer, M. Hirscher, Appl. Phys. A 121 (2015) 1417–24. [13] J.P. Ortmann, N.S. Kaisare, Int. J. Hydrogen Energy 41 (2016) 342–54. [14] B. Hardy, C. Corgnale, R. Chahine, S. Garrison, D. Tamburello, D. Cossement, D. Anton, Int. J. Hydrogen Energy 37 (2012) 5691–705.

Jo

[15] V.S. Kumar, Int. J. Hydrogen Energy 36 (2011) 15239–49. [16] J. Xiao, R. Peng, D. Cossement, P. Bénard, R. Chahine, Int. J. Hydrogen Energy 38 (2013) 10871–9.

[17] G. Momen, R. Jafari, K. Hassouni, Int. J. Therm. Sci. 49 (2010) 1468–76. [18] H. Nguyen-Xuan, K.N. Chau, K.N. Chau, Comput. Methods Appl. Mech. Eng. 355 (2019)

405–37. [19] H. Nguyen-Xuan, S. Nguyen-Hoang, T. Rabczuk, K. Hackl, Comput. Methods Appl. Mech. Eng. 313 (2017) 1006–39. [20] K.N. Chau, K.N. Chau, T. Ngo, K. Hackl, H. Nguyen-Xuan, Comput. Methods Appl. Mech. Eng. 332 (2018) 712–39. [21] T. Vu-Huu, P. Phung-Van, H. Nguyen-Xuan, M. Abdel Wahab, Comput. Math. with Appl.

of

76 (2018) 1198–218.

ro

[22] R. Chahine, J. Xiao, J. Wang, D. Cossement, P. Benard, Int. J. Hydrogen Energy 37 (2011) 802–10.

-p

[23] S. Ubaid, R. Zacharia, J. Xiao, R. Chahine, P. Bénard, P. Tessier, Int. J. Hydrogen Energy

re

40 (2015) 9314–25.

[24] V.S. Kumar, S. Kumar, Int. J. Hydrogen Energy 35 (2010) 3598–609.

lP

[25] A. Delahaye, A. Aoufi, A. Gicquel, I. Pentchev, AIChE J. 48 (2002) 2061–73. [26] Z. Shi, X. Wang, Proceeding COMSOL Conf. (2007).

ur na

[27] G. Momen, G. Hermosilla, A. Michau, M. Pons, M. Firdaous, P. Marty, K. Hassouni, Int. J. Heat Mass Transf. 52 (2009) 1495–503. [28] A.L. Myers, AIChE J. 48 (2002) 145–60. [29] J. Xiao, L. Tong, D. Cossement, P. Bénard, R. Chahine, Int. J. Hydrogen Energy 37 (2012)

Jo

12893–904.

[30] J.P.B. Mota, A.E. Rodrigues, E. Saatdjian, D. Tondeur, Carbon N. Y. 35 (1997) 1259–70. [31] G. Hermosilla-Lara, G. Momen, P.H. Marty, B. Le Neindre, K. Hassouni, Int. J. Hydrogen Energy 32 (2007) 1542–53. [32] W. Zhou, H. Wu, M.R. Hartman, T. Yildirim, J. Phys. Chem. C 111 (2007) 16131–7.

[33] Y. Wang, C. Ercan, A. Khawajah, R. Othman, AIChE J. 58 (2012) 782–8. [34] M. Polanyi, Science. 141 (1963) 1010–3. [35] A.M. Czerny, P. Bénard, R. Chahine, Langmuir 21 (2005) 2871–5. [36] M.A. Richard, P. Benard, R. Chahine, Adsorption 15 (2009) 43–51. [37] M.M. Dubinin, V.A. Astakhov, Bull. Acad. Sci. USSR Div. Chem. Sci. 20 (1971) 3–7. [38] P.B. Whittaker, X. Wang, W. Zimmermann, K. Regenauer-Lieb, H.T. Chua, J. Phys. Chem.

of

C 118 (2014) 8350–8.

ro

[39] TOTH, J., Acta Chim. Hung. 69 (1971) 311–28.

[40] P.B. Whittaker, X. Wang, K. Regenauer-Lieb, H.T. Chua, Phys. Chem. Chem. Phys. 15

-p

(2013) 473–82.

Jo

ur na

lP

re

[41] M. Lamari, A. Aoufi, P. Malbrunot, AIChE J. 46 (2000) 632–46.

Figure 1: Schematic of the hydrogen storage tanks with axial symmetry for (a) end-flow and (b)

re

-p

ro

of

flow-through systems modeled in this work.

lP

Figure 2: A comparison of temperature contours in the adsorbent bed, predicted for the three different isotherm models, at various times during refueling of hydrogen in an end-flow system under cryogenic operating conditions. The insulated bed is initially at 140 K and 1.1 bar, and

ur na

hydrogen enters the bed at 80 K temperature. The three columns show results with Unilan (left), modified DA (middle) and Tóth (right) isotherm models. The full model equations are solved and

Jo

all properties are taken as temperature-dependent.

of ro

Figure 3: The effect of the three isotherm models on predicting the amount of hydrogen adsorbed

-p

in the bed for refueling of the end-flow system at cryogenic temperatures (with hydrogen inflow at

ur na

lP

re

80 K). The parameters are the same as that for Figure 2.

Figure 4: Comparison of the heat of adsorption for the various isotherm models in the end-flow

Jo

system for the conditions shown in Figure 3.

of

ro

Figure 5: A comparison of the refueling behavior of the end-flow system for the three isotherms models at room temperature (with real-gas model and temperature-dependent physical properties).

-p

Hydrogen gas is introduced into the vessel at 297 K, with the vessel kept in an ice bath (282 K). Modified D-A and Toth models significantly over-predict the amount of adsorbed hydrogen,

re

compared to the Unilan isotherm.

15 s

ur na

lP

50 s

Unilan model

Modified DA model

100 s 150 s 200 s Tóth model

Figure 6: A comparison of the temperature contours, using the Unilan (left column), modified DA (middle column) and Tóth (left column) isotherm models, at various times during flow-through

Jo

refueling under cryogenic conditions. The fixed-bed is initially at 140 K, hydrogen enters the system at 80 K and mass flux of 8.049 kg/m2.s, and the downstream pressure is kept at 20 bar.

of

ro

Figure 7: Profiles of adsorbed hydrogen in the MOF-5 adsorbent bed for flow-through refueling at cryogenic conditions (for conditions shown in Figure 6) for the three isotherm models. In the first

-p

part of refueling, there is no outflow from the system until the pressure in the bed reaches 20 bar.

ur na

lP

re

The shaded region indicates this initial phase for Unilan isotherm (17 seconds).

Jo

Figure 8: The isosteric heat of adsorption at various times for the three isotherm models, computed using Clausius-Clapeyron equation for the flow-through system under cryogenic refueling conditions.

of ro

Figure 9: The role of using ΔH computed from Clausius-Clapeyron equation (solid line) and

-p

constant ΔH (dashed line) in the prediction of the adsorbed hydrogen profile for flow-through refueling of the bed under cryogenic conditions. The solid line represents the same case as Unilan

Jo

ur na

lP

re

model in Figure 7, whereas a fixed value of (– Δ𝐻) = 4 kJ/mol is used for comparison.

Figure 10: The role of the choice of momentum balance equation on predictions of the amount of hydrogen adsorbed for flow-through refueling of the adsorbent bed under cryogenic conditions. A “simplified” Darcy model approach (dashed line) provides similar refueling results as the Navier Stokes model (solid line).

Figure 11: Contours of gas velocity within the adsorbent bed computed using Darcy’s law (left) and Navier Stokes equation (right), at various times during refueling of the flow-through system.

of

Contours are shown only for the fixed bed for clarity and not the inlet and outlet regions. Both the

lP

re

-p

ro

assumptions provide similar velocity profiles in the bed at all times during refueling.

ur na

Figure 12: Effect of the choice of the equation of state indicates that ideal gas assumption (dashed lines) reasonably predicts the dynamics of adsorbent bed refueling under cryogenic conditions. The solid line, representing the real gas equation (modeled using the virial equation of state, Appendix

Jo

B), is same as that seen in Figure 7 for Unilan model.

of ro

Figure 13: The effect of temperature-dependent thermophysical properties of the adsorbent and

-p

hydrogen on the dynamical response for flow-through refueling under cryogenic conditions (solidline: “Base-case”; dotted-line: thermo-physical properties of the wall and MOF-5 were kept

Jo

ur na

lP

re

constant; dashed-line: properties of hydrogen were kept constant).

Figure A.1: Comparison of adsorption isotherm models in range of 77-298K with experimental data of MOF-5. Circles indicate the experimental data[32], solid blue lines indicate Unilan model, dotted black lines indicate Tóth model, and dashed red lines indicate modified DA model.

Table 1: Material properties of the MOF-5, stainless steel container and hydrogen gas used in this work. †Density of hydrogen is computed either using ideal gas law or virial equation of state (Appendix B). *Temperature-dependent material properties are given in Appendix C. MOF-5

SS

Hydrogen

Density (𝜌, kg/m3)

130

7830

Eqn. of state†

Specific heat* (𝑐𝑝 , J/kg. K)

780

468

14700

Thermal conductivity* (𝜆, W/m.K)

0.088

13

0.206

Viscosity* (Pa-s)

_

_

8.41e-6

Bed porosity, 𝜀𝑏

0.246

_

Jo

ur na

lP

re

-p

ro

_

of

Property

Table 2: Summary of typical mass, momentum and energy balance equations used in recent literature on 2D and 3D modeling of adsorbent bed. The model equations are rearranged to the standard form and presented with consistent notations used in this paper.

𝜕 𝜖𝜌 𝜕𝑡 𝑡 𝑔

𝑣=−

0]

2∇𝑃 𝛼 + √𝛼 2 + 4𝛽𝑐∇𝑃

𝐶𝑝𝑔

𝜕 𝑇(𝜖𝑡 𝜌𝑔 + 𝜌𝑏 𝑞) 𝜕𝑡 + 𝜌𝑏 𝐶𝑝𝑠

= −𝑀H2 𝑆𝑚

Flat radial and linear axial 𝜕 (𝜖𝜌𝐶𝑝 𝑇 + 𝜌𝑏 𝐶𝑝𝑠 𝑇) 𝜕𝑡 velocity profile ∂ 𝑣 𝜕 + (𝜌𝑔 ) + (𝜌𝑢𝐶𝑝 𝑇) = ∂z 𝜀 𝜕𝑧

𝜕𝜌𝑔 𝜕𝑡

lP

𝑀𝐻2 𝑆𝑚 𝜖

−

ur na

1]

𝜕𝑃 𝜕𝑡

re

𝜕𝜌 𝜕𝑡

−

[3

𝜕𝑇 𝜕𝑡 Solver:

-p

+ (−∆𝐻)𝑆m

= ∇. 𝒟∇𝜌𝑔

𝜕 (𝜌 𝑣 ) + ∇. (𝜌𝑔 𝑣𝑣𝑖 ) 𝜕𝑡 𝑔 𝑖

Jo

+ ∇. (𝜌𝑔 𝑣𝑠 ) = −𝑀𝐻2

𝑆𝑚 𝜀

in

activated

In-house

+𝐶𝑝𝑔 ∇. (𝑇𝑐𝑣) code = ∇. (𝑘eff ∇𝑇) − 𝜖𝑡

5]

NG

carbon

+ ∇. (𝜌𝑔 𝑣)

[2

Remarks

Isotherm: Langmuir

ro

[3

Energy Balance

of

Mass Balance Momentum Balance

H2

in

activated

carbon Solver:

In-house

∇. (Λ∇𝑇) + (−∆𝐻)𝑆𝑚 code 𝜕 𝑅𝑇 (𝜌𝑔 ) 𝜕𝑡 𝑀H2

𝜕 (𝜀𝜌𝑔 𝐸𝑔 + 𝜌𝑏 𝐸𝑠 ) 𝜕𝑡

Isotherm:

Radke-

Prausnitz H2

in

activated

carbon

𝜕𝑝 =− + [∇. (𝜏)]𝑖 𝜕𝑥

+ ∇. (𝑣𝜌𝑔 𝐸𝑔 ) =

𝜇 −( 𝑣 𝑘 2

+ 𝐶2 𝑣 )

∇. 𝜆𝑒 ∇𝑇 + 𝛻. 𝜏. 𝑣 − ∇. 𝑣P + (−∆𝐻)𝑆𝑚

Solver: FLUENT Isotherm: Modified DA

[1

𝜀𝑛𝑒𝑡

7]

𝑘 𝑈 = − ∇𝑃 𝜇

𝜕𝜌𝑔 𝜕𝑡

[𝜀𝑛𝑒𝑡 𝜌𝑔 𝑐𝑣 + 𝑞𝜌𝑏 𝑐𝑣

H2

𝜕𝑇 carbon 𝜕𝑡 Solver: + ∇. 𝐸𝑔 𝑣 difference

= ∇. (𝐷𝑒 ∇𝜌𝑔 )

= −∇. (𝑃𝑣) + ∇(𝑘eff ∇𝑇)

+ 𝑀𝐻2 𝑆𝑚

+ 𝑆𝑚 (−𝛥𝐻)

+ 𝜀𝑏 ∇. [𝜌𝑔 𝑣] = −𝑀𝐻2 𝑆𝑚

𝜕𝑣 𝜇 𝑆0 + ( + )𝑣 𝜕𝑡 𝑘 𝜀𝑏

re

lP

(𝜀𝑏 𝜌𝑔 𝐶𝑝𝑔 + 𝜌𝑏 𝑞𝐶𝑝𝑔 + 𝜌𝑏 𝐶𝑝𝑠 )

ur na

= −𝑀H2 𝑆𝑚

𝑘 𝑣 = 𝜌𝑔 (− 𝛻𝑃) µ

𝜕(𝜀𝑏 𝜌𝑔 ) 𝜕𝑡

[2 9]

Jo

+ ∇. [𝜌𝑔 𝑣]

= −𝑀𝐻2 𝑆𝑚

1

H2 in AX-21

Isotherm: Modified

𝑇 )]

2𝜇 − ∇. [( − ɳ𝑑 ) (∇. 𝑣)𝐼] 3 [2 𝜕(𝜀𝑏 𝜌𝑔 ) 𝜕𝑡 2] + ∇. [𝜌𝑔 𝑣] =

DA

Solver: COMSOL

= −∇𝑝 + ∇. [𝜇(∇𝑣 + ∇𝑣

Isotherm: Modified

of

Equation is given below4

Finite

ro

4]

𝜌𝑔

-p

𝜕(𝜀𝑏 𝜌𝑔 ) 𝜕𝑡

activated

+ 𝜌𝑏 𝑐𝑠 ]

+ ∇. 𝜌𝑔 𝑣

[1

in

𝜕𝑇 𝜕𝑡

DA

H2 in MOF-5 bed Solver: COMSOL Isotherm: Modified

+𝜌𝑔 𝐶𝑝𝑔 𝑣. 𝛻𝑇

DA

= ∇. (𝑘eff ∇𝑇) + 𝑄

𝜕(𝜌𝑔 𝑣) + ∇. (𝜌𝑔 𝑣𝑣) 𝜕𝑡

𝜕(𝜀𝑏 𝜌𝑔 𝐸𝑔 + 𝜌𝑏 𝐸𝑠 ) 𝜕𝑡

= −∇𝑝

+ ∇. 𝑣𝜌𝑔 𝐸𝑔

+ ∇. 𝜏

= ∇. (𝑘eff ∇𝑇) − ∇. 𝑣𝑃

𝜇 1 +𝜌𝑔 𝑔 − ( 𝑣 + 𝐶2 𝜌𝑔 𝑣 2 ) 𝑘 2

+ ∇. (𝜏. 𝑣) + (−Δ𝐻)𝑆𝑚

H2

in

activated

carbon Solver: FLUENT Isotherm: Modified DA

The equations are rearranged to use the adsorption source term, 𝑆𝑚 , given by 𝑆𝑚 =

𝜌𝑏

𝜕𝑛𝑎 𝜕𝑡

molH2 /m3bed . Also note that (1 − 𝜀𝑏 )𝜌𝑠

𝜕𝑛𝑎 𝜕𝑡

is replaced as 𝑆𝑚 .

𝜕𝑞

2

Adsorption source term in kg H2 /(m3bed . 𝑠) is written as 𝑀H2 𝑆𝑚 = 𝜌𝑏 𝜕𝑡

3

For microporous adsorbent, Momen et al. [17] defined the net porosity as the difference between

total and micro-porosity 𝜀𝑛𝑒𝑡 = 𝜀𝑡 − 𝜀𝜇 4

𝜕𝑇

𝑇 𝜕𝜌𝑔 𝜕𝑃 ( 𝜕𝑡 𝑔 𝜕𝑇

𝜕𝑇

𝜌𝑔 𝑐𝑝𝑔 (𝜀 𝜕𝑡 + 𝑣. 𝛻𝑇) + 𝜌𝑏 𝑐𝑝𝑠 𝜕𝑡 = ∇. 𝑘∇𝑇 − 𝜌 𝜕∆𝑈𝑎 𝜕𝑡

+

𝜕(𝑢0 𝑛𝑎 ) )− 𝜕𝑡

ℎ𝑆𝑚

Jo

ur na

lP

re

-p

ro

of

ɳ𝑑 ) (∇. 𝑣)𝐼] : ∇𝑣 + 𝜌𝑏 (

2𝜇

+ 𝑣. ∇𝑃) + 𝜇 [(∇𝑣 + ∇𝑣 𝑇 ) − ( 3 −

Table 3: Parameters for Unilan, modified D-A and Tóth isotherm models for adsorption of hydrogen on MOF-5 pellets, re-fitted to the data of Zhou et al. [32] in the temperature range of 77 K to 300 K, and pressures of 1 to 60 bar.

Unilan

Modified D-A

Toth

Value

Parameter

Value

Parameter

Value

nmax (mol/kg)

58

nmax (mol/kg)

151.8

nmax (mol/kg)

65.76

Emax (J/mol)

4834

𝛼 (J/mol)

1941

𝑛0

0.4609

Emin (J/mol)

1604

𝛽 (J/mol. K)

19.2

𝑞𝑇 (J/mol)

3032.1

∆S (J/mol. K)

-46.07

𝑃0 (MPa)

1246

𝑏0 (bar-1)

5.70e-3

𝑚

2

α

-p

ro

of

Parameter

5.08e-3

1.348

𝑣𝑎 (cm3/g)

2.84

𝑣𝑎 (cm3/g)

2.151

RMSE

0.386

RMSE

0.678

RMSE

1.45

Jo

ur na

lP

re

𝑣𝑎 (cm3/g)

Table B.1: Virial parameters for compressibility factor of H2 𝑋 = 𝐵𝑧

𝑋 = 𝐶𝑧

𝑋 = 𝐷𝑧

–1.1501e-4

𝑐𝑧1

4.4713e-6

𝑑𝑧1

–6.4090e-8

𝑏𝑧2

9.2824e-3

𝑐𝑧2

–7.2264e-5

𝑑𝑧2

1.0079e-6

𝑏𝑧3

–2.5577e-2

𝑐𝑧3

3.4550e-4

𝑑𝑧3

–4.9093e-6

𝑏𝑧4

–1.0695e-2

𝑐𝑧4

–3.2069e-4

𝑑𝑧4

7.7846e-6

Jo

ur na

lP

re

-p

ro

of

𝑏𝑧1

Table C.1: Parameters for computing the specific heat (in J/ (kg-K)) of hydrogen 𝑋 = 𝐴𝑐

𝑋 = 𝐵𝑐

𝑋 = 𝐶𝑐

𝑋 = 𝐷𝑐

1.9503e4

𝑏𝑐1

–1.1655e6

𝑐𝑐1

3.9207e7

𝑑𝑐1

-5.7349

𝑎𝑐2

–8.3051

𝑏𝑐2

6.5414e2

𝑐𝑐2

3.233e5

𝑑𝑐2

1.8313e-2

𝑎𝑐3

2.8837

𝑏𝑐3

–4.8813e2

𝑐𝑐3

2.4894e4

𝑑𝑐3

-5.2438e-3

𝑎𝑐4

–3.4213e–2

𝑏𝑐4

5.8284

𝑐𝑐4

–3.0842e2

𝑑𝑐4

6.2135e-5

Jo

ur na

lP

re

-p

ro

of

𝑎𝑐1

Table C.2: Parameters for computing viscosity of hydrogen (in Pa-s) as a function of temperature (K) and pressure (bar). 𝑋 = 𝐴𝑣

𝑋 = 𝐵𝑣

𝑋 = 𝐶𝑣

–9.6555e-6

𝑏𝑣1

8.7136e-6

𝑐𝑣1

3.6370e-6

𝑎𝑣2

1.1302e-8

𝑏𝑣2

–4.5954e-9

𝑐𝑣2

4.3738e-9

𝑎𝑣3

–6.7193e-10

𝑏𝑣3

2.7575e-10

𝑐𝑣3

9.7752e-10

Jo

ur na

lP

re

-p

ro

of

𝑎𝑣1

Table C.3: Coefficients for temperature-dependent specific heat[23] of MOF-5 and steel wall Material

𝑎0

𝑎1 –8.885e–3

0.524

Steel wall

2.24567 –1.71591

𝑎3

𝑎4

9.624e–5

–3.469e–7

4.417e–10

0.11268

0

4.4132e–6

𝑎5

8.74117e–9

𝑎6

6.69455e–12

Jo

ur na

lP

re

-p

ro

of

MOF-5

𝑎2