Modeling and Simulation of the Hydrogen Value Chain with ThermoSysPro and Modelica

Hydrogen will play a key role in the global energy transition if we are able to produce it with low-carbon emissions. However, clean hydrogen production today is mostly limited to demonstration projects ranging from 2 up to 20 MW. Modeling the way low-carbon hydrogen is produced, stored and used will allow us to significantly improve our understanding of how clean hydrogen could be produced and thus increase the production efficiency. In this paper, we show that the newest version of Ther-moSysPro (TSP), an open-source Modelica library for modeling energy systems, provides a suitable framework to model and simulate the hydrogen production, storage and consumption. The model presented in this paper contains three electrolyzers, a storage station and a vehicle station. We present how the model was built, which components were adapted and how, and show that its simulation can be useful in the design phase, as well as for diagnosis purposes. In the future, a complete validation of these developments will be performed when experimental data is publicly available.


Introduction
The global energy system is undergoing major changes, (International Energy Agency 2022b).This transformation aims mostly at reducing the oil and other fossil fuels dependency in order to cut as much as possible CO 2 emissions and thus the consequences of global warming.In this context, hydrogen plays a crucial role in accelerating this energy transition by allowing the decarbonisation of key sectors like industry, aviation, shipping or heavy-duty transport, (International Energy Agency 2022a; Green hydrogen cost reduction 2020).
Hydrogen has already been identified as indispensable in the strategy to build a climate-neutral Europe (European Commission 2020).However, it is still mostly produced with natural gas, which represents 96% of the total production (European Commission 2023).For this reason, and in order to accomplish the ambitious climate objectives, it is extremely important to develop technologies able of producing clean hydrogen (i.e. with no or low CO 2 emissions) at a large scale and in a competitive way.
To successfully develop these hydrogen related technologies and integrate them in the future energy mix, the development of flexible and adequate modeling tools is also fundamental.This will allow to contribute and accelerate the ongoing energy transition.Hydrogen based systems can be modeled with different types of modeling tools.Without attempting to be exhaustive, in the literature it is possible to find Matlab (Khan and Iqbal 2005;Eriksson and Gray 2017;Bhuyan, Hota, and Panda 2018), Python (Kuroki et al. 2021) and Modelica (Migoni et al. 2016) based models.These models aim at representing different aspects of a hydrogen based system (electrolyzers, vehicle hydrogen fueling, storage, etc.) and for different purposes (optimization of the management of hydrogen based technologies, better understanding of the physical behavior of such systems, etc.).
Hydrogen based systems usually require different types of physics to be described and this is why the Modelica language can be suitable to describe such systems.In particular, it can be interesting to adapt existing Modelica libraries that have already shown relevant and trustful results in similar domains.This is the case of the ThermoSysPro (TSP) library1 (El Hefni and Bouskela 2019), an open-source Modelica library developed by EDF R&D2 , which has already been used to successfully model other energy systems such as power plants, industrial processes, energy conversion systems etc. (El Hefni and Bouskela 2017).To achieve these results, different physics are already available in TSP such as thermalhydraulics, combustion, solar radiation and neutronics.In addition, to these already existing features, the latest version of ThermoSysPro3 , the so-called V4, allows to easily modify the fluid considered for modeling in addition of giving the possibility to handle efficiently zero-mass flowrates configurations (based on the principles explained in (Bouskela and El Hefni 2014b).For the reasons previously mentioned, the ThermoSysPro library appears as suitable starting point to model and simulate hydrogen based systems.
The goal of this paper is therefore to expose and show how the latest version of the ThermoSysPro library has been adapted to efficiently model and simulate hydrogen based systems that already exist or that are going to exist in the future.To illustrate the modification of Ther-moSysPro previously described, this paper focus on the modeling and simulation of a low-carbon hydrogen platform in which the whole value chain of hydrogen is considered (production, storage and distribution, and consumption/application (TÜV SÜD 2023)).
This paper is structured as follows.Section 2 describes the model considered in this study as well as the equations implemented in the new developed modules.Section 3 presents the construction of the model as well as the simulations performed with this model.Finally, section 4 presents the discussion as well as the future work that can be initiated with the promising results presented in this paper.

Model Description and Equations
In this section, the proposed low-carbon hydrogen platform is presented as well as the components that constitute a library dedicated to hydrogen which is compatible with ThermoSysPro.As mentioned previously, TSP is already available and thus provides certain advantages like the possibility of reusing some of its components in this hydrogen platform.

Gas Plate Scheme
The proposed platform is denominated Gas Plate and it serves to connect the components of the system.Such components include mainly, as seen in the simplified model scheme on Figure 1, electrolyzers and a set of blocks as valves, compressors and heat exchangers to compose two branches: Storage and Station, each representing one of the considered final applications.First of all, the block corresponding to the electrolyzer allows to estimate hydrogen production through water electrolysis at given conditions.Such hydrogen gas will then go through a check valve, ensuring one-directional flow in the platform.Afterwards a compressor and a heat exchanger block will allow to reach the desired conditions in terms of temperature and pressure depending on each final use: either storing it directly or manipulating it in a filling station.Finally, to ensure the correct flow of hydrogen, a control valve will be used for each one of these so-called branches.

Notations
Ideal gas isochoric heat capacity ∆ r g 0

J mol
Reference specific Gibbs free energy difference for the reaction ∆ r h

Pure Real Gas Properties
The hydrogen platform operates at a pressure ranging from a few decades to several hundred bars.At high pressure, the ideal gas assumption is not precise enough and the use of an Equation of State for real gas is preferable.
The so-called Peng-Robinson Equation of State (Equation 1)4 can be found in the classical literature of thermodynamics and process engineering (Zohuri 2018): (1) where "a" accounts for inter-molecular attractive forces (it decreases the pressure that molecules exert on the reservoir (Zohuri 2018)), while "b" represents the volume occupied by the gas molecules.The molar volume is represented by v, P sat is the saturation pressure of the fluid (to find the acentric factor ω we use the value of P sat considering T r = T T c = 0.7) and R is the universal gas constant.Taking the total derivative of the internal energy (u = f (T, v)) and after development using Maxwell's third relation, we obtain the following expression verified by (Cengel and Boles 2015) where the pressure P is defined by the PR-EOS and the isochoric heat capacity C V is unknown.
Concerning Modelica, a function called PR_InternalEnergy (containing the Equation 4) is created after having developed the Equation 2 considering the state 1 as v 0 = ∞ and using the PR-EOS Equation 1 (Trujillo, O'Rourke, and Torres n.d.): Where u 0 = C 0 V T , and C 0 V corresponds to the isochoric heat capacity of perfect gases.Then, according to (Leachman et al. 2009) one can express such isochoric heat capacity according to the following expression: The parameters u k and v k are taken from (Leachman et al. 2009) for normal hydrogen (n-H 2 ) because in the temperature range of interest (≈ 260 − 360K) one will always have n-H 2 .Moreover, n-H 2 is defined as hydrogen's equilibrium concentration at room temperature, which corresponds to 75% orthohydrogen and 25% parahydrogen.
Finally, as there is interest in estimating enthalpy and entropy at different states, we take advantage of previous developments to use the following expressions in order to calculate such properties: However, Equation 8 is expressed more specifically as follows, where the integrals are to be solved: The derived formulas for internal energy, enthalpy and entropy can be found in the Appendix subsection 4.1.

Electrolyzer
To model the electrolyzer we start from the classical expression (Kuroki et al. 2021;Ulleberg 1998) 5 which includes the over-potentials in the electrodes of the cell: The term U rev concerns the reversible tension expressed by the Nernst equation (Atkins and Paula 2006).It is obtained by considering that two opposite forces of equal magnitude act on the ion at equilibrium, electrical forces and diffusion forces are equal in magnitude but of opposite sign.From such expression, while making certain assumptions6 , the Equation 11 is obtained to approximate the value of the reversible voltage.
where the vapor pressure P vH 2 O is defined by with pressures expressed in bar and T in Kelvin.The approximation ( 12) is valid between 25-250°C (Patterson et al. 2019).
To determine the last two terms of Equation 10(namely the over-potentials U res and U act , referring to the electrical and the activation resistance of the electrodes, respectively) we consider semi-empirical models found in the literature for alkaline electrolyzers.For the time being, two semi-empirical models were coded and an ulterior and more detailed work will compare both models as it is beyond the scope of the present work.Therefore, only the one found in literature will be presented here.The model in question is taken from the work of (Ulleberg 1998;Ulleberg 2003) which models a stand-alone photovoltaichydrogen power plant (PHOEBUS) in Jülich (Germany) equipped with an alkaline electrolyzer with circular bipolar cells and a 30 wt% KOH solution as electrolyte.
Where s i , t i , and r i are known parameters 7 , A cell is the area of a cell in m 2 , I is the current in amperes and the temperature T is expressed in °C.We can observe that Equation 14is an expression derived from the Butler-Volmer relation.
Then, for all that concerns the Faraday efficiency, we consider the model developed by (Ulleberg 1998) and already used in (Khan and Iqbal 2005) where the temperature T is expressed in °C: This calculation allows to estimate the real production of hydrogen in the electrolyzer considering the current losses and the number of cells N cell in the stack.
In addition, we want to estimate the heat losses since they could potentially be used as an energy source for the heating of the water inlet.This is obtained by an energy balance: Finally, we can also estimate the energy efficiency of our electrolyzer by using the lower heating value of hydrogen.
In Figure 2, we observe that there are inputs as working temperature (T re f ), the supplied current (I) and the water inlet at a given pressure (P re f ); as well as the outputs of oxygen and hydrogen flows on the top side of the electrolyzer model. 7The empirical parameters (s i , t i , and r i ) can be found numerically using non-linear regression techniques.Consequently one can fit the model for an alkaline electrolyzer in particular After having modeled the electrolyzer, it is pertinent to obtain the polarization curve (which represents the voltage as a function of current density) to verify that the cell voltage calculated using Ulleberg (1998) model has been well coded.From its shape and values, Figure 3 corresponds to the expected behaviour for alkaline electrolysis at two working temperatures (80 • C in red line, and 25 • C in blue line).The cell voltage U cell is represented by the continuous lines and the reversible voltage U rev is represented by the dotted lines.By observing the continuous lines (cell voltage), one could also specifically verify that when placed around low current densities, the activation over-potential (Equation 14) has a more important impact on the total voltage of the cell and at higher densities the behavior is rather linear due to the shape of the resistance over-potential (Equation 13).One can also verify that the reversible potential (Equation 11) is independent of the current density.

Compressor
The compressor is modeled using a isentropic efficiency coefficient.From Equation 8 isentropic processes may then be introduced.They refer to the assumption of equal entropy between two different states (1 − → 2).Namely, there is not any change when moving from state 1 to state 2, and then the left side of Equation 8 becomes equal to zero.Knowing the properties in such a state allows us to use the isentropic efficiency of the compressor defined by Equation 19.A graphical representation of the compressor model is given in Figure 4a.
Starting from the fact that the state 1 and the value of η is,comp is known, and considering ∆s = 0 in Equation 8, h 2is , the enthalpy in the state 2 is where s 2is = s 1 is calculated, and then h 2 for state 2 is obtained from Equation 19.
Moreover, when looking at the safety aspect, one must not forget that there are limits concerning the temperature and pressure that the gas should reach.For example, the SAE J2601 standard mentioned in the section "SAE2020" indicates these upper limits during hydrogen refueling.So for compression a maximum temperature of 135°C (Institute 2007) may be established (recommendation for high pressure hydrogen rich services) when using a piston compressor.This safety measure consequently sets a maximum compression ratio and therefore probably the need to compress in several stages.

Heat Exchanger
For the heat exchanger used to cool the hydrogen, the concept of the Logarithmic Mean Temperature Difference (LMTD) should be introduced first.According to Newton's law of cooling, the rate of heat transfer is related to the instantaneous temperature difference between the hot and cold media (in a heat transfer process, the temperature difference varies with position and time).As a result, the temperature variation is non-linear and can best be represented by a logarithmic calculation, hence the LMTD (Engineering ToolBox 2003).
Where ∆T a is the temperature difference between the two flows at the A end, and ∆T b is the temperature difference between the two flows at the B end (sides of the exchanger whether if its design is that of a counter-current heat exchanger or not), Q [W ] is the heat service exchanged, U [ W m 2 K ] is the overall heat transfer coefficient, and is the area of exchange.After developing Equation 20we obtain an expression according to the specific heat capacities of "cold" and "hot" fluids (C p, f and C p,c respectively): As the main interest of this first work is not particularly the heat transfer phenomenon, the heat exchanger model is quite simple which entails to having certain limitations that must be indicated to not make physical mistakes while using the component.Such limitations and assumptions will be discussed in Section 4

Pressure Control Valve
The valve flow coefficient (C v ) is defined as the flow capacity of a control valve under fully open conditions relative to the pressure drop across the valve.It is defined as the volume of water (GPM in the U.S.) at 60°F that will flow through a fully open valve with a pressure differential of 1 psi across the valve.It is useful to know how to calculate C v because it is the standard method of sizing and selecting control valves used in the industry.In addition, an oversized valve can lead to control problems such as flushing or poor heat transfer (∆T ) through a coil due to overflow.Conversely, an undersized valve may not provide sufficient flow and exceed the available ∆P.In general, the volume flow rate q and C v are related by an expression of type Where G is the specific gravity of the fluid (G = ρ/ρ air ).However, the functions to describe the behavior of the gas passing through the valve are generally more complex.
According to (Swagelok 2007) they can be expressed with the following two equations, where p = [bar] , T = [K] , q = std L min : N 1 and N 2 are parameters provided in (Swagelok 2007) that depend on the units used and the volumetric flow rate q is expressed under standard conditions (25°C and 1 atm).On the other hand, it is also necessary to illustrate this valve from a thermodynamic point of view to understand the physical process.Choke valves are generally small devices, and the flow through them can be assumed to be adiabatic (Q ∼ = 0) since there is not enough time or area for efficient heat transfer to occur.There is also no work done (w ∼ = 0), and the change in potential energy, if there is any, is very small (e p ∼ = 0).Even though the output velocity is often considerably larger than the input velocity, in many cases the increase in kinetic energy is negligible (Cengel and Boles 2015).Then the energy conservation equation for this single-stream steady-state flow device reduces to: Thus, the end result of a throttling process depends on how much a certain property has increased during the process.If the flow of energy increases during the process (P 2 v 2 > P 1 v 1 ), this can be at the expense of the internal energy.As a result, the internal energy decreases, which is usually accompanied by a decrease in temperature.If Pv of the product decreases, the internal energy and temperature of a fluid will increase during a throttling process 8 .In these cases, the magnitude of the ∆T is governed by a property called the Joule-Thomson coefficient : The Joule-Thomson coefficient is a measure of the variation of temperature versus pressure during a process with constant enthalpy, which corresponds to the situation of the valve.So if this coefficient is negative, the temperature increases during an expansion.The sign of the coefficient will therefore depend on the conditions in which the gas is, and the so-called inversion temperature (Figure 6) at these conditions.Thus, considering such phenomenon, a valve model is adapted (Figure 4b) from TSP to take into account the effects on temperature of a throttling process with hydrogen as a fluid.From such development, and through a comparative simulation, one can observe the difference in properties at the end of the expansion process using the TSP valve 9 and the hydrogen-adapted Pressure Control Valve (PCV).As a simple study case, a simulation is carried out with the following parameters for both of the valves: • Inlet pressure: 400bar • Hydrogen flow rate: 0.0038885674 kg s • Inlet temperature: 65°C 8 In the case of an ideal gas, h = h(T ); therefore, the temperature must remain constant during a throttling process, which is not the case for a real gas where h = h(T, P) 9 The valve available on TSP was initially designed for water and steam flows applications In this comparative simulation, the valve opening is controlled to obtain a certain outlet pressure that will vary from 150 to 350 bar.The goal is to analyze the change of the outlet temperature as a function of pressure and while maintaining the same inlet pressure, thus generating a different pressure drop value at each moment.
From the main result (Figure 7), it is observed that the outlet temperature of the PCV (PCV _realsgas.T ) evolves in function of the outlet pressure (sensorP.C2.P), which is not the case when using the TSP default valve (PCV _idealgas.T ) that considers an ideal fluid.Specifically, it is noted that for the operating conditions of temperature and pressure, the Joule-Thomson coefficient is always negative, causing a temperature rise during the isenthalpic expansion process, which can be verified in Figure 6.Therefore, it is proved that the modified valve does include the Joule-Thomson phenomenon through the use of PR-EOS which describes a real gas.Even though the results have not been compared to real values from experimentation, the PR-EOS is widely accepted as an approximation in the process engineering field, however the use of a hydrogen-dedicated equation of state could be envisaged (Sakoda et al. 2012).

Gaseous Storage
The modeling of the pressurized gas storage is summarized in (Migoni et al. 2016) where an equation of state is implemented to estimate the pressure.However, unlike such work which uses the ideal gas equation, one can improve the estimates by using the PR-EOS.For that the Equation 1 which allows to model the pressure is used.Then it is necessary to know the mass stored in a given moment.For that we use a mass balance (Equation 26) which allows to calculate the quantity of hydrogen in each bottle and finally to estimate the pressure according to this mass and the volume of the tank.The number of bottles in the storage is a parameter that can be modified and the way in which the bottles are filled is directly managed by the model itself (a maximum and a minimum pressure in each bottle can be provided in order to manage the switches).
3 Gas Plate modeling

Model construction
The case study model (Figure 1) was briefly presented in subsection 2.1 where the general idea of the so-called Gas Plate was introduced alongside the two branches (Storage and Station) which illustrate the final usages of the lowcarbon hydrogen produced upstream through water electrolysis.At this time the theoretical fundamentals of each of the components (or blocks) in the platform has been settled and so a more detailed version of the model may be presented.Such is the case for Figure 8, where three electrolyzers are placed upstream on the left, and then a series of components such as pressure loss blocks at the outlet of each electrolyzer, valves 10 , volumes, sensors (to read temperature, pressure and mass flow) compose the hydrogen pipeline towards two branches or usages.On the right upper side (Storage Branch) a pressure control valve, a compressor, and counter-current heat exchanger are used to regulate the flow, to obtain a desired pressure and to refrigerate the hydrogen flow to a safe temperature, respectively.This, to finally arrive to the gaseous storage system at the end of the branch.On the other hand, the charging station branch in our chain is composed of a compressor to increase the pressure to the desired values according to the type of vehicle, a valve (PCV) that allows to regulate the output pressure according to the normative (SAE 2020), and a heat exchanger between those previous components in order to limit the temperature rise due to the Joule-Thomson effect in the valve.
Finally, it must be noted that some of the components used in this platform are taken from the TSP library.Such 10 The Check Valve right before the volume that serves to divide the flow into the two branches is found on TSP (El Hefni and Bouskela 2019) and the PCVs used once on the Storage branch and twice on the Station branch is adapted from the Control Valve model of the same TSP library is the case for the Singular pressure loss component at the outlet of each electrolyzer, the volumes, the check valve and other common blocks as signals and sources.Components like the sensors, the compressor and the pressure control valve where adapted from TSP; other blocks like the electrolyzer and the storage tanks where developed for a new hydrogen-dedicated library compatible with TSP.

Model simulation
The model presented in Figure 8 and described in the previous subsection has been used to perform various dynamic simulations that have allowed to validate and better understand the physical behavior of the hydrogen platform.These simulations give as well an overview of the different types of predictions that can be performed with this model.
For this purpose, three different scenarios have been considered.These dynamic simulations in which several events occur are described below together with graphs showing the evolution of key variables for each scenario.
Scenario 1: the three electrolyzers are working during the whole simulation and at the beginning the hydrogen storage is being filled up.At 500s, the outlet mass flow rate of the storage is temporarily increased (see dotted red line in Figure 9).During this simulation it is possible to observe the filling of gaseous storage while there is no increase in the outlet mass flow rate, see continuous blue line in Figure 9 in which the evolution of the hydrogen mass in the first bottle (and only bottle in this simulation) of the gaseous storage is presented .When the outlet mass flow rate increases, it can be observed how the bottle is being emptied before resuming the filling at the end of the simulation.
Scenario 2: as in the previous scenario, the three electrolyzers are working during the whole simulation and at 1200 seconds the valve of the storage branch of the platform is partially closed (as shown by Figure 10a).The consequence of closing this valve is directly observed in the storage inlet mass flow-rate evolution, see Figure 10b where a remarkable decrease is shown after 1200s.In addition to this important decrease, it is also possible to notice in the previous graph a small perturbation of the inlet mass flow-rate at around 800s.This reduction can be explained by the graph plotted in Figure 10c and in which the evolution of the hydrogen mass in the storage is represented: the blue line corresponds to the total hydrogen mass of the first bottle and the red one to the total hydrogen mass in the second bottle.The above-mentioned perturbation of the inlet mass flow-rate occurs when the switch between the filling of the first and second bottle occurs (corresponding to the moment at which the maximum allowed pressure in the first bottle is reached).The valve closing occurs therefore when the second bottle is being filled and the consequences can be observed in the red curve where the increase of the hydrogen mass in the second bottle is reduced (slope decrease) after 1200s.
Scenario 3: contrary to the previous scenarios, in this

Discussion and Further Work
In this paper, we illustrated how to use the latest version of ThermoSysPro Modelica-based library, the so-called V-4, in order to model a hydrogen platform where hydrogen is produced by three electrolyzers, stored in a storage station and consumed by a vehicle station.It turns out that one advantage of using TSP is that some components of TSP library can be reused and easily adapted to hydrogen modeling.In our case study, we created two new components in TSP (the electrolyzer and the storage station) and adapted the other components (the compressor, the pressure control valve and the heat exchanger).
Even if the results obtained so far are already extremely encouraging and useful for the current studies, it is important to keep in mind that the modeled components have some limitations.For instance, the dynamics were neglected for the pipes, the heat exchanger and the electrolyzer.Also, the heat transfer equations might not be adapted for fast transients.Moreover, even though the last version of TSP can handle zero flow rate, this feature has not been tested in the present work and fast closing of valves may require solutions like multi-mode simulation for instance (Bouskela and El Hefni 2014a;Bouskela 2016).In addition, some simplifications have been made in the architecture of the gas plate with respect to a real installation.For example, here we only modeled the electrolyzer behaviour without considering the water purification or gas separation.We also merged the various compression and cooling stages into one single stage.
Considering the above mentioned limitations, we argue that the model presented in this paper is good enough for what it was meant for: model a hydrogen platform at a large-scale with Modelica and show its potential applications.The model presented in this work shows how this modeling approach can be used in the design phase of an hydrogen platform (determine an adequate architecture, study the dynamic response of the platform, etc.) and could be used as well in the operation phase of the platform.In this later case, the model developed would correspond to an existing platform, and the results provided by this model could be combined to on-site measurements in order to provide accurate diagnosis.This diagnosis could be performed using advanced mathematical methods that allow to combine simulated and measured and that have recently been adapted to Modelica models such as Data Proceedings of the Modelica Conference 2023 October 9-11, 2023, Aachen, Germany DOI 10.3384/ecp204229 Assimilation (Corona Mesa-Moles, Argaud, et al. 2019;Corona Mesa-Moles, Henningsson, et al. 2021) or Data Reconciliation (Bouskela, Jardin, et al. 2021).
In addition of the above-mentioned advanced uses of the models, in future work, the model developed could be used to optimize other systems of the hydrogen platform such as the associated instrumentation and control (using a linearized version of the model) or to perform technicaleconomic optimization.This kind of approach could be used as well to optimize the overall architecture of the hydrogen platform (number of electrolyzers, storage volume, number of branches for the fuelling station, etc.).In addition, this model can be easily enriched to include for example the source of energy used for the electrolyzers (for instance wind or solar energy sources).It is however important to keep in mind the possible simulation difficulties that may occur when increasing the size of the model or combining different kind of physics.Furthermore, the TSP developments to correctly model hydrogen based systems could be continued.For example, the library can be completed to model more detailed components and physical phenomena as needed or required to correctly model the behaviour of the real system.To reach this goal it is important to validate the models with experimental data, however, as far as we know, there is no publicly available data on real experiments to validate the developments presented in this article as a whole.For now, only validated models in the literature can be considered, such as the filling stations models (Kuroki et al. 2021) which have been validated with experimental data (SAE 2020).

Figure 2 .
Figure 2. Model of Alkaline Electrolyzer with corresponding inputs and outputs

Figure 4 .
Figure 4. Example of model components.

Figure 5 .
Figure 5. Test model of the heat exchanger component

Figure 6 .
Figure 6.Inversion temperatures for three real gases: nitrogen, hydrogen and helium (Central Michigan University 2013)

Figure 7 .
Figure 7. Outlet temperature comparison for two valve models at variable outlet pressure