CET vol 99 DOI: 10.3303/CET2399073 Paper Received: 5 January 2023; Revised: 29 March 2023; Accepted: 17 April 2023 Please cite this article as: Restelli F., Spatolisano E., Pellegrini L.A., 2023, Hydrogen Liquefaction: a Systematic Approach to Its Thermodynamic Modeling, Chemical Engineering Transactions, 99, 433-438 DOI:10.3303/CET2399073 CHEMICAL ENGINEERING TRANSACTIONS VOL. 99, 2023 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Flavio Manenti Copyright © 2023, AIDIC Servizi S.r.l. ISBN 978-88-95608-98-3; ISSN 2283-9216 Hydrogen Liquefaction: a Systematic Approach to its Thermodynamic Modeling Federica Restelli, Elvira Spatolisano, Laura A. Pellegrini* GASP - Group on Advanced Separation Processes & GAS Processing, Dipartimento di Chimica, Materiali e Ingegneria Chimica “G. Natta”, Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milano, Italy laura.pellegrini@polimi.it In the present work, a thermodynamic approach capable of describing the hydrogen behavior during its cooling and liquefaction is proposed both for the case of catalytic ortho to para conversion occurring inside dedicated reactors and for the case of continuous conversion inside heat exchangers where the catalyst is packed on the hydrogen side. The state-of-the-art Equation of State to describe the properties of normal-, para- and ortho- hydrogen is the Helmholtz free energy explicit equation. However, it can only describe pure components and not mixtures. The novelty of the proposed approach is that it is based on the widespread Peng Robinson Equation of State and that it allows to accurately describe the calorimetric and volumetric properties of the different forms of hydrogen and their mixtures. Furthermore, it can be easily implemented in the Aspen Plus® process simulator, resulting to be useful in view of design and optimization of the hydrogen liquefaction process. 1. Introduction Hydrogen (H2) is getting attention as energy carrier since its utilization does not produce greenhouse gas emissions and it can be used to store the renewable energy. The most cost-competitive renewable electricity hubs, for example the wind farms of the North Sea or the solar parks of the Middle East, are far away from the large productive sites in Europe. Therefore, supply chain infrastructures need to be developed to move H2 through long distances. The main obstacle to gaseous H2 delivery is its low volumetric density, which would lead to excessive dimensions of the storage and transport equipment. To increase the volumetric density, liquefaction is a promising option for large-scale delivery. H2 liquefaction occurs at about 20 K and atmospheric pressure. It is the most energy-intensive process of the supply chain and, hence, special attention must be paid in its design. H2 is a quantum gas and it exists in two isomeric forms, named ortho-hydrogen and para-hydrogen, whose relative composition at equilibrium changes with temperature. The isomerization reaction is very slow and must be accelerated using a catalyst in order to achieve the equilibrium ortho-para composition before storage, otherwise the exothermicity of the reaction will lead to the evaporation of a great portion of the stored liquid. Iron-based catalysts are the most used to this purpose. This ortho-para conversion (OPC) can occur into dedicated catalytic reactors placed between two heat exchangers along the cooling train, or continuously by packing the catalyst on the hydrogen side of the heat exchangers. Existing plants, like the ones constructed by Linde in Germany and by Praxair in the United States, are all based on a liquid nitrogen (LN2)-precooled Claude cycle. The Linde liquefier in Ingolstadt, now decommissioned, included reactors for the OPC. The most recent plants, like the Linde liquefier in Leuna, use catalyst-packed heat exchangers (Al Ghafri et al., 2022). Various conceptual designs of liquefaction processes have been proposed in recent years and their performances are discussed by analyzing the results from commercial simulators. Valenti and Macchi (2008) proposed a liquefaction process in which refrigeration is provided via four helium recuperative Joule-Brayton cycles, with the assumption of continuous OPC. They adopted Aspen Plus® as process simulator, with the Benedict-Webb- Rubin-Starling equation of state (EoS). Krasae-in et al. (2014) simulated a mixed refrigerant (MR) precooled liquefaction process, in which the cooling at low temperature is provided via four hydrogen recuperative Joule-Brayton cycles. The OPC to equilibrium is achieved in five reaction stages. The Soave-Redlich-Kwong EoS was selected in the simulation, carried out using PRO/II together with a C language program. Sadaghiani 433 and Mehrpooya (2017) proposed a MR cascade liquefaction process in which OPC is performed in two reactors. They adopted Aspen HYSYS®, selecting the Peng-Robinson EoS. A MR precooled dual-pressure Claude cycle and a MR cascade cycle were simulated by Cardella et al. (2017), assuming continuous OPC. They carried out the simulations using UniSim® Design software and calculating the pure H2 properties with the Helmholtz energy explicit EoS implemented in REFPROP. What emerges from the literature scouting is the lack of a systematic method to approach the simulation of hydrogen liquefaction. This arises the need to define a method, to be easily implemented in commercial simulators, accurate in predicting the calorimetric and volumetric properties of H2 during liquefaction and versatile in being used for different processes at different operating conditions. 2. Thermodynamic model definition The H2 isomers are distinguished by the orientation of the nuclei spins: the form in which the spins are in the same direction (antisymmetric) is called ortho-hydrogen (o-H2), the one in which the spins are in the opposite direction (symmetric) is called para-hydrogen (p-H2), as shown in Figure 1a. Beside these two forms, normal-hydrogen (n-H2) and equilibrium-hydrogen (e-H2) are commonly defined, representing mixtures of o-H2 and p-H2. e-H2 corresponds to the equilibrium mixture of o-H2 and p-H2 whose composition depends on temperature T, as shown in Figure 1b. Above the ambient temperature, the equilibrium composition remains fixed at 75% o-H2 and 25% p-H2: this mixture is called n-H2. The equilibrium between the isomers is due to the lower energy state of p-H2 compared to o-H2. Therefore, the conversion from the ortho to the para form is an exothermic reaction (ΔHr < 0) and the reaction enthalpy ΔHr is a function of T, as shown in Figure 1b. a) b) Figure 1: a) Orientation of nuclei spins in ortho- and para-hydrogen; b) Percentage content of p-H2 and enthalpy of the isomerization reaction -ΔHr [kJ/kmol] at equilibrium as functions of temperature. In this section, a method to approach H2 thermodynamics during liquefaction is presented, which is able to represent the effective behavior of the o-H2 and p-H2 mixture at any temperature level along the process, considering that in the case of continuous OPC, obtained by packing the catalyst on the H2 side of the heat exchangers, their relative content can be assumed to follow the one reported in Figure 1b. Instead, in the case of OPC occurring inside reactors placed at given positions along the cooling train, the composition remains constant along the heat exchangers and change only at the catalyzed reaction stages, reaching the equilibrium composition at the stage temperature. To the knowledge of the authors, o-H2 is not present in any commercial simulator. Aspen Plus® (AspenTech, 2016) is selected as process simulator because n-H2 and p-H2 are built-in components and because it is easy to customize. To simulation purposes, to represent the mixture of o-H2 and p-H2, a mixture of n-H2 and p-H2 can be used, considering that n-H2 consists of 25% p-H2 and 75% o-H2 in order to adjust the composition in terms of n-H2 and p-H2. This method is useful to represent the case of stepwise conversion. In the case of continuous conversion, a new component representing the e-H2 can be created. Considering that n-H2 and p-H2 have similar volumetric properties, except in the cryogenic zone (McCarty et al., 1981), and that at cryogenic temperatures e-H2 is mainly composed of p-H2, it is possible to assume that its volumetric behavior is the same as that of p-H2. Therefore, to create the e-H2 component, the p-H2 component is chosen, whose expression for the isobaric specific heat is modified, as explained in Section 2.1. In Aspen Plus® V.11, enthalpy and entropy at a given temperature and pressure are calculated as the sum of three quantities: the contribution due to the formation of the species starting from the elements at standard conditions, the contribution involved in bringing the species from standard conditions to the system temperature at ideal gas conditions, ΔHIG and ΔSIG, and the contribution involved in bringing the species to system pressure and state, ΔHdep and ΔSdep. These latter quantities are called departure enthalpy and entropy, respectively, and the method of calculation varies depending on the thermodynamic model used to represent the vapor and liquid phases. The terms referring to the properties of the ideal gas are calculated starting from the ideal specific heat o-H2 p-H2 0 250 500 750 1000 1250 0 20 40 60 80 100 0 50 100 150 200 250 300 -Δ H r [k J /k m o l] p -H 2 c o n te n t [% ] T [K] 434 at constant pressure and, therefore, depend on the calorimetric properties of the species. The terms referring to the departure properties, instead, depend on the volumetric properties of the species. The calorimetric and volumetric properties of the defined components, n-H2, p-H2 and e-H2, are detailed below. 2.1 Calorimetric properties The calorimetric property of interest is the ideal gas specific heat at constant pressure cp 0, which can be calculated according to the Mayer relation (cp 0= R+cv 0). From the kinetic theory of gases, the specific heat at constant volume cv 0 is the sum of different contributions, corresponding to the degrees of freedom of the molecule, classified as: translational, rotational and vibrational. For the H2, the vibrational degrees of freedom are negligible up to 600 K (Valenti et al., 2012) and therefore these contributions will be considered equal to zero in the present discussion. The expression of the rotational specific heat c0 v,rot was firstly derived by Hund (1927) for the ortho and para isomers on the basis of the wave mechanics. Dennison (1927) used the results of Hund to compute the c0 v,rot of n-H2, considering the ideal mixture of ortho and para in proportion of 3:1. Leachman et al. (2009) proposed the functional form of Eq(1), which takes up the theoretical results of Hund and Dennison in terms of dependence on T, but with the coefficients uk and vk regressed against experimental data. The same functional form is adopted by Valenti et al. (2012), who regressed the coefficients uk and vk against ab initio data, to compute the c0 v,rot of e-H2. In the present work, the correlation DIPPR 127 (Eq(2)), already present in Aspen Plus®, is proposed as an alternative to Eq(1). The coefficients are found by minimizing the root mean square deviation of cp 0 calculated with DIPPR 127 with respect to the values obtained by Leachman et al. and by Valenti et al. for different T in the range 20 – 300 K and are reported in Table 1. 22N 0 k k k v,rot k k 1 v v v c R u exp exp 1 T T T − =          =     −                 (1) ( ) ( ) ( ) 22 k 1 i k 1 i k 1 i0 p 1i ki k 2,4,6 C C C c C C exp exp 1 T T T − + + + =           = +    −                       (2) Table 1: Estimated coefficients for calculating the ideal specific heat at constant pressure cp 0 [kJ/(kmol·K)] of n-H2, p-H2 and e-H2 using the correlation DIPPR 127. C1 C2 C3 C4 C5 C6 C7 n-H2 20.7862 11.3784 513.39 -3.24626 1242.75 6.0525 4812.55 p-H2 20.7862 75.1267 561.48 66.4278 1037.49 -131.678 810.13 e-H2 20.7862 673.413 211 -797.49 230.13 133.449 326.14 Figure 2 shows the trend of cp 0 for n-H2, p-H2 and e-H2, calculated with DIPPR 127, as a function of T. The points highlighted as circles and triangles refer to the literature data presented by McCarty et al. (1981) and Le Roy et al. (1990), respectively. Since the Average Absolute Deviation (AAD) is 1.00 %, 1.23 %, 0.33 % for n-H2, p-H2 and e-H2, respectively, it is possible to assert that the correlation predicts the data correctly. Figure 2: Ideal specific heat at constant pressure cp 0 [kJ/(kmol·K)] as a function of temperature: comparison between the values obtained with the correlation DIPPR 127 and the literature data in for n-H2, p-H2 and e-H2. 2.2 Volumetric properties The volumetric property of interest is the molar volume v, which can be evaluated through Equations of States or correlations. The most accurate EoS available in the literature to describe the properties of normal- and para- hydrogen are the Helmholtz free energy explicit equation, proposed by Leachman et al. (2009), and the modified 0 5 10 15 20 25 30 35 40 0 50 100 150 200 250 300 T [K] McCarty et al. (1981) - McCarty et al. (1981) - Le-Roy et al. (1990) - Le-Roy et al. (1990) - c p 0 [k J /( k m o l∙ K )] n-H2 p-H2 p-H2 e-H2 cp 0 e-H2 cp 0 p-H2 cp 0 n-H2 435 Benedict-Webb-Rubin EoS. They are both implemented in the Aspen Plus® REFPROP package. However, they can only describe pure components and not mixtures. As stated before, it is possible to assume for e-H2 the same volumetric behavior as p-H2. Since in the case of a liquefaction process, which involves OPC reactors, hydrogen is represented in Aspen Plus® as a mixture of n-H2 and p-H2, an EoS is required that is accurate in describing both hydrogen as a pure compound and as a mixture. The Peng-Robinson EoS is selected for its simplicity, accuracy and popularity in the process industry. Since the cubic EoS is not accurate in predicting the liquid molar volume, the following two approaches are considered for its calculation. • the Costald model (Eq(3)) (Hankinson and Thomson,1979), in which the compressed liquid molar volume is computed starting from the value at saturation, calculated with the constant parameters vCTD, ωCTD and the temperature-dependent parameters vR,0, vR,δ, and correcting it by the pressure P using the Tait correlation (Thomson et al.,1982). B, C are functions of T, ω and critical parameters. Psat is the saturated pressure at T. ( )   + = − −   +    CTD R,0 CTD R, sat B P v v v 1 v 1 Cln B P (3) • the Chueh-Prausnitz model (Eq(4)) (Chueh and Prausnitz,1969), in which the saturated liquid molar volume vsat, calculated using the correlation DIPPR 105 for pure components and the Rackett model for mixtures, is corrected by the pressure to compute the compressed liquid molar volume, adopting as correlation parameters the critical compressibility factor ZC, the critical pressure PC and N, which is a function of the acentric factor ω and T. ( ) −  − =  +    1/9 satsat C C P P v v 1 9Z N P (4) a) b) c) d) Figure 3: Molar volume v [m3/kmol] as a function of temperature: comparison between the values obtained with the Peng-Robinson EoS, the Costald correlation, the Chueh-Prausnitz correlation and the literature data, for n- H2 at constant pressure equal to: a) 1 bar, b) 20 bar; for p-H2 at constant pressure equal to: c) 1 bar, d) 20 bar. Figure 3 shows the molar volume v computed with the Peng-Robinson EoS, Costald and Chueh-Prausnitz correlation (these latter only for liquid phase) in the range 14 - 300 K and at two pressure levels equal to 1 and 20 bar, together with the data collected by McCarty et al. (1981). The AAD values of the two correlations with respect to the literature data are equal to 12.5% for the Costald model in the case n-H2 at 20 bar and less than 3% in all the cases for the Chueh-Prausnitz model. Since this latter presents lower AAD values than the other equation, it is selected in the definition of the model. 0 5 10 15 20 25 30 0 50 100 150 200 250 300 v [ m ³/ k m o l] T [K] 0.01 0.03 0.05 10 15 20 25 30 350.01 0.03 0.05 10 15 20 25 30 35 McCarty et al. (1981) - n-H2 Peng-Robinson - n-H2 Costald - n-H2 Chueh-Prausnitz - n-H2 n-H2 n-H2 n-H2 n-H2 0.00 0.25 0.50 0.75 1.00 1.25 1.50 0 50 100 150 200 250 300 v [ m ³/ k m o l] T [K] 0.01 0.03 0.05 10 15 20 25 30 350.01 0.03 0.05 10 15 20 25 30 35 McCarty et al. (1981) - n-H2 Peng-Robinson - n-H2 Costald - n-H2 Chueh-Prausnitz - n-H2 n-H2 n-H2 n-H2 n-H2 0 5 10 15 20 25 30 0 50 100 150 200 250 300 v [ m ³/ k m o l] T [K] 0.01 0.03 0.05 10 15 20 25 30 350.01 0.03 0.05 10 15 20 25 30 35 McCarty et al. (1981) - p-H2 Peng-Robinson - p-H2 Costald - p-H2 Chueh-Prausnitz - p-H2 p-H2 p-H2 p-H2 p-H2 0.00 0.25 0.50 0.75 1.00 1.25 1.50 0 50 100 150 200 250 300 v [ m ³/ k m o l] T [K] 0.01 0.03 0.05 10 15 20 25 30 350.01 0.03 0.05 10 15 20 25 30 35 McCarty et al. (1981) - p-H2 Peng-Robinson - p-H2 Costald - p-H2 Chueh-Prausnitz - p-H2 p-H2 p-H2 p-H2 p-H2 436 3. Application of H2 thermodynamic modeling to a H2 liquefaction cycle Once defined an approach to hydrogen thermodynamics during liquefaction, in this section the method is tested by performing the simulation of a LN2-precooled Claude cycle considering both the alternatives for the catalytic OPC (Figure 4). a) b) Figure 4: Process scheme of the proposed LN2-precooled liquefaction process, taken with minor modification by Crawford (1963), in the case of OPC occurring: a) inside two reactors or b) along the heat exchangers’ tubes. Considering the scheme in Figure 4a, the inlet stream 1 is n-H2 at 20 bar pressure and ambient temperature and is supposed to be free of moisture and condensable gases. This stream is mixed with the recycle streams 16 and 20 of high p-H2 content to form stream 2, which is cooled using service water and then compressed to 103 bar in an intercooled two-stage compressor C-1. The outlet stream 4 from C-1 is cooled in multi-pass heat exchangers HX-1 and HX-2, to 64.8 K and 45.8 K respectively, by warming the recycle gaseous H2 streams at low and intermediate pressure. The heat exchanger HX-1 include a pool of LN2 boiling at a pressure of 1.1 bar. Both heat exchangers are characterized by a minimum temperature approach of 2 K. The process stream is sent to the reactor R-1, equipped with a suitable catalyst to promote the OPC, at the outlet of which it is assumed to reach the equilibrium composition, obtained in the simulation with a mixture of the components p-H2 and n-H2, as defined in Section 2. The reactor is interested by a thermal duty qr, as defined in Eq(5), because of the exothermicity of the ortho-para isomerization reaction: ( )( ) ( ) 2 2r r p H ,in p H ,outq H T x x − − = −  − (5) where xp-H2,in and xp-H2,out are the molar fraction of p-H2 at the inlet and at the equilibrium conditions at the temperature of the stage. The stream 6a at the outlet of R-1 is cooled in HX-2 before reducing its pressure to 5.5 bar in Joule-Thompson valve VLV-1. The resulting biphasic stream is separated in flash vessel V-1. The vapor is recycled back to cool the process stream in HX-1 and HX-2. The liquid is passed through the process- process heat exchanger HX-3, which is characterized by a minimum temperature approach of 1 K. The stream 9 is sent to the reactor R-2, equipped with a suitable catalyst for OPC, and then is cooled passing through the Joule-Thompson valve VLV-2, reaching the storage pressure of 1.3 bar. The flash vessel V-2 separates the vapor, which is recycled back passing as cold stream in HX-3, HX-2 and HX-1, and the liquid, which leaves the process to be stored in an insulated tank at 1.3 bar pressure. For sake of simplicity, pressure drops inside the heat exchangers are neglected. If the same process is carried out by packing the heat exchangers tubes with the catalyst (Figure 4b), it is possible to assume that the equilibrium composition is obtained throughout the cooling process. To represent this situation, the processed hydrogen stream is pure e-H2, as defined in Section 2. For comparison purposes the same pressure levels are considered, as well as the same minimum temperature approach inside the heat exchangers. The two described processes are compared in terms of specific electricity consumption (SEC) of the compressors and exergy efficiency (ηexe) of the heat exchangers. Table 2: SEC [kWh/kg] of each compressor and of the entire process, and ηexe [%] of each heat exchanger in the case of stepwise (Figure 4a) and continuous (Figure 4b) OPC. SEC [kWh/kg] ηexe [%] Type of OPC C-1 C-2 C-3 TOT HX-1 HX-2 HX-3 Stepwise 8.72 1.47 5.13 15.31 82.35 70.59 84.69 Continuous 8.23 0.54 5.13 13.90 82.99 75.61 84.96 437 The results are reported in Table 2. The SEC of compressor C-2 in the case of continuous OPC is about one third with respect to the one in the case of stepwise conversion. This is explained by the fact that the absence of reactor R-2 leads to a lower temperature of the stream 9 (Figure 4) entering VLV-2 and, thus, to a lower vaporization ratio inside the flash V-2. The lower electricity consumption of C-2 is a consequence of the lower vapor flowrate, exiting from the top of V-2 and flowing inside the compressor. The total SEC of the two alternatives highlights the advantage, quantified in 9% reduction of electricity consumption, of inserting the catalyst into the tubes of the heat exchangers. In HX-2, the ηexe is lower in the case of stepwise conversion. In fact, the addition of heat, due to the presence of the reactor R-1, during the cooling process causes a greater difference in temperature between hot and cold streams at the cold side of the exchanger. 4. Conclusions This study presents a novel approach to the simulation of H2 liquefaction, which can be easily implemented in Aspen Plus® process simulator. The method accurately describes the calorimetric and volumetric properties of n-H2, p-H2 and e-H2. The developed approach is used to assess the performances of a liquefaction process based on a LN2-precooled Claude cycle, in the case of OPC occurring stepwise inside dedicated reactors placed between the heat exchangers and of continuous conversion obtained by packing the catalyst on the hydrogen side of the heat exchangers. The results show that a reduction in the specific electricity consumption of 9% and a higher exergy efficiency is achieved in the case of continuous conversion with respect to stepwise case, as expected considering that the reaction enthalpy is lower, in absolute value, at high temperatures. This work is useful in view of the development of a green H2 value chain involving transport and storage on a large-scale. In fact, since liquefaction is the most energy-intensive step of the value chain, a reliable thermodynamic framework is of paramount importance to design the optimum process that minimizes the costs. References Al Ghafri S., Munro S., Cardella U., Funke T., Notardonato W., Trusler J.P.M., Leachman J., Span R., Kamiya S., Pearce G., Swanger A., Rodriguez E.D., Bajada P., Jiao F., Peng K., Siahvashi A., Johns M.L., May E.F., 2022, Hydrogen liquefaction: a review of the fundamental physics, engineering practice and future opportunities, Energy & Environmental Science, 15, 2690–2713. AspenTech, 2016, Aspen Plus®, Burlington (MA), United States. Cardella U., Decker L., Sundberg J., Klein H., 2017, Process optimization for large-scale hydrogen liquefaction. International Journal of Hydrogen Energy, 42, 12339–12354. Chueh P.L., Prausnitz J.M., 1969, A generalized correlation for the compressibilities of normal liquids. AIChE Journal, 15, 471–472. Crawford D.B., 1963, Hydrogen liquefaction and conversion systems, Patent US3095274A, Allentown, Pa, USA. Dennison D.M., 1927, A note on the specific heat of the hydrogen molecule, Proceedings of the Royal Society of London, Series A, 115, 483–486. Hankinson R.W., Thomson G.H., 1979, A new correlation for saturated densities of liquids and their mixtures, AIChE Journal, 25, 653–663. Hund F., 1927, On the interpretation of the molecular spectra (in Deutsch). II. Zeitschrift für Physik, 42, 93–120. Krasae-In S., 2014, Optimal operation of a large-scale liquid hydrogen plant utilizing mixed fluid refrigeration system, International Journal of Hydrogen Energy, 39, 7015–7029. Leachman J.W., Jacobsen R.T., Penoncello S.G., Lemmon E.W., 2009, Fundamental equations of state for parahydrogen, normal hydrogen, and orthohydrogen, Journal of Physical and Chemical Reference Data, 38, 721–748. Le Roy R.J., Chapman S.G., McCourt F.R.W., 1990, Accurate thermodynamic properties of the six isotopomers of diatomic hydrogen, Journal of Physical Chemistry, 94, 923–929. McCarty R.D., Hord J., Roder H.M., 1981, Selected properties of hydrogen (Engineering Design Data), US Department of Commerce, National Bureau of Standards, Boulder, Co, USA. Sadaghiani M.S., Mehrpooya M., 2017, Introducing and energy analysis of a novel cryogenic hydrogen liquefaction process configuration. International Journal of Hydrogen Energy, 42, 6033–6050. Thomson G.H., Brobst K.R., Hankinson R.W., 1982, An improved correlation for densities of compressed liquids and liquid mixtures, AIChE Journal, 28, 671–676. Valenti G., Macchi E., 2008, Proposal of an innovative, high-efficiency, large-scale hydrogen liquefier. International Journal of Hydrogen Energy, 33, 3116–3121. Valenti G., Macchi E., Brioschi S., 2012, The influence of the thermodynamic model of equilibrium-hydrogen on the simulation of its liquefaction, International Journal of Hydrogen Energy, 37, 10779–10788. 438 444restelli.pdf Hydrogen Liquefaction: a Systematic Approach to its Thermodynamic Modeling