Authors: Ahmadreza Rahbari, Julio C. Garcia-Navarro, Mahinder Ramdin, Leo J. P. van den Broeke, Othonas A. Moultos, David Dubbeldam, Thijs J. H. Vlugt
Categories: Article
Source: Journal of Chemical and Engineering Data
Thermodynamic Properties of Compressed Hydrogen
Force field-based molecular simulations were used to calculate thermal expansivities, heat capacities, and Joule–Thomson coefficients of binary (standard) hydrogen–water mixtures for temperatures between 366.15 and 423.15 K and pressures between 50 and 1000 bar. The mole fraction of water in saturated hydrogen–water mixtures in the gas phase ranges from 0.004 to 0.138. The same properties were calculated for pure hydrogen at 323.15 K and pressures between 100 and 1000 bar. Simulations were performed using the TIP3P and a modified TIP4P force field for water and the Marx, Vrabec, Cracknell, Buch, and Hirschfelder force fields for hydrogen. The vapor–liquid equilibria of hydrogen–water mixtures were calculated along the melting line of ice Ih, corresponding to temperatures between 264.21 and 272.4 K, using the TIP3P force field for water and the Marx force field for hydrogen. In this temperature range, the solubilities and the chemical potentials of hydrogen and water were obtained. Based on the computed solubility data of hydrogen in water, the freezing-point depression of water was computed ranging from 264.21 to 272.4 K. The modified TIP4P and Marx force fields were used to improve the solubility calculations of hydrogen–water mixtures reported in our previous study [Rahbari A.;et al. J. Chem. Eng. Data 2019, 64, 4103−4115] for temperatures between 323 and 423 K and pressures ranging from 100 to 1000 bar. The chemical potentials of ice Ih were calculated as a function of pressure between 100 and 1000 bar, along the melting line for temperatures between 264.21 and 272.4 K, using the IAPWS equation of state for ice Ih. We show that at low pressures, the presence of water has a large effect on the thermodynamic properties of compressed hydrogen. Our conclusions may have consequences for the energetics of a hydrogen refueling station using electrochemical hydrogen compressors.
To supply the energy demand for a growing population worldwide
and to reduce carbon emissions, efforts are being made to switch from
fossil-based energy production to sustainable energy production using
renewable energy sources.^1−5^ Current estimations show that in 2050 ca. 75% of
the global energy production will be based on fossil fuels.^1^ In recent years, hydrogen is considered as one
of the most promising alternatives in the energy transition for combating
global warming and reducing the carbon footprint from anthropogenic
activities.^4,6−11^ Hydrogen is one of the most abundant elements found on earth mainly
in water and organic compounds.^3^ To produce
hydrogen as a free gas, energy is required.^2,12^ Hydrogen
can be a carbon-free renewable energy source depending on the production
pathway,^12^e.g., green
hydrogen can be produced from water electrolysis in which the required
energy is supplied from renewable sources. The high energy density
of hydrogen per unit of mass (higher heating value of 141.8 MJ/kg
at 298 K)^12^ makes it a flexible off-peak
energy carrier, especially for storing intermittent renewable energy
at peak times, or as a fuel for transportation.^6,13−15^ To date, most of the required hydrogen in refineries
is produced from steam methane reforming (SMR) with CO2 as a byproduct.^2,4,8,12,16^ About 50%
of the produced hydrogen worldwide is used in ammonia synthesis plants.
Hydrogen is also used as a hydrocarbon upgrader or feedstock in refineries, e.g., reducing nitrogen oxides, sulfur oxides, and other
particles causing smog.^6,17^
The heating value of hydrogen (per mass) is higher than that of most common fossil fuels.^2,12^ However, due to the low volumetric density of hydrogen at standard temperature and pressure, 273.15 K and 1 bar (10.79 MJ/m^3^),^9^ and to save space, gaseous hydrogen should be compressed, liquefied, or absorbed by metal hydrides or nanoporous materials such as metal–organic frameworks (MOFs).^2,3,7,17−19^ Compressing hydrogen increases the heating value (per unit of volume) compared to the heating value of conventional liquid fossil fuels.^12^ Compressed hydrogen can be used in fuel-cell electric vehicles or hydrogen-fueled combustion engines.^1,17,20^ To facilitate and increase the use of hydrogen-driven vehicles, hydrogen refueling stations (HRSs) should provide compressed hydrogen to fuel-cell electric vehicles or internal combustion hydrogen engines.^21^ Current vehicle technologies allow direct on-board hydrogen storage with pressures up to 700 bar inside the tank.^9,18,20,22,23^ This allows a driving range of at least 500 km.^1,9,24−26^ To meet the target refueling time (3–4 min),^22,23^ a large pressure difference is required when filling the hydrogen tank, e.g., the pressure of compressed hydrogen before entering the vehicle should be as high as 875 bar.^21^ Therefore, energy- and cost-efficient hydrogen compression technologies are required for operating refueling stations at these high pressures.
Hydrogen compression technologies can be divided into mechanical and nonmechanical categories.^9^ There are well-known limitations associated with mechanical compression such as high mechanical losses, high maintenance, high levels of noise, etc.^20^ In addition, the compression work can be as high as 30% of the energy of the stored hydrogen,^27^ and the efficiency is low when operating at small scales.^9,28^ Alternatives to mechanical compression are (1) cryogenic compression, (2) metal hydride or thermally powered compressors,^29^ and (3) electrochemical hydrogen compressors (EHCs), also known as electrochemical hydrogen pumps.^9,20,30^ A detailed overview of conventional and state-of-the-art technologies for hydrogen compression technologies is provided in ref (9).
Recently, EHCs have
become cost-competitive alternatives to their
mechanical counterparts.^16,30,31^ Due to the advantages over mechanical compression, EHCs are on the
way to be commercialized on a large scale. EHCs can generate a high
output pressure (up to 1000 bar), which is suitable for refueling
hydrogen vehicles.^20,30^ EHCs have a higher efficiency
compared to mechanical compression, especially for low volumetric
flows.^9,32^ Since there are no moving parts in an EHC,
there are no maintenance costs associated with rotary parts, and the
compressor operates silently. The working principle of an EHC is based
on the proton-exchange membrane fuel-cell technology.^4,9,30^ A schematic representation of
the working principle of an EHC is shown in Figure 1. On the anode side, electricity is used
to split the hydrogen molecule into protons and electrons, and the
protons pass through the membrane. Hydrated protons (H2n+1O~n^+^ in which n is the electro-osmosis coefficient) are transferred to
the cathode side.^4,9,13^ At
the cathode side, the protons are combined with electrons to form
hydrogen molecules.^20^ As the number of
hydrogen molecules increases on the cathode side, the discharge pressure
at the outlet of the EHC increases. In EHCs, the driving force for
compression is the electrical potential difference between the anode
and the cathode. This means that the discharge pressure can be as
large as allowed by material limitations. It is important to note
that generating a large differential pressure over a single compression
cell lowers the efficiency of the EHC and may cause back-diffusion
of hydrogen from the cathode to the anode side.^9,20^ To
circumvent this issue, one can raise the discharge pressure using
a cascaded arrangement of single compression cells.^30,33^ The highest differential pressure reported in the literature for
a single compression cell is 168 bar.^20^ Two types of proton conducting membranes are typically used for
electrochemical ^16^^16^ (1) proton electrolyte membranes such as Nafion-type
membranes,^16,28,34^ which are suitable for low-temperature applications; and (2) protonic-ceramic
electrolyte membranes, which are suitable for high-temperature applications
(e.g., 600 °C).^16^ These membranes are highly selective to permeation of protons as
impurities cannot pass the membrane.^33^ The
proton-exchange membrane requires hydration for its functionality.
Therefore, it is important that the membrane inside the compressor
remains hydrated with water when operating. While compressing the
hydrogen stream, the EHC can simultaneously purify the hydrogen stream
such that the compressed hydrogen at the outlet is saturated with
water, and no other impurities are present. Recently, HyET BV^30,31^ has developed state-of-the-art EHCs operating based on the proton-exchange
membrane technology. Currently, each compression stack can compress
50 kg H2~/day up to 1000 bar. The EHC capacity can be scaled
up using several compression stacks.^31,33^ The EHC by
HyET BV can achieve the same pressure difference as in mechanical
compression with fewer steps and silently. The compressed hydrogen
is saturated with water and free from other contaminants.
Figure 1 Schematic representation of the working principle of electrochemical compression of hydrogen.^30^ Hydrogen molecules are split at the low-pressure side of the membrane (i.e., using platinum-alloy catalysts). Protons are forced through the membrane and form hydrogen molecules at the high-pressure side.
Here, we consider the application of an EHC to generate compressed hydrogen for refueling a hydrogen vehicle. Although the solubility of water in compressed hydrogen at the outlet of the EHC is very low, this affects the thermodynamic behavior of the gas phase (depending on the water content). As the compressed hydrogen gas enters the on-board storage tank, it follows directly from the first law of thermodynamics of an open system that the temperature of hydrogen inside the tank rises. A similar effect happens for the isenthalpic expansion of hydrogen (e.g., using a throttle) due to the negative Joule–Thomson coefficient of hydrogen at high pressures.^23,35^ To protect the mechanical integrity of the tank, the maximum temperature limit of the hydrogen in the process of filling should not exceed 85 °C.^26,35,36^ To avoid high temperatures inside the hydrogen tank, precooling is required before refueling to bring the temperature of compressed hydrogen down between −33 and −40 °C.^23,26,35^ A complete hydrogen fueling protocol is compiled by the Society of Automotive Engineers (SAE) in “Fueling Protocols for Light Duty Gaseous Hydrogen Surface Vehicles (SAE J2601)”.^22,37^
Due to the presence of water content at the output of an EHC, ice formation may take place during precooling, which may clog the hydrogen delivery system. Experimental data on the solubility and other thermodynamic properties of water–hydrogen systems in the published literature are limited, especially at high pressures.^38,39^ Published experimental data for solubilities of water–hydrogen mixtures in the gas or/and liquid phase at high pressures can be found in refs (40−49). A summary of the experimental data in the tabulated form can be found in ref (38). However, little is known about the influence of water on the thermodynamic properties of hydrogen at high pressures. Performing experiments for obtaining thermodynamic data at high pressures is challenging due to safety protocols, material limitations, high costs, etc.^50−54^ A robust alternative for high-pressure experiments is force-field-based molecular simulation.^55,56^ In this work, we use force-field-based molecular simulation to calculate properties of compressed hydrogen by varying the water content at different temperatures and pressures. Molecular simulation is an approach for calculating macroscopic as well microscopic properties of materials based on the knowledge of their constituent molecules and atoms as well as their interactions. Molecular simulation allows the computation of material properties at high pressures, which would be difficult or expensive to access using experiments. By modeling the VLE of water–hydrogen systems, it is possible to compute the amount of water present in compressed hydrogen at temperatures and pressures corresponding to the freezing point of water. In this way, one can predict water/ice formation when precooling the hydrogen stream. The VLE calculations can also be used to predict the freezing-point depression of water in the presence of hydrogen. In this work, the VLE calculations of hydrogen–water are performed at temperatures and pressures corresponding to the water–ice Ih equilibrium in the pressure range between 100 and 1000 bar. Using molecular simulation, we also calculate thermodynamic properties such as thermal expansivity, heat capacity, and the Joule–Thomson coefficient of the water–hydrogen mixtures.
This paper is organized as follows. Thermodynamic properties of compressed hydrogen–water mixtures are investigated using force-field-based molecular simulation. An overview of the molecular models used in this paper is provided in Table 1. Thermodynamic properties of mixtures are related to the derivatives of extensive properties of the mixtures.^57^ One can calculate the appropriate thermodynamic derivatives (i.e., for calculating heat capacities, thermal expansivities, and Joule–Thomson coefficients) from ensemble fluctuations.^55,56,58^ In Section 2, we provide the appropriate expressions for calculating thermodynamic derivatives to obtain the required thermodynamic properties. In Section 3, we explain how vapor–liquid equilibria of hydrogen–water mixtures along the melting line of ice can be computed. This allows one to estimate how much water would exit the saturated hydrogen stream during deep cooling at low temperatures. To the best of our knowledge, high-pressure experimental data of water–hydrogen systems at temperatures between 264.2 and 272.4 K are not available in the literature. Volumetric data from the literature are used to compute the chemical potential of ice for a large pressure range. Simulation details are provided in Section 4. All simulation results are presented and discussed in Section 5. Our conclusions are summarized in Section 6. The results show that small concentrations of water dissolved in hydrogen in the gas phase can significantly change the properties of the gas phase (i.e., Joule–Thomson coefficient, heat capacity, and thermal expansivity) compared to pure hydrogen at the same conditions. Our results will therefore have consequences for the energetics of a hydrogen refueling station using EHCs including drying or precooling of the compressed hydrogen stream.
of Mixtures Obtained from Ensemble Fluctuations
To compute thermodynamic properties
of compressed hydrogen with
and without traces of water using molecular simulations, we use derivatives
of volume, internal energy, and enthalpy with respect to temperature
and pressure. These derivatives are required to calculate properties
such as thermal expansivity, heat capacity, and the Joule–Thomson
coefficient.^57,59,60^ These thermodynamic derivatives are directly obtained from ensemble
fluctuations at constant composition.^56,58,61,62^ Lagache et al.^58^ showed that the derivative of an extensive property X with respect to β = 1/(kBT) (in which kB is the
Boltzmann constant and T the absolute temperature)
in the NPT ensemble can be obtained from the ensemble
fluctuations as follows
where Ĥ = U + PV is the configurational enthalpy of the system,^58^P is the imposed pressure, and U is the potential energy of the system consisting of an intermolecular contribution U^ext^ and an intramolecular contribution U^int^. The mathematical proof for this is provided in the Supporting Information. In a similar manner, one can obtain an expression for the derivative of the ensemble average of an extensive property X with respect to pressure^58^
The last term on the right-hand side of eq 2 is nonzero if the extensive property X includes a pressure-dependent contribution, for instance, for the configurational enthalpy Ĥ. The term was not explicitly provided in the corresponding equation in ref (58). The derivation of eq 2 is also provided in the Supporting Information. The extensive property X may be replaced by an intensive property x. Using eqs 1 and 2, one can obtain various thermodynamic properties of mixtures from ensemble fluctuations. For instance, combining the definition of thermal expansivity α~P~ with eq 2(58) leads to
To compute the Joule–Thomson coefficient, one needs to compute the heat capacity of the mixture as well. The extensive heat capacity follows from combining the definition of the heat capacity with eq 1
where the enthalpy H also includes the kinetic energy term (K): H = U^int^ + U^ext^ + K + PV = Ĥ + K. Jorgensen et al.^63^ and Lagache et al.^58^ emphasize that the expression for the heat capacity should be split into an ideal contribution and a residual contribution.
Including
the ideal gas contribution of heat
capacity from the force field usually leads to overestimating the
heat capacity.^64,65^ To calculate CP^id^(T) and CP^res^(T, P), Lagache et al.^58^ proposed the following split for the enthalpy of the system
in which N is the number
of molecules in the systems. The ideal gas heat capacity for a pure component can be obtained from
tabulated thermodynamic reference tables.^59,66^ Alternatively, CP^id^ can be calculated by taking into
account the contribution of translational, rotational, vibrational,
and electronic energy levels of an isolated molecule, using quantum
mechanical calculations or experimental data.^67,68^ The residual heat capacity CP^res^ is obtained
by combining eqs 1 and 6
The Joule–Thomson coefficient μJT is obtained
using
Equilibrium of Hydrogen–Water Mixtures along the Melting Line of Ice
of Solubilities in the Liquid and Gas Phases
In molecular simulation, VLEs are conveniently computed by Monte Carlo simulations^55^ in the Gibbs ensemble.^69,70^ In the Gibbs ensemble, the liquid and gas phases are simulated in different simulation boxes. During the simulation, standard Monte Carlo trial moves are used to generate configurations in each phase and to exchange molecules between the phases. In simulations of multicomponent systems in the Gibbs ensemble, at each temperature and pressure, chemical equilibrium is obtained when the chemical potentials of both phases are equal. From the few published experimental VLE data of hydrogen–water mixtures at high pressures,^38,42,44,47−49^ it is known that the solubility of hydrogen in the liquid phase and the solubility of water in the gas phase are very low. The experimental data in refs (40−49) are tabulated in the Supporting Information of ref (38). Simulating the VLE of hydrogen–water mixtures at such dilute concentrations in the Gibbs ensemble^69,70^ would be computationally expensive (requiring large system sizes), and the simulation time would be prohibitively long. It is possible to simulate the VLE of hydrogen–water mixtures using a smaller system size; however, one needs to compensate by performing much longer simulations to obtain reliable statistics.^38^ Despite its reliability and small finite-size effects (except in the vicinity of the critical point),^71,72^ the Gibbs ensemble is not the preferred method here. The liquid phase is almost pure water at high pressures, and the gas phase is almost pure hydrogen.^38^ Therefore, one can simulate the liquid phase and the gas phase independently in the isobaric–isothermal ensemble and obtain the chemical potentials of water and hydrogen at phase equilibrium. Since the concentration of the solute is very low in each phase, the chemical potential of the solvent will not be different from that of the pure solvent.
In Monte Carlo simulations, one of the methods to efficiently compute the chemical potentials, especially in dense systems such as liquid water, is the continuous fractional component Monte Carlo (CFCMC) method.^73−75^ To compute the chemical potentials of water and hydrogen in the liquid and the gas phase, we use the CFCMC in the NPT ensemble.^73,76^ In the CFCMC method, a so-called fractional molecule (one per component type) is added to the system. The interaction potential of the fractional molecule of type i is scaled using a coupling parameter λ~i~ ∈ [0, 1]. By constantly performing a random walk between λ~i~ = 0 and 1 during the simulation, one can mimic gradual insertion or deletion of a molecule and compute the chemical potential μ~i~ of a certain molecule type^77^
where
μ~i^0^(T)
is the reference chemical potential, which contains the intramolecular
contributions (rotation, vibration, and translation) and the bond
dissociation energy of component i.^67^ ρi~ is the number density
of i (i.e., the ratio between the
number of molecules of i and the volume of the system).
ρ0 = 1 Å^–3^ is the reference
number density to make the argument inside the logarithm dimensionless.
The reference chemical potential is important when simulating chemical
reactions, taking into account free energy differences of bond breaking/formation.
In this work, the choice of μ~i^0^(T) is not important
since chemical reactions are not considered. The first and the second
term on the right-hand side of eq 9 are the ideal gas parts of the chemical potential
(μi^id^), and the last term on the right-hand side
of eq 9 is the excess
chemical potential (μi^ex^). The excess chemical potential
is related to the interactions between a molecule and the surrounding
molecules and is independent of the reference chemical potential μi^0^. For an ideal gas, βμi^ex^ = 0. p(λi) is the Boltzmann probability distribution
of λi. The excess chemical potential
can be directly computed using the values of p(λi) at λi~ = 0 and 1 or by performing extrapolation or interpolation; see eq 9.^73,78^ Further details regarding the CFCMC method and recent advances of
this method are provided in ref (73) and references in there.
By including
a fractional molecule of water in the liquid phase,
it is straightforward to calculate the chemical potential of water
in the liquid phase (μH2O(l)).^73,76^ To calculate the excess chemical potential
of water in the gas phase, a fractional molecule of water is added
to the gas phase (containing hydrogen molecules). The solubility of
water in the gas phase is obtained by imposing equal chemical potentials
of water in the gas and liquid phases (μ = μH2O(l)H2O(g)~~)
and solving the ideal gas contribution of the chemical potential of
water in the gas phase
In a similar way, one can compute the solubility
of hydrogen in the liquid phase. Since the solubility of hydrogen
is very low in the liquid phase, the excess chemical potential of
hydrogen can be computed by adding a single fractional molecule of
hydrogen in the liquid phase. The ideal gas chemical potential of
hydrogen (related to the number density) in the liquid phase is initially
unknown. By imposing equal chemical potentials of hydrogen in the
gas phase (almost pure hydrogen) and the liquid phase (μH2(g)~~ = μH2(l)~~), the
ideal gas chemical potential of hydrogen and the number density of
hydrogen in the liquid phase at equilibrium are solved using
Equations 10 and 11 are derived in the Supporting Information.
of the Chemical Potential of the Most Stable Phase
Precooling compressed hydrogen containing water may lead to ice formation during the refueling process for a hydrogen-driven vehicle. To predict ice formation, one can use an empirical equation of state for ice and water or compare the chemical potentials of water and ice at different temperatures and pressures. At equilibrium, the chemical potentials of ice and water are equal, which corresponds to the minimum Gibbs energy of the system. By departing from the equilibrium pressure, the most stable phase (liquid or ice) will be the phase that lowers the Gibbs energy of the system the most.^59,79^
Simulating the molecular structure of ice,^80−84^ solid–liquid equilibria, and computing the chemical potential of ice require special simulation techniques such as anisotropic NPT simulations, the Einstein crystal method,^85^ the Einstein molecule approach,^86^ and/or the Parrinello–Rahman-type sampling.^83^ For details on free energy calculations of ice, the reader is referred to the works of Vega et al.^81,86−88^ An alternative and a simpler pathway for computing the chemical potential of ice is to compute the chemical potential of liquid water along the experimental solid–liquid equilibrium. It is straightforward to compute the chemical potential of liquid water from simulations in the CFCNPT ensemble.^73,76^ Alternatively, thermodynamic tables, or the empirical equation of state for ice and water,^89,90^ can be used to compute the chemical potentials of ice/water along the melting line. In this work, we use both molecular simulation of liquid water and the IAPWS equation of state for water,^91^ which is implemented in REFPROP.^92^ The equations of state used in REFPROP for pure hydrogen, liquid water, and ice are fitted to experimental data. Since the compressibility of both water and ice is low, it is possible to calculate the chemical potentials of liquid water and ice at different melting temperatures by varying the pressure, using
where Tm is the
melting temperature of ice. The notation for melting pressure Pm(T) is used to emphasize that
the melting pressure is a function of the melting temperature. From eqs 12 and 13, it is clear that knowledge of the molar volume of ice and
water is required to calculate the chemical potential of the most
stable phase at different pressures for each melting temperature.
Marion et al.^93^ reported experimental values
for the molar volume υice of pure ice at P^0^ = 1 atm, without cracks or gas bubbles, based
on crystal lattice parameters for ice Ih^94^
where T is in units of kelvin. To calculate the molar volume of ice at P ≠ P^0^, the isothermal expansivity of ice is required. Marion et al.^93^ estimated the isothermal compressibility β~T~ of ice using
where T is in units of kelvin, and the molar volume of ice at pressure P follows from
The data for molar volume of liquid water are available in REFPROP.^91,92^ By combining eqs 12 and 14–16, one can estimate the chemical potentials of ice and water at different pressures as a function of melting temperature.
Depression of Water with Dissolved Hydrogen
The solubility of hydrogen in liquid water affects the freezing temperature of ice. In the Supporting Information, we show that the freezing-point depression is obtained from^57^
where Δh̅fus =h̅ice – h̅water is the enthalpy
of fusion of ice, R is the gas constant, and xH2~~ is the mole fraction of hydrogen
dissolved in liquid
water. In eq 17, it
is assumed that the solution is ideal and that the solvent is nonvolatile.
The difference between the enthalpy of ice and water is a function
of temperature and does not vary significantly with pressure.^95^ Osborne and Dickinson^95^ performed an experiment to
obtain the enthalpy of fusion of ice, which is Δh̅fus = 6010.44 J/mol. Osborne and Dickinson^95^ also tabulated the difference in enthalpy between
water and ice at different temperatures. As a service to the reader,
in Table S1 of the Supporting Information,
this enthalpy table is provided.
Details
In this work, two types of simulations are simulations with fractional molecules and simulations without fractional molecules. To calculate the chemical potentials of water or hydrogen from simulations, fractional molecules are used. For computation of thermodynamic derivatives from ensemble fluctuations, no fractional molecules are used. To calculate the thermodynamic derivatives of water–hydrogen mixtures from ensemble fluctuations, Monte Carlo simulations are performed in the NPT ensemble, as explained in Section 2. To compute the chemical potentials of water and hydrogen in the gas and liquid phases, independent simulations are performed by expanding the NPT ensemble with fractional molecules of water and hydrogen. All simulations are performed using open-source Brick-CFCMC software.^96^ All molecules are considered rigid objects. Intermolecular interactions consist of pairwise Lennard-Jones (LJ) and Coulomb interactions. Force-field parameters for hydrogen and water models are provided in the Supporting Information of ref (38). The interaction parameters of the modified TIP4P force field are provided in the Supporting Information. The Ewald summation with a precision of 10^–6^ is used to handle electrostatic interactions. Lennard-Jones interactions are truncated at 12 Å, and analytic tail corrections are applied. Periodic boundary conditions are used. The Lorentz–Berthelot mixing rules are applied for cross-interactions between different interaction sites.^55,56^ For the temperature range considered in this study, no quantum corrections are needed for hydrogen.^97,98^
For all simulations of pure water and pure hydrogen, 730 molecules of water and 600 molecules of hydrogen are used. The mixture compositions of the water–hydrogen mixtures are based on the VLE calculations of ref (38), using the Marx^99^ force field for hydrogen and the TIP3P^100^ force field for water. All mixture compositions are provided in Table 2. In the VLE calculations of water–hydrogen systems, simulations using the TIP3P force field predict the composition of the gas phase in good agreement with available experimental data.^38^ Also, compared to other commonly used force fields, the chemical potential of the TIP3P water is in good agreement with the chemical potential of water obtained from the IAPWS equation of state.^91^ It is important to note that while the performance of the TIP3P force field is not the best overall among the different force fields of water,^101^ it performs reasonably well for calculating the chemical potential of water. As shown in ref (38), this property is crucial for VLE calculations of water–hydrogen systems. An in-depth discussion on the role of the chemical potential of water can also be found in Section S5 of the Supporting Information. The hydrogen force field by Vrabec et al.^102^ was also used to investigate the thermodynamic properties of water–hydrogen mixtures.^102^
For simulations including fractional molecules, equilibration
runs
were performed between 2 × 10^5^ and 5 × 10^5^ cycles until the optimum weight function is obtained. In
each cycle, k trial moves are performed to generate
new configurations (to be accepted or rejected based on the Metropolis
acceptance rules^103^). In Brick-CFCMC software, k is set as the maximum between the number of molecules
in the system and the number 20.^96^ In the
CFCMC algorithm, the scaling of the interactions of the fractional
molecules is performed as follows. Starting from the value of λ
= 1, first the electrostatic interactions are switched off. This is
followed by switching off the van der Waals interactions. Details
about scaling the interactions of the fractional molecule can be found
in the manual of Brick-CFCMC.^96^ A weight
function W(λ~i)
was used in the sampling scheme to ensure that the observed probability
distribution piobs(λ) is approximately flat. The Boltzmann probability
distribution p(λi) is recovered by p(λi) ∼ exp[−W(λi)] × pi~). The weight function is computed using
the Wang–Landau algorithm^104^ followed
by an iterative scheme. Applying a weight function results in a simulation
where all different λ states are visited with the same probability.
This improves the accuracy of the computed chemical potentials.^105^ The probabilities p(λ
= 1) and p(λ = 0) are directly sampled using
the scheme in ref (105) without using interpolation or extrapolation. For each temperature
and pressure, 10^7^ production runs are performed, and the
uncertainties of ensemble averages are calculated using block averaging
(five blocks). Simulations in the CFCNPT ensemble
are performed using trial moves with the following
1% volume changes, 35% translations, 30% rotations, 17% λ changes,
8.5% random reinsertions, and 8.5% identity changes. The random reinsertions
and identity changes are combined in a hybrid trial move.^76,96^obs(λ
For simulations of hydrogen–water mixtures without fractional molecules, 2 × 10^5^ equilibration cycles are performed followed by 8 × 10^6^ production cycles. Simulations in the NPT ensemble are performed using trial moves with the following 1% volume changes, 49.5% translations, and 49.5% rotations. The uncertainties of ensemble averages are calculated using block averaging (five blocks). The confidence level of all reported uncertainties in this work is 95%.
and Discussion
of Pure Hydrogen at High Pressures
Thermodynamic properties of pure hydrogen
were computed using different force fields. The Cracknell,^97^ Buch,^106^ Hirschfelder,^107^ Vrabec,^102^ and Marx^99^ force fields for hydrogen were used to compute
the thermal expansivity α~P, specific
heat capacity, cP, and
the Joule–Thomson coefficient μJT~. For pure
hydrogen, these thermodynamic properties were also obtained from an
empirical equation of state for hydrogen in REFPROP.^92^ The results are shown in Figure 2, and the raw data are provided in Table 3. The Vrabec force
field predicts α~P~ and cP, in good agreement with the empirical
equation of state obtained from REFPROP.^92^ This quantitative agreement is consistent with the fact that the
hydrogen force field by Vrabec et al. is optimized for thermal properties
of hydrogen and the speed of sound.^102^ None
of the force fields could very accurately predict the Joule–Thompson
coefficient from experimental data; however, the force fields can
capture the behavior of the Joule–Thomson coefficient qualitatively, i.e., the Joule–Thomson coefficient decreases with
the pressure as it is the case for the Joule–Thomson coefficient
obtained from experimental values. The results from the Marx and Vrabec
force fields show the least deviation. The Joule–Thomson coefficient
of TIP3P water is also calculated, and the results are compared to
those obtained from REFPROP.^92^ For pure
water, the simulations were performed along the VLE line at pressures
between 10 and 100 bar using experimental densities. The results are
shown in Figure S1 of the Supporting Information.
Excellent agreement is found between the experiments and the simulations.
The Marx and Vrabec force fields, together with the TIP3P force field,
are used in this work to further investigate thermal properties of
hydrogen–water mixtures.
Figure 2 Comparison of different force fields to predict (a) thermal expansivity, (b) constant pressure heat capacity, and (c) the JT coefficient of pure hydrogen in the gas phase at T = 323 K and pressures ranging between P = 10 and 1000 bar. Experimental data from REFPROP^92^ (lines), molecular force Cracknell^97^ (circles), Buch^106^ (upward-pointing triangles), Hirschfelder^107^ (squares), Vrabec^102^ (downward-pointing triangles), and Marx^99^ (diamonds). Raw simulation data are provided in Table 3. Error bars are smaller than symbols.
Mixtures at Low Water Concentrations in the Gas Phase
In Figure 3, the thermal expansivities
α~T~ are shown for force-field combinations
TIP3P–Marx and TIP3P–Vrabec at temperatures T = 366 and 423 K. The composition of the mixture at each
temperature and pressure, provided in Table 2, is based on the solubility calculations
of ref (38). The mixture
compositions, and raw data of Figure 3, obtained from molecular simulations are provided
in Tables 4 and 5. Thermodynamic properties of hydrogen–water
mixtures in Table 4 computed using force-field combinations TIP3P–Marx and TIP3P–Vrabec
were used for thermodynamic properties provided in Table 5. Thermodynamic properties of
pure hydrogen in the gas phase obtained from REFPROP are provided
in Table S2 of the Supporting Information.
To better see the influence of the water content on the properties
of compressed hydrogen, the corresponding properties of pure hydrogen
from simulations and empirical data are also provided in this figure.
The simulation results show that the values of thermal expansivity
of the mixture are very similar to the behavior of pure hydrogen,
with deviations of ca. 1% or less for pressures higher
than 50 bar, at T = 366 K, and for pressures higher
than 800 bar, at T = 423 K. In both cases, the mole
fraction of water is ca. yH2O = 0.01. At pressures below 50 bar, the influence of
water content on the properties of the mixture becomes significant.
Figure 3 Thermal expansivities of pure hydrogen in the gas phase and hydrogen in the gas phase that is saturated with water. Triangles, pure hydrogen; circles, hydrogen–water mixtures; lines, empirical data for pure hydrogen from REFPROP.^92^ Thermal expansivities of mixtures computed using the TIP3P^100^ and Marx^99^ force fields are shown in subfigures (a) and (b) for T = 366 and 423 K, respectively. Thermal expansivities of mixtures computed using the TIP3P^100^ and Vrabec^102^ force fields are shown in subfigures (c) and (d) for T = 366 and 423 K, respectively. The composition in the gas phase is obtained from VLE simulations of the water–hydrogen mixture at T = 423 and 366 K. Raw data obtained from molecular simulations (including the composition in the gas phase) are provided in Tables 4 and 5. Thermodynamic properties of pure hydrogen in the gas phase, obtained from REFPROP, are provided in Table S2 of the Supporting Information. Error bars are smaller than symbols.
A similar behavior is observed for the heat capacity of hydrogen–water
mixtures at T = 366 and 423 K; see Figure 4. Compared to pure hydrogen,
the heat capacity of hydrogen–water mixtures deviates around
1% or less for pressures higher than 100 bar, at T = 366 K, and for pressures around 800 bar, at T = 423 K. As can be seen in Figure 4, the presence of water at pressures lower than 50
bar significantly changes the heat capacity of the mixture compared
to pure hydrogen. For instance, based on the raw data in Table 4, one can observe
that the presence of around 13% water (molar) in the gas phase results
in a 26% deviation of the heat capacity at P = 50
bar and T = 423 K. A similar observation is made
based on the raw data in Table 5. Based on eq 8, it is expected that the Joule–Thomson coefficient shows
a similar behavior as the thermal expansivity and heat capacity. In Figure 5, it can be seen
that at around 9% water content in the gas phase (P = 80 bar), the deviation of the Joule–Thomson coefficient
is around 255%, at T = 423 K. Based on the raw data
in Table 5, the deviation
of the Joule–Thomson coefficient is around 189%. Due to the
presence of water, the Joule–Thomson coefficient has a positive
sign for pressures below 50 bar. The deviation from the Joule–Thomson
coefficient of pure hydrogen is reduced to 2% at P = 1000 bar and T = 423 K, where the mole fraction
of water is around yH2O = 0.008.
At T = 366 K and P = 300 bar, the
mole fraction of water is around 0.004 and the Joule–Thomson
coefficient of the mixture deviates around 3.6% from the Joule–Thomson
coefficient of pure hydrogen. It can be observed that at relatively
low pressures, the effect of water on the properties of the gas phase
(i.e., the Joule–Thomson coefficient, heat
capacity, and thermal expansivity) is significant. However, at pressures
higher than 700 bar, the effect of water content on the properties
of hydrogen gas is almost negligible. Since the Joule–Thomson
coefficient is also a function of molar volume, additional plots of
the Joule–Thomson coefficients as a function of molar volume
of the mixture are shown in Figure S2 of
the Supporting Information.
Figure 4 Heat capacities of hydrogen in the gas phase and hydrogen in the gas phase that is saturated with water. Triangles, pure hydrogen; circles, hydrogen–water mixtures; lines, empirical data for pure hydrogen from REFPROP.^92^ Heat capacities of mixtures computed using the TIP3P^100^ and Marx^99^ force fields are shown in subfigures (a) and (b) for T = 366 and 423 K, respectively. Heat capacities of mixtures computed using the TIP3P^100^ and Vrabec^102^ force fields are shown in subfigures (c) and (d) for T = 366 and 423 K, respectively. The composition in the gas phase is obtained from VLE simulations of the water–hydrogen mixture at T = 423 and 366 K. Raw data obtained from molecular simulations (including the composition in the gas phase) are provided in Tables 4 and 5. Thermodynamic properties of pure hydrogen in the gas phase, obtained from REFPROP, are provided in Table S2 of the Supporting Information. Error bars are smaller than symbols.
Figure 5 Joule–Thomson coefficients of hydrogen in the gas phase and hydrogen that is saturated with water. Triangles, pure hydrogen; circles, hydrogen–water mixtures; lines, empirical data for pure hydrogen obtained from REFPROP.^92^ Joule–Thomson coefficients of mixtures computed using the TIP3P^100^ and Marx^99^ force fields are shown in subfigures (a) and (b) for T = 366 and 423 K, respectively. Joule–Thomson coefficients of mixtures computed using the TIP3P^100^ and Vrabec^102^ force fields are shown in subfigures (c) and (d) for T = 366 and 423 K, respectively. The composition in the gas phase is obtained from VLE simulations of the water–hydrogen mixture at T = 423 and 366 K. Raw data obtained from molecular simulations (including the composition in the gas phase) are provided in Tables 4 and 5. Thermodynamic properties of pure hydrogen in the gas phase, obtained from REFPROP, are provided in Table S2 of the Supporting Information. Error bars are smaller than symbols.
at the Melting Temperature
In Figure 6, the phase boundary diagram of water and
ice Ih is shown. The corresponding data for this figure are obtained
from experimental data.^91,92,108^ The melting temperatures and pressures in Figure 6 are used for calculating chemical potentials
of ice and water, from molecular simulation and empirical data. The
IAPWS equation of state in REFPROP^91^ is
used to calculate the chemical potential of ice and water along the
melting line (Tm, Pm(Tm)). The results are shown in Figure 7, and the raw data
are provided in Table 6. For hydrogen, excellent agreement is observed between the chemical
potentials of the Marx force field and the chemical potentials obtained
from REFPROP. The chemical potential computed with the TIP3P force
field differs by a maximum of ca. 1 kJ/mol compared
to the chemical potential obtained from the empirical equation of
state of water.^91^ This difference does
not cause significant deviation between the computed solubilities
and the experimental solubilities of water in compressed hydrogen,
as was shown in ref (38).
Figure 6 Phase boundary diagram of water and ice Ih taken from experimental literature data.^91,92,108^ The open circle indicates the solid–liquid–vapor triple point: P = 611.657 Pa and T = 273.16 K.^92^
Figure 7 Chemical potentials of pure hydrogen and water along the melting line of ice (T
m, Pm) computed from CFCMC simulations using the TIP3P^100^ and Marx^99^ force fields. Lines indicate the chemical potentials calculated from the empirical equation of states using REFPROP.^92^ Raw data are provided in Table 6. Error bars are smaller than symbols.
By changing the pressure, the equilibrium between ice and water shifts toward the most stable phase depending on the pressure. At different melting temperatures of ice, one can compute the chemical potentials of the most stable phase (liquid water or ice) as a function of pressure using eqs 12 and 14–16. This allows one to calculate the chemical potential of ice, based on empirical data without explicitly simulating the crystal structure of ice. The chemical potential of the most stable phase at each temperature can be calculated as a function of pressure using eqs 12 and 13. An example is provided in Figure 8 where the melting point of ice is at T = 266.73 K and P = 750 bar. It can be seen in Figure 8 that by increasing the pressure, the chemical potential of water would be lower than that of ice, which is the most stable phase for P > 750 bar. Decreasing the pressure would lead to having ice as the most stable phase.
Figure 8 Chemical potentials of ice and water at T
m= 266.73 K as a function of pressure. Lines, crystalline ice; dashed lines, liquid phase. At this temperature, ice and water are in equilibrium at Pm= 750 bar. At P > Pm, the chemical potential of water is smaller than the chemical potential of ice, which indicates the most stable phase for that pressure range. For P < Pm, ice is the most stable phase. Raw data are provided in Table S3 of the Supporting Information, including different melting temperatures 264.21 K < Tm< 272.4 K. Equations 12 and 14–16 and the IAPWS equation of state for water^91^ were used to generate the data.
Mixtures at Freezing Temperature of Water
The VLE calculations of water–hydrogen mixtures are performed along the melting line of ice; see Figure 6. The corresponding melting temperatures and pressures for ice Ih are directly obtained from REFPROP.^91,92^ One can obtain the same melting curve from ref (108). The melting temperatures selected in this work are between T = 272.4 and 264.21 K, which correspond to pressures between P = 100 and 1000 bar. It is crucial to consider the phase diagram of Figure 6 in precooling of the compressed hydrogen mixed with water (from the output of an EHC) and to calculate the amount of water dissolved in hydrogen at low temperatures. By performing independent simulations of the liquid phase (water) and the gas phase (hydrogen), the chemical potentials of water and hydrogen are computed along the melting line of ice. In Table 7, we provide an example to emphasize that the chemical potential of the solvent (water) is not affected by the chemical potential of the solute (hydrogen) at low concentrations. Due to the application of the EHC, the pressure range is limited to the range between P = 100 and 1000 bar. Using eqs 10 and 11, the VLEs of water–hydrogen mixtures were computed at the melting temperature of ice as a function of pressure. At each melting temperature, different pressures were considered for the VLE calculations (while remaining in the liquid water region). The results are shown in Figure 9. The raw data corresponding to the liquid phase are provided in Table 8, and the raw data corresponding to the gas phase are provided in Table 9. In Figure 9, only simulation results are provided, as to the best of our knowledge, no experimental VLE data of water and hydrogen, along the melting line of ice, are reported in the published literature. From simulations, the solubility of water is the highest at ca. 90 ppm (molar) at 272.4 K (and 100 bar) and the lowest at ca. 7 ppm (molar) at 264.21 K (and 1000 bar). Such low solubilities justify our choice of not using the Gibbs ensemble. The solubilities of hydrogen in the liquid phase vary between 0.01 and 0.001, corresponding to 272.4 and 264.21 K, respectively. It is well-known that to obtain solubilities in good agreement with experimental data, the accuracy of the chemical potentials of water and hydrogen obtained from the force fields is crucial. In the Supporting Information (Section S5), we provide an extended discussion on the importance of the chemical potentials (focusing on water) regarding the accuracy of solubility calculations for water–hydrogen systems. Additional VLE calculations of water–hydrogen mixtures were performed using the TIP3P and TIP4P/2005 force fields. For the purpose of this discussion, a modified force field based on the TIP4P/2005 force field is introduced in the Supporting Information, which improves the accuracy of solubility calculations simultaneously in the liquid and the gas phase.
Figure 9 VLE of water–hydrogen at low temperatures as a function of pressure, obtained from eqs 10 and 11. (a) Mole fraction of water (y
H2O) in the gas phase. (b) Mole fraction of hydrogen (xH2~~) in the liquid phase. For each temperature (Tm), the points on the dashed line correspond to the composition calculated at the melting pressure (Pm) of ice; see Table 6. For each value of Tm, solubility calculations are performed for pressures in the range between Pmand 1000 bar. The melting temperature and pressure were selected based on the values from the IAPWS equation of state.^91,92^ Raw data, including the melting temperatures and pressures, are provided in Tables 8 and 9. Error bars are smaller than symbols.
Based
on the solubility of hydrogen in liquid water along the melting
line, one can calculate the freezing-point depression of water using eq 17. The results are shown
in Table 10. The largest
change in the freezing point of ice is around ΔTF = 1 K, corresponding to 1.1% dissolved hydrogen in water
at P = 1000 bar. At lower pressures, the freezing-point
depression is less than 1 K.
In this work, we investigated the thermodynamic properties of hydrogen–water mixtures using molecular simulation. Using molecular-based modeling, we obtained important thermochemical properties of hydrogen–water mixtures, which are not readily available from experiments. Ensemble fluctuations were used to calculate thermodynamic derivatives from molecular simulation. The chemical potentials of water and ice were calculated along the melting line of ice. Different force fields for hydrogen (Cracknell,^97^ Buch,^106^ Hirschfelder,^107^ Vrabec,^102^ and Marx^99^) were used to calculate thermal expansivities, heat capacities, and the Joule–Thomson coefficients at high pressures. The heat capacities and thermal expansivities calculated using the Vrabec force field were in good agreement with experimental data. The results from other force fields agreed qualitatively well with the empirical results (e.g., following similar trends). No force field could very accurately calculate the Joule–Thomson coefficient of hydrogen; however, the qualitative behavior of the thermodynamic properties of hydrogen is correctly captured. The results obtained from the Vrabec and Marx force fields are closest to the experimental data. The TIP3P force field for water was used for simulating hydrogen–water mixtures and performing the VLE calculations, at temperatures and pressures corresponding to the melting line of ice Ih. From the commonly used rigid water models with point charges, the TIP3P force field is one of the best force fields for calculating the solubility calculations without using any binary interaction parameters. To examine the effect of the chemical potential on the VLE of water–hydrogen mixtures, a modified force field based on the TIP4P/2005 force field was introduced. Compared to the TIP4P/2005, the accuracy of the computed chemical potentials was significantly improved for the modified force field, which consequently led to some loss of accuracy for density predictions in the liquid phase. The computed solubilities of water in the gas phase using the modified TIP4P force field were similar to those obtained using the TIP3P force field. However, the modified TIP4P force field outperformed the TIP3P force field for the solubility of hydrogen in liquid water. Simulation results showed that thermodynamic properties of compressed hydrogen significantly change as the mole fraction of water in the gas phase is increased. For example, in sharp contrast to pure hydrogen, the Joule–Thomson coefficient of a saturated hydrogen–water mixture, containing 0.7% water (at T = 423.15 K and P = 100 bar), has a positive value, clearly showing that the effect of water cannot be neglected. At pressures above 700 bar, the solubility of water in hydrogen is low so that the behavior of the hydrogen–water mixture becomes very similar to the behavior of pure hydrogen. Experimental validation of the effect of water on thermodynamic properties of hydrogen would be very helpful and complementary to the simulations results. Unfortunately, to the best of our knowledge, such experimental data of the exact system at similar conditions are not available in the literature. The chemical potentials of water and ice along the melting of ice were obtained from molecular simulation and the IAPWS equation of state. For each melting temperature, the chemical potentials of water and ice were calculated as a function of pressure using experimental volumetric data. From the VLE of water–hydrogen, the compositions in the gas and liquid phase were obtained at low temperatures. The solubility calculations were performed in the pressure range between 100 and 1000 bar. The solubility of water in the gas phase is the highest at 100 bar (272.4 K), which is around 90 ppm (molar), and the solubility of hydrogen in the liquid phase is the highest at 1000 bar (264.21 K), which is around 1%. Based on the solubility data in the liquid phase, the freezing-point depression of water was calculated. The largest freezing-point depression is around 1.1 K corresponding to a pressure of 1000 bar and temperature of 264.21 K. For lower pressures, the freezing-point depression is very small. Our simulation results may have consequences for the energetics of a hydrogen refueling station using EHCs. This includes drying or precooling of the compressed hydrogen stream.
This work was sponsored by NWO Exacte Wetenschappen (Physical Sciences) for the use of supercomputer facilities, with financial support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (The Netherlands Organization for Scientific Research, NWO). T.J.H.V. acknowledges NWO-CW for a VICI grant.
The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jced.1c00020.
The authors declare the following competing financial interest(s): One of the authors (Julio C. Garcia-Navarro) was an employee of HyET Hydrogen BV, a Dutch company that is developing Electrochemical Hydrogen Compressors.