\pi\_{\rm CA} = \frac{3}{2} \frac{D\_{\rm CA}}{(\rho c)\_{\rm CA} \cdot c^2 \delta\_l} + 0.5 \tag{12}
$$

and

$$\tau\_{\rm NCA} = \frac{3}{2} \frac{D\_{\rm NCA}}{(\rho c)\_{\rm NCA} \cdot c^2 \delta\_{\rm I}} + 0.5\tag{13}$$

where *c* is a pseudo sound speed, *DCA* and *DNCA* are the diffusivity coefficient in crack area and non-cracked area, respectively. To ensure the continuities of the flux and density distribution functions at the interface of different phases, the ρ*c* is set to 1 in this model [16]. Then the calculation of macroscopic quantities of the density and flux can be obtained by

$$T = \sum\_{i} f\_{i},\tag{14}$$

$$J = \left(\sum\_{i} c\_{i} f\_{i}\right) \frac{\pi - 0.5}{\pi}.\tag{15}$$

At the end of each computational loop, the effective diffusivity coefficient *Deff* can be obtained by Fick's law [17, 18]:

$$D^{\rm eff} = \frac{L\_{\chi} \int J dA}{(C\_2 - C\_1) \int dA}.\tag{17}$$

where *C*2 and *C*1 are the boundary values of inlet and outlet boundaries, respectively. *Lx* is the distance between these two boundaries, and *A* is the computational width of the domain.

The normalized effective diffusivity coefficient *Dn* is defined as [19],

$$D^n = D^{\rm eff} / D\_{NCA} \,. \tag{18}$$

The*DNCA* is equal to the*DO*2 , which is the oxygen diffusivity in bulk. This normalized method makes the crack effect obviously on the diffusivity change of the CLs.

### **4 Results and Discussion**

#### **4.1 Model Validation**

Two typical structures, the "parallel" and the "series" model which are introduced in the previous work are used to validate the LBM model [20]. In the model validation section, the CR is set to 0.5. The *DCA*/*DNCA* value is set to 1–1000. Meanwhile, the domain side length effect is also taken into consideration which is set from 300 to 5000 µm at the resolution of 10 µm. The model validation result is shown in Fig. 2. When the domain side length is small, the error value is big, and the error is sensitive to the ratio between two conductivity values. To balance the computational consumption and an error less than 1%, the domain side length is set to 3000 µm, and the conductivity ratio is limited at 1:100. In another word, the representative element volume (REV) in this case is 3000 µm × 3000 µm.

**Fig. 2.** Comparison between model prediction value and theoretical value.

#### **4.2 IP Effective Transfer Properties Influenced by Cracks**

**Normalized Concentration with Diffusion Time**.  The porosity of NCA is 0.3, and five different crack quantities are selected from 20 to 180. The triangle crack with 25 µm crack width are reconstructed, and every crack arm length is random from 0 to 200 µm. As shown in Fig. 3, the average temperature at each coordinate with 180 cracks is higher than that of 20 cracks. The *DCA* is significant than the *DNCA*, thus a higher *CR* is good for mass transport in the CLs. After 15e5 lattice steps, the computing converges, but none of the concentration reaches the theoretical linear line. This is because the existence of cracks changes the real tortuosity of the mass

**Fig. 3.** The normalized concentration distribution with different time steps.

transfer. Furthermore, the more crack quantity, the farther the concentration deviates from the theoretical value.

**Normalized Concentration with Crack Shape**. To evaluate the effect of crack shape on the on normalized concentration in the crack area, three kinds of domains with different crack shape, which are linear cracks, circle cracks and triangle cracks were employed. The porosity of NCA is 0.3, and the *CR* were set to 0.15 ± 0.05. The triangle cracks case shows a little higher normalized concentration than that of linear crack and circle crack cases. However, when the calculation converges, all three cases shows similar normalized concentration distribution. In addition, the deviation of the computational value from the theoretical value is caused by the cracks affecting the tortuosity of the computational domain (Fig. 4).

**Crack Ratio Effect on Diffusivity**. To understand the effect of crack ratio on effective diffusivity, the relationship between crack ratio and normalized diffusivity was carried out. The NCA porosity is set to 0.2 and the *DNCA*/*Dbulk* is 0.04111. Different crack width, length and quantity cases are carried out. Meanwhile, the parallel model and series model are also used as the upper and lower limit of the normalized diffusivity. With the CR increase from 0 to 0.15, the normalized diffusivity of crack width case, crack length case and crack quantity case increase from 0.0411 to 0.0682, 0.0780 and 0.0866, respectively. Thus, the crack length shows the highest influence on the IP mass transfer enhancement. The computational model result is between the parallel model and series model results (Fig. 5).

**NCA Porosity Effect on Diffusivity**. A lower porosity will cause a lower *DNCA*, and the relationship between these two parameters is shown in Eqs. (2)–(6). When the porosities are 0.2, 0.3, 0.4 and 0.5, the *DNCA*/*Dbulk* are 0.0411, 0.089, 0.1579 and 0.2452, respectively. And the simulation is carried out for these four cases. Result shows, with the crack ratio increase from 0 to 0.17, the IP diffusivity enhancements of 0.2, 0.3, 0.4

**Fig. 4.** The normalized concentration distribution of different crack shapes.

**Fig. 5.** Diffusivity change caused by crack ratio change.

and 0.5 porosity cases are about 96.93%, 63.86%, 43.40% and 30.08%, respectively. Considering the relative low porosity of CL structure, this enhancement shows great potential to enhance the IP mass transportation of oxygen phase and water vapor phase (Fig. 6).

**Fig. 6.** Diffusivity change caused by NCA porosity change.

## **5 Conclusions**

In order to understand the crack effect on the PEMFC catalyst layer IP diffusivity, a lattice Boltzmann method D2Q9 model was employed in this research. The normalized concentration of mass transport was acquired, and the effective diffusivity of the complex crack zone was calculated by the flux result. The following conclusions were mainly obtained:

## **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Brief Review and Technical Insight of Liquefied Hydrogen Carriers Development**

Wei Song, Zilong Zhu(B) , Zhong Wan, Chan Long, and Tao Liu

Hudong-Zhonghua Shipbuilding (Group) Co., Ltd, Shanghai 200129, China<EMAIL_ADDRESS>
**Abstract.** Countries worldwide are shifting to a hydrogen economy in response to stringent environmental regulations, and hydrogen transport between countries is expected to increase in the medium to long term. Although hydrogen is traded between countries in different forms such as ammonia, liquid hydrogen, and methanol, from the perspective of volume density and production/demand area without a separate process, the transportation of hydrogen in liquid form is the potential way for large-scale transportation of hydrogen in the future. This article aims to highlight the opportunities and challenges technical for the ocean-going liquid hydrogen carriers. An overview of development state-of-theart and key technical challenges of liquid hydrogen carrier ships are summarized, including regulation, the cargo containment structure and insulation, boil off ratio (BOR) evaluation, boil-off gas (BOG) handling system and propulsion system. Finally, detailed technical route of the key technology required by future liquid hydrogen carrier is extrapolated, and securing a possible design through various technological alternatives.

**Keywords:** Liquid hydrogen carrier · Insulation · Boil off ratio · BOG handling · Propulsion system · Technical insight

## **1 Background and Motivation**

With the strengthening of international environmental regulations and the geopolitics of the Russia-Ukraine war, the energy market will undergo rapid changes in the future. The Paris Agreement adopted on December 12, 2015 is not only a substitute for emotion in Kyoto, but also a consensus of the international community on greenhouse gases (GHGs) [1]. Most advanced and developing countries have participated in the agreement, 186 countries put forward the goal and contribution plan of reducing greenhouse gas emissions. Investment in new renewable energy has being increased in recent years, but the fluctuations in new renewable energy production and difficulties in trading renewable energy resources across countries are considered the biggest hurdles obstacles to the utilization of renewable energy. To break through the limitations of this new renewable energy, hydrogen energy has attracted attention. Hydrogen uses fuel cells to generate carbon free electric energy and heat energy that can be easily converted into renewable energy [4]. Hydrogen energy produced by water electrolysis with renewable energy is best solution to balance the problem of renewable energy fluctuation and imbalance in various countries. China is today the largest hydrogen consumer in the world, at about 24 MtH2/year in 2020 [5].

Niermann et al. investigated hydrogen exports from Algeria to Hamburg [6] and analyzed the feasibility of different liquid organic hydrogen carriers (LOHC), pipeline transportation and bulk transportation of liquid hydrogen. In a similar study, chain energy efficiency and costs for ammonia and LH2 sea transport from northern Norway to Rotterdam and global (Tokyo) were estimated by Ishimoto et al. [7]. The energy efficiency and life-cycle costs of energy transported by submarine high-voltage cables are compared with pipeline hydrogen transport, compressed hydrogen ship and liquid hydrogen with different distances in [8]. Hydrogen can be transported by ship in the form of LOHC, ammonia or liquefied hydrogen and ammonia carriers have been commercialized and widely used in industrial. For LOHC carriers, ordinary chemical carriers can be used, so it is not necessary to further develop new technologies [9]. The above three recently published overlap papers, they all agree that liquid hydrogen (LH2) is the most promising option for long-distance seaborne hydrogen transport. Coincidentally, the liquefied hydrogen transport ship aims to the demonstration of liquefied hydrogen transport between Australia and Japan. At the end of 2019, Kawasaki Heavy Industry of Japan built a world's first liquid hydrogen carrier with two 1250 m3 double-shell vacuum liquid hydrogen storage tank, which horizontal cylindrical pressure vessel freely enable thermal shrinkage for transporting LH2 [10].

Several economic analyzes on the hydrogen supply chain including maritime transportation of liquefied hydrogen have been performed, but few studies has been done on the technical characteristics of liquefied hydrogen carriers. Although Japan has completed the construction of liquid hydrogen ships, this is only an experimental attempt. Therefore, this is a comprehensive and difficult task, because the volumes are vast and data sets, and the necessary methodology statement may not be fully available. A comprehensive comparison of theories and assumptions, methodological choices and levels of technical abstraction is beyond scope of this paper. Therefore, we review and discuss the most obvious differences in the existing liquid hydrogen carrier or conceptual design assumptions, and extrapolate the preliminary technical appearance of the future liquid hydrogen ship. While identifying detailed technologies necessary for the development of liquefied hydrogen carriers in the future and securing various technical alternatives, the technical feasibility of liquefied hydrogen carriers is analyzed.

### **2 Classification of Liquid Hydrogen Carrier Tanks**

The capacity of the liquefied hydrogen carrier depends on the economics of the liquefied hydrogen supply chain which have a significant impact on the analysis. Therefore, referring to the LNG ship type, the possible capacity of each tank type of the liquefied hydrogen carrier is analyzed. As shown in Fig. 1, membrane and type B tanks are applied to large LNG carriers, and Type C tanks are mainly applied to small LNG carriers or bunkering ships.

The Type C tank has the advantage of higher design pressure and ability to store BOG in the tank due to its relatively pressure build-up (accumulation) system locking

**Fig. 1.** Classification of cargo tanks for ships

BOG into the cargo tank, but has the disadvantage of being difficult to enlarge due to its shape and having low space efficiency. Hyundai Heavy Industries Group completed the conceptual design of a 20 K class liquid hydrogen carrier and received Approval in Principle (AIP) certification from the Korean Register of Shipping. Based on this, it is expected that it will be possible to manufacture up to 5000 ~ 6000 m3 per type C tank, and it is judged that the capacity of the liquefied hydrogen carrier can be up to about 20 K depending on the number of tanks.

Membrane type tanks have high space efficiency and can be manufactured with a capacity of 160 K or higher based on LNG carriers. However, since this is a judgment from the viewpoint of manufacturability of the tank, it is also necessary to consider the possibility of applying an appropriate insulation system according to the capacity.

## **3 State of the Art and Extrapolations of Hydrogen Tanks Boil off Ratio**

#### **3.1 Type-C Tank**

Liquid hydrogen storage and handling technology was first developed and applied by NASA aerospace projects [11]. The world's largest spherical LH2 storage tank with approximately 3200 m3 LH2 capacity was built in the 1960s. The absolute value of the boil off gas is reported about 530 gal/day, which corresponds to boil off ratio of approximately 0.0625% per day [12]. Different types of LH2 tanks have been designed and brought to market by Linde, the BOR of which is decided by size, shape, insulation, environment and usage pattern. For example, a cylindrical tank with a capacity of 300 m3 has a boil off ratio of 0.3% per day, while a spherical tank with a capacity of 1100– 2300 m3 has a boil off ratio of less than 0.1% per day [13].

In the 1980s KHI built spherical LH2 tanks with a volume of 600 m3 and 540 m3 LH2 capacity which achieves a boil off ratio of 0.18% per day [14]. After 30 years of operation, no degradation in insulation performance has been detected by KHI [15]. By 2020, 1250 cubic meters of horizontal cylindrical seaborne tank have been put in operation [16], the spherical onshore terminal tank achieves thermal insulation performance with ≤ 0.1% per day boil off ratio [17]. Recently, a basic principle design of 11,200 m3 volume spherical LH2 tank with ≤ 0.1% per day boil off ratio performance have been completed by KHI and certified by the classification society [18].

Most hydrogen storage vessels are double-layered with a vacuum in the middle. The space between these can also contain other materials, such as aluminum-coated polyester sheets, alternating layers of aluminum foil as well as fiberglass. The vacuum aims to reduce losses by conduction and convection, while the alternated layer aims to reduce losses by radiation [19]. Heat loss can also be reduced by reducing the ratio of the exposed surface to the volume of the tank, which is the reason why spherical tanks are often used to store the liquid hydrogen. That can be seen as a trade-off with cost. Although the surface-volume ratio of cylindrical tanks is higher than that in spherical tanks, they are easier to manufacture due to the low costs, thus making them more common. The boil-off ratio depends on the size of the tank and the intended pattern of use. For example, a 300 m3 small tank has boil off ratio of 0.3% per day, while tanks of 1100–2300 m3 can already achieve boil off ratio of less than 0.1% per day [20].

Based on the brief review above, it is reasonable to assume that low boil off ratio can also be achieved when LH2 cargo tanks are sized close to the capacity of today's Liquefied Natural Gas (LNG) carriers. Since boil off ratio of 0.06–0.25% per day is already achieved for tanks between 100 and 4000 m3, a natural question is that whether it will facilitate or impede achievable performance of low boil off ratio if further scalingup tanks volume. Under the assumption of well insulated tank with a uniform temperature distribution inside, the heat flow into tank LH2 will proportional to the difference between the ambient temperature and the LH2 temperature (K), overall heat transfer coefficien t(W m−2 K−1) and the tank surface area (m2), as shown in Eq. (1).

$$\mathcal{Q}\_{\rm in} = (T\_{\infty, ambient} - T\_{\infty, \rm LH\_2}) U\_{\rm overall} A\_{\rm surface} \tag{1}$$

The boil off ratio is defined as the ratio of the amount of evaporated gas produced per unit time to the full tank inventory, but [% per day] is a more commonly used unit for cryogenic storage tanks. Qin is heat (kW), *V*tank is volume (m3), ρLH2 is the density of liquid hydrogen (kg/m3), Levap,LH2 is the latent heat of vaporization (kJ/kg), *A*surface is the surface area (m2), and *U*overall is overall heat transfer coefficient (W m−2 K−1). The boil off ratio (% per day) can be expressed as:

$$BOR(\%/\text{day}) = \frac{\mathcal{Q}\text{in}}{L\_{\text{evap,LH}\_2} \cdot \rho\_{\text{LH}\_2} \cdot V\_{\text{tan k}}} \cdot 24 \cdot 3600 \cdot 100\% \tag{2}$$

where the units in the equation are international standard.

From (1) and (2), it can be found that the boil off ratio is proportional to the ratio of surface area to volume, which is usually called specific surface area. For spherical and cylindrical tanks, Fig. 2 gives how the specific surface area changes with size. It can be observed that it will decline sharply with the increase of trading volume, while on the other hand, it will rise sharply for sufficiently low trading volume.

To further illustrate the effect of size on the boil off ratio at the achievable design point, Fig. 3 shows the change of daily the boil off ratio with the total heat transfer coefficient and spherical tank volume. The ambient temperature is set at 288 K. For simplicity, it is assumed that the storage tank contains 1.2 bar (a) of pure saturated liquid to hydrogen, and the air-fuel rate is 10%. The liquid density and temperature can be

**Fig. 2.** Relation between internal surface area and volume for cylindrical and spherical tanks

retrieved respectively from the thermophysical properties in REFPROP [21], and the specific heat of evaporation (*h*evap. LH2) is 443.17 kJ/kg. For any fixed value of U, the daily the boil off ratio decreases significantly with the increase of tank size. The daily boil off ratio decreases by about 54% when the volume increases by 10 times (e.g. for 1000 m3 relative to 100 m3). Apart from these correlations based on equation, the rough data of three spherical tanks existed of different sizes and vintage [22], and the tank performance of other two indicators [13, 18] are also plotted in the same chart. Reference [6–8] is also included. The estimated value of the total heat transfer coefficient U of each tank can be read from the curve intersected with the data point, and it indicates that [12, 13, 17] is about 0.004 W m−2 K−1. In theory, if the storage tank [13, 14] or [17] can be further scaled with U unchanged, the boil off ratio would be reduced to about 0.07% per day with a volume of 10,000 m3. The corresponding figure of [22] would eventually reach about 0.17% per day.

**Fig. 3.** Evaluation for the relation between spherical LH2 tank volume, overall heat transfer coefficient and daily BOR.

As these illustrative examples show, due to the reduction of surface-to-volume ratio, the increase of size and diameter is usually beneficial to the low-voltage insulation system [23]. It is still to be identified what is the best technical and economic boil off ratio of large LH2 transport carriers in the future [24], and its size is equivalent to the current LNG tanker. This is an overall research and development task. It is necessary to weigh tank design and thermal insulation layout with a series of trade-offs, some of which may relax some restrictions, such as balancing boil off ratio and energy demand of propulsion and auxiliary systems. Multidisciplinary capabilities at different levels are required, including construction and materials technology, thermodynamics and fluid dynamics, mass and heat transfer, thermal process, naval architecture, and power and propulsion systems. In the research project, "LH2 Pioneer—Super Insulated Marine Containment System for Global LH2 Ship Transportation" led by SINTEF, LH2 containment and thermal insulation, cargo loading and marine hydrogen re-liquefaction process are the key research topics [25].

#### **3.2 Membrane Type Cargo Tank**

Therefore, in order to determine the thermal insulation performance of membrane tank insulation, it is assumed well insulated membrane cargo tank with a uniform temperature distribution both inside and outside of insulation layer. Excluding insulation mounting members, piping, support members and cargo handling systems, only the thermal insulation and heat transfer around cargo tank which volume is assumed to be cubic are considered. The heat flowing into the cargo tank will be proportional to the difference between the outer surface temperature of the insulation layer and the LH2 temperature (K), the outer surface area and the thermal conductivity, as represented in the Eq. (3).

$$Q\_{\rm in} = \lambda \cdot \frac{(T\_{\infty, \rm out} - T\_{\infty, \rm LH\_2})}{D\_{\rm thickness}} \cdot A\_{\rm surface} \tag{3}$$

$$\text{BOR}(\%/\text{day}) = \frac{\text{Qin}}{L\_{\text{evap,LH\_2}} \cdot \rho\_{\text{LH\_2}} \cdot V\_{\text{tan k}}} \cdot 24 \cdot 3600 \cdot 100\% \tag{4}$$

This paper only considers heat conduction for the outer surface temperature of the insulation layer of 0 °C and the inner temperature of − 253 °C, the volume of the cargo tank varies from 1000 to 60,000 m3, and the thickness of the insulation layer is between 0.5 and 2.5 m. Due to the different areas of the hot and cold ends, the average value of the insulation outer area is considered as the outer area of the cargo volume for the conduction heat input calculations and following calculation formula is shown in (4). *D*thickness and λ represent the thickness (m) of the insulation layer and thermal conductivity (W m−1 K−1).

Through the summary of the above literature, it is reasonable to assume that the boil off ratio of the liquid hydrogen cargo tank which is supposed to cube cargo tank is set at 0.2% per day. Figure 4 shows the requirements for the thermal conductivity of the thermal insulation layer in liquid hydrogen ships with different capacities and thicknesses of the thermal insulation layer. As mentioned above, since this is a calculation for a fully enclosed cargo tank that does not include other heat input components, it is necessary to apply an insulation with a lower thermal conductivity than that shown in this result when other heat leaking components are considered. It is obvious from the figure that with the increase of the volume of the cargo hold, the thermal conductivity gradually increases but its increase rate gradually decreases, which means that when the volume of the cargo tank reaches a certain value, the benefits obtained by improving the thermal insulation performance of the insulation material are significantly reduced and the benefits are more pronounced with the thickness of the insulation. However, increasing the thickness of the insulation layer will affect the usable space of the cargo tank and bring about the problem of lower loading economy. Therefore, this is a pair of contradictions and need to find the best balance point.

**Fig. 4.** Thermal conductivity of thermal insulation materials required for liquefied hydrogen tanks with different volumes

Table 1 shows the conductivities of insulation materials applicable to tanks of cryogenic cargo which are commonly used in current industry. Among the proposed materials, VIP is the only material with λ < 0.01 W m−1 K−1 and MLI exists as a better heat-insulating property than VIP. But considering the vacuum structure transformation of the membrane insulation system, it is judged to be technically difficult to apply.

**Table 1.** Thermal conductivity of insulation material

When VIP is applied to the inter barrier space (IBS) of the liquefied hydrogen membrane tank, the temperature is lower than that of the IBS of the existing LNG cargo tank, so there is a difference in operating concept. The IBS of the existing LNG cargo tank was operated in the form of nitrogen purging, but the IBS of the hydrogen cargo hold cannot perform nitrogen purging because there is an area where the temperature is lower than the freezing point of nitrogen (− 210 °C). A method of applying helium instead of nitrogen is also possible.

It is possible to apply a vacuum to the IBS space, but it is necessary to develop and verify the technology for applying/maintaining a vacuum in a large space. In the case of an LNG membrane tank, a vacuum of − 800 mbar is applied through a global test to confirm the airtightness of the primary barrier, so it is expected that a certain level of vacuum application is structurally possible. However, in a liquefied hydrogen carrier, it is necessary to verify whether the vacuum can be maintained during the life cycle of the ship and whether freezing of nitrogen or the like occurs at the corresponding level of vacuum. When vacuum is applied to the IBS, the effects of air condensation and oxygen enrichment due to air inflow in case of vacuum loss should be considered. In addition, the correlation between the reduction in vacuum and the amount of BOG generated due to VIP aging should be considered. Figure 5 is an example of an insulation system arrangement for a membrane tank.

**Fig. 5.** Example of an insulation system for LH2 membrane tank in the future.

## **4 BOG Handling System and Propulsion System**

The development of BOG handling and propulsion systems needs to be considered together with the tank insulation system. In other words, it is anticipated that it will be necessary to develop an insulation system that will control BOG generation to levels required for propulsion/power generation. Furthermore, re-liquefaction of BOG and combustion via Gas Combustion Unit (GCU) is also possible. However, in the case of re-liquefaction, a large amount of energy is consumed compared with the re-liquefaction of LNG, so its application is limited and it is unrealistic to apply it in the short to medium term [28]. Example of BOG handling system and propulsion system in illustrated in Fig. 6.

LNG-based rules prohibit the release of BOG into the atmosphere except in emergency situations because air pollution is closely related to the greenhouse effect and methane emissions [29]. In the case of hydrogen, there is no risk of air pollution, so it is considered necessary to establish a certain area around the exhaust point as a safety

**Fig. 6.** Example of BOG handling system and propulsion system.

zone to eliminate the risk of fire/explosion, and to examine whether natural discharge after dilution can be performed without a separate BOG treatment facility.

Hydrogen can be usually used in the field of fuel cells (see Table 2). Due to the absence of moving parts, fuel cells have some advantages like very low noise, low vibration and low pollutant emissions. However, the tolerance to impurities are their challenges, especially the proton exchange membrane fuel cells (PEMFCs) and shock resistance. PEMFCs and solid oxide fuel cells (SOFCs) are both effective at about 60%. Besides, SOFCs has an additional advantage [30]. They can operate under high temperature (700– 1000 °C), which implies that they can be used to generate steam and power needed in steam turbine. The overall system efficiency can be increased to about 80% [31]. As for the both technologies, the additional ancillary components (e.g. plant balancing) reduce the overall efficiency of a few percent. These losses become higher as long as the fuel cells become larger. The specific power (kg/kW) of SOFCs is lower than that of PEMFCs. They have high operating temperatures, long start-up times and poor tolerance to load variations. Internal combustion engines (ICEs) become more efficient owning to the larger sizes, which have a higher average power density, lower costs as well as more tolerance to load changes, and lasts longer. Some disadvantages of them can be noticed such as noise, vibration and low efficiency.

## **5 Conclusions**

Trade of hydrogen between countries will be carried out in the form of ammonia, liquid hydrogen, LOHC, etc., taking into account the renewable energy resources of the exporting country and the hydrogen usage type and technological maturity of the importing country, and will not be traded in only one form. Therefore, in this article, the key technology or potential demand technology are listed.

(1) Development of efficient insulation system.

• For type C tank, vacuum and MLI/glass bubble insulation technology need to develop application technology.

**Table. 2.** Comparison between direct hydrogen use in fuel cells and ICEs

	- Large-capacity fuel cell (PEM or SOFC) needs to be developed.
	- Hydrogen BOG dilution emission method needs further verification and inspection.

At the beginning of the introduction of liquid hydrogen carriers, about 20,000– 40,000 m3 capacity liquid hydrogen carriers equipped with type C tanks with relatively low technical hurdles are suitable, and in the long term it is considered necessary to develop membrane-type liquid hydrogen carriers. In order to succeed in the commercialization of liquefied hydrogen carriers, the reliability and safety of the technology have to be ensured, and there are still technical challenges to be overcome. It is expected that domestic shipyards will successfully enter into the technology of liquefied hydrogen carriers by utilizing their know-how in developing LNG carriers, which is their strength, and that liquefied hydrogen will play an important role in the hydrogen supply chain.

## **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Influence of Blade Geometry on Performance of Hydrogen Vortex Blower in Fuel Cell System**

Hongshuai Li, Lei Tan(B) , and Huanxin Zhao

State Key Laboratory of Hydroscience and Engineering, Department of Energy and Power Engineering, Tsinghua University, Beijing 100084, China<EMAIL_ADDRESS>
**Abstract.** Hydrogen fuel cell has great potential in replacement of traditional fossil energy systems to decrease carbon dioxide emission. Vortex blower is a key device in the hydrogen recirculation system, which need to be studied deeply to improve the performance of the whole fuel cell. In this paper, the steady internal flow of a hydrogen vortex bower was numerical simulated, and the effect of blade number and blade flapping angle on the performance was studied. The simulation results were compared with experimental data, and the deviation of simulation in choking condition was observed. With the validated simulation method, the influence of blade number and flapping angle was studied. Higher blade number causes more friction, and less blade number leads to flow separation. The negative flapping angle also has the effect on depressing low-pressure region. This research illuminates the simulation method can be further applied to the aerodynamic study and structure optimization of vortex blower.

**Keywords:** Hydrogen vortex blower · Experiment · Blade number · Angle

## **1 Introduction**

Hydrogen fuel cell is considered as a promising option in replacement of traditional fossil energy systems with high efficiency, low emission, and short refueling time. In the automotive industry, polymer electrolyte membrane fuel cell (PEMFC) is used widely. To ensure the operation stability and the duration of membrane, adequate hydrogen is critical for the anode, and same to the cathode with air. The hydrogen ratio defines the ratio of hydrogen inlet mole flow and mole flow consumed by fuel cell stack and it is usually greater than 1. Therefore, hydrogen recirculation system is required to reuse the excess hydrogen from anode exhaust gas, which improves the total efficiency of PEMFC and minimizes environmental pollution [1]. There are many researchers [2, 3] focusing the control strategy and control system simulation of hydrogen recirculation system.

Vortex blower is the key device in hydrogen recirculation system, which inputs work to anode exhaust hydrogen, and it is one of the most energy-consuming devices in PEMFC. Hence, a high efficiency vortex blower is desired for a superior performance PEMFC. Many researchers analyze vortex blower through experimental analysis, theoretical analysis, and simulation analysis. Song [4] discussed about the experimental characteristics of forward and backward bending impeller, which showed little difference in large flow condition. Badami [5] proposed a theoretical model with the consideration of clearance to study leakage. Han [6] analyzed the effect of different blade thickness and blade number, and also predicted fatigue life of blade through simulation. Cai and Zhang [7, 8] studied the influence of blade bending direction and bending angle on flow characteristic of vortex blower with simulation method. In addition, many researchers focus on aero-acoustic and noise simulation [9–11]. This study simulated the flow characteristic of a compact hydrogen vortex blower with experimental data.

## **2 Numerical Simulation**

### **2.1 Geometry Modelling**

As shown in Figs. 1 and 2, the vortex blower consists of two main parts, stational flow path with inlet and outlet, and impeller with rotational flow path. The rotational speed of impeller is 10,000–23,000 rpm. The diameter of inlet and outlet is 5 mm. The impeller inner diameter is 27 mm and outer diameter is 76 mm. The impeller has 36 blades. The thickness of blade is 0.8 mm. The fluid domain was calculated by interference operation. To ensure the full development of the flow, extension is added to both inlet and outlet, with length as same as the impeller outer diameter of 76 mm.

**Fig. 1.** Stator of vortex blower.

## **2.2 Mesh Generation**

Poly-Hexcore meshes were generated separately on stational flow domain and impeller domain by Fluent Meshing 2020 R1, Fig. 3 shows the Poly-Hexcore meshes of impeller domain and stational flow domain. Mesh independency was checked by refining the impeller domain and checking the pressure difference between inlet and outlet. The inlet flow of mesh independency checking state was set as 12.60 m3/h with rotational speed 17,000 rpm. Four groups of impeller domain mesh were generated with 921,221, 1,228,781, 3,872,419, and 5,853,783. The static pressure differences between inlet and outlet were simulated, and compared with the result from experiment (see Table 1). By consideration of both accuracy and computational resource, mesh group 1 was finally chosen for following numerical simulation.

**Fig. 2.** Impeller of vortex blower.

**Fig. 3.** Poly-hexcore meshes of impeller domain and stational flow domain.

#### **2.3 Solver Setup**

In this study, all numerical simulations were conducted by Fluent 2020 R1. To achieve temporal accuracy and numerical stability in transient simulation, Courant number *Co* is required less than 1. The Courant number is defined as:

$$Co = \frac{\delta t |\mathbf{U}|}{\delta x} \tag{l}$$

where δ*t* is the where is the time step, |**U**| is the magnitude of the velocity and δ*x* is the cell size. Because the rotational speed of the vortex blower is at the scale of ten thousand, according to *Co*, the time step should be chosen little than 10−5 s, which will cause enormous computational effort. Therefore, in this study, simulations were completed with MRF method. The interface between the stator and the impeller is not conformal. In Fluent, there are several methods dealing with non-conformal interface, periodic boundary condition, periodic repeats, coupled wall, matching, and mapped. Because the couple wall interface has no flow exchange between two sides of interface, it is not suitable for the simulation in this study, and neither is mapped method, which is always combined with coupled wall.

The viscous model was selected as realizable *k-*ε two equation model, which is suitable for high Reynolds number flow. With *k-*ε model, near-wall treatment method should be selected. Scalable wall functions, non-equilibrium wall functions, and enhanced wall treatment are all suitable for this kind of simulation. By comparison with experimental data, enhanced wall treatment was finally chosen. Coupled scheme was used for pressure velocity coupling, with all velocity-relative terms in second order upwind scheme. For better convergence, pseudo transient option was chosen. As to boundary conditions, mass flow inlet and pressure outlet were selected.

## **3 Results and Discussion**

#### **3.1 Internal Flow Simulation**

The flow characteristic of the vortex blower was obtained by four groups experiments on a standard blower test rig. In every group of experiments, the valve opening was set as a constant, and by increasing impeller revolution speed from 10,000 to 23,000 rpm, the mass flow and pressure difference were recorded.

The internal flow simulation was conducted based on the experiments. The rotational speed and inlet mass flow were set as same as experimental data, and the pressure difference was computed numerically. The error was acceptable in engineering with less than 20%. Hence the validation of internal flow simulation was checked. In Fig. 4, four groups of simulation results and experimental data were plotted. From group 1 to group 4, the valve opening was increased. Black squares represent experimental data and red circles are simulation results. From left to right on x axis, the rotational speed is increased from 10,000 to 23,000 rpm.

From Fig. 4(a)–(d), the deviation between simulation results and experimental data is growing. The reason for this deviation is that the mass flow was increasing from group 1 to group 4, and the vortex blower was operated with over inlet flow outside the design region, which caused choking condition. This is also a characteristic phenomenon of vortex blower, that it is more sensitive to choking condition than surging condition.

In Fig. 5, the pressure contours of interface between stator and rotator was plotted. The fluid is pressured sequentially by blades with the high-speed rotation of impeller, and the blade structure has influence of the flow pattern on the interface. The flow pathlines were shown in Fig. 6 and a series complex eddy structures near inlet were observed with low velocity. The low-energy flow was gathered around the inner circle of impeller, which will cause flow separation. Because of the centrifugal force, the flow separation

**Fig. 4.** Experimental data and simulation results of the vortex blower. Black squares represent experimental data. Red circles represent simulation results. The rotational speed is increased from 10,000 to 23,000 rpm from left to right on x axis.

near outer circle region was depressed. In the Fig. 7, the pressure contour of impeller also showed pressure gradient both from inlet to outlet and inner circle to outer circle.

**Fig. 5.** Pressure contours of interface between stator and rotator.

**Fig. 6.** Pathlines colored by velocity magnitude.

**Fig. 7.** Pressure contour of impeller.

#### **3.2 Influence of Blade Parameters**

In following simulation, the operation point with rotational speed of 17,000 rpm and inlet mass flow of 12.60 m3/h was selected to study the influence of blade parameters to the performance of the vortex blower.

Firstly, the effect of blade number was studied. With the blade number changing evenly from 34 to 42 by 2, six new impeller models were generated. Under the selected operation point, the pressure difference between inlet and outlet was simulated as shown in Fig. 8. With the decrease of blade number, the pressure rise of the vortex blower is increasing. With the same revolution speed, the more blade number means the more friction between fluid flow and blades, which explain the tendency. Besides, the more blade number makes choking condition more easily to happen. So, higher blade number does not benefit the internal flow of vortex blower.

The pressure contours of the interface between stator and rotator are shown in Fig. 9. Blade number is 34 on the left. Blade number is 42 on the right. The main flow patterns with blade structure are observed on both sides. However, the pressure contour of 34 blade-number impeller is much more uneven compared with 42-blade-number impeller, though the pressure rise is higher.

**Fig. 8.** The pressure rises of vortex blower with different blade number at rotational speed of 17,000 rpm and mass flow of 12.60 m3/h.

**Fig. 9.** Pressure contours of interface between stator and rotator. Blade number is 34 on the left. Blade number is 42 on the right.

Figure 10 shows that, in the 34-blade-number impeller, there is a large low-pressure region at the suction surface near inlet, which means sever flow separation. With flow separation, larger secondary flow and other eddy structures will lead to higher energy loss. To ensure the performance of vortex blower with less blade number, further experiments should be done, and there will be an optimal number of blades for specific vortex blower structure.

Secondly, the influence of blade flapping angle was simulated. Blade flapping angle is defined as the space angle between rotation axis and blade plate. In addition, the blade flapping angle is positive when the flapping direction is the same as the rotational direction. In the baseline impeller, all blade plates are parallel to rotation axis, and the blade flapping angle is 0°. New impeller models with blade flapping angle of 5° and − 5° were shown in Fig. 11 left and right. The rotational direction of impeller is counterclockwise face to the paper.

These two impeller models were simulated and compared with baseline impeller as shown in Table 2. The positive flapping angle got higher pressure rise than baseline, and the negative flapping angle got lower pressure rise than baseline. By observing the

**Fig. 10.** Pressure contours of suction surface near the inlet. Blade number is 34 on the left. Blade number is 42 on the right.

**Fig. 11.** Impeller models with blade flapping angle. Left: 5°. Right: − 5°.

pressure distribution on interface in Fig. 12, the positive blade flapping angle impeller has similar flow structure as low blade number, and the negative blade flapping angle influences flow like high blade number. The reason for negative flapping angle impeller has more smooth pressure contour is that the structure of the suction surface depressed the flow separation relatively. For further optimization of vortex blower, the optimal blade number combined with proper blade flapping angle is worth to discuss.

**Table 2.** Simulation results of impellers with different blade flapping angle.

(b) Blade Flapping angle 5° (c) Blade Flapping angle −5°

**Fig. 12.** Pressure contours of interface between stator and rotator.

## **4 Conclusions**

In this paper, a simulation method of the steady internal flow in a vortex blower for hydrogen recirculation system was detailly illustrated. Multiple reference frame method was used to dealing with the relative motion between rotor and stator. The energy characteristic of the vortex blower was recorded through experiments. With the validated simulation method, the influence of blade number and blade flapping angle was discussed. The following conclusions can be drawn:

172 H. Li et al.

• Positive blade flapping angle improves the pressure rise of the vortex blower. Negative blade flapping angle has the influence of depressing the flow separation at the impeller and make pressure field more even.

**Acknowledgments.** This work has been supported by the National Key Research and Development Program of China [2020YFB1901401], the Huaneng Group science and technology research project [KTHT-U22YYJC07], the Guoqiang Institute of Tsinghua University [2021GQG1003], the State Key Laboratory of Hydroscience and Engineering [2021-KY-04, sklhse-2022-Iow06], the Creative Seed Fund of Shanxi Research Institute for Clean Energy of Tsinghua University.

## **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Relationship Between Stress Distribution and Current Density Distribution on Commercial Proton Exchange Membrane Fuel Cells**

Weichen Ji and Rui Lin(B)

School of Automotive Studies, Tongji University, 4800 Cao'an Road, Shanghai, China<EMAIL_ADDRESS>
**Abstract.** The current density of proton exchange membrane fuel cells (PEM-FCs) is directly linked to their electrochemical reaction. Its distribution over the active area can give the local performance of the cells, which is significant for exploration of internal process and optimization of performance. In this paper, segmented cell technology is applied to investigate the current density distribution for a commercial PEMFC with different clamping strategies. The stress distribution and current density distribution as well as the overall performance of the cell are tested under the same operating conditions. The results show that a more uniform stress distribution can lead to a more uniform reaction current density distribution and the good uniformity of the stress distribution and current density distribution has a positive impact on the improvement of the cell overall performance. Thus, it is significant to improve the clamping strategy in order to improve the uniformity of the stress distribution and reaction current density distribution, which ultimately improves the cell overall performance.

**Keywords:** Commercial proton exchange membrane fuel cells · Clamping strategies · Stress distribution · Contact form · Current density distribution

## **1 Introduction**

The PEMFC is regarded as one of the most promising power sources for vehicles and various portable electronic applications due to its high energy conversion rate, high reliability and low pollution emission in operation [1]. Currently, commercial fuel cells usually need hundreds of cells and large active area for each single cell to realize high power applications [2]. Scaling up the active area from no more than 50 cm2 in laboratory level to over 150 cm2 for commercial use encounters challenges in reaction gas flow distribution, water-heat management, and reaction uniformity [3]. In order to optimize the design of commercial PEMFC, it is necessary to carefully consider the coupling relationship of various uneven parameters such as temperature, reaction gas concentration, stress distribution and current density distribution [4].

Conventional fuel cell diagnostic methods such as polarization curve can give an overall performance but lack of spatial information at different positions for the large reaction active area. As current density is directly linked to electrochemical reaction of fuel cells, its distribution over the active area can give the local performance of the cells [5]. Thus, relevant diagnostic methods need to be developed to obtain the current density distribution in the reaction process of commercial PEMFCs with large active area. So far, various researches have been conducted to investigate the current density distribution inside PEMFCs with large active area. In order to study the influence of relative humidity (RH) on the performance and durability of reversal-tolerant-anodes (RTAs) during hydrogen starvation, Wang et al. [6] employed an advanced segmented technique to examine the coupling reactions by simultaneously measuring current density, RH and temperature in a PEMFC with large active area. It is found that the inlet of RTAs undergoes degradation earlier than the outlet and the membrane electrode assembly (MEA) with a RTA has an optimal humidity during cell reversal. Lin et al. [7] studied the changes of the overall performance and current density distribution of commercial PEMFCs under dynamic gas operating parameters especially transient process through fast segmented cell technology. The results showed that the cathode stoichiometric ratio has a more significant effect on the cell performance than the back pressure. The work leads to suggestions on the control strategy during dynamic load demands and lays a foundation for further study into the durability under variable gas operation parameters of the fuel cell operation.

Moreover, some researches have been conducted to study the relationship between stress distribution and current density distribution as well as the cell overall performance. Vijayakumar et al. [8] measured the current distribution profile along the cathode flow field channel using a segmented current measurement plate to understand the influence of uneven clamping pressure distribution on the MEA. Although scholars have done some research on the relationship between the stress distribution and the current density distribution affected by the clamping force, few studies on the stress distribution and its relationship with the current density distribution caused by the cell structure and clamping strategy have been conducted so far inside commercial PEMFCs with large active area.

In this work, two different clamping strategies are applied to a commercial PEMFC with the active area of 360 cm2. The current density distribution is tested by segmented cell technology (from Yuanzhu technology, www.szfcet.com) and the stress distribution is obtained. The results of both strategies are compared not only to analyze the influence of different structures on stress distribution but also to examine the relationship between stress distribution and current density distribution as well as the cell overall performance.

### **2 Experimental Part**

The research is carried out on a commercial single PEMFC with the active area of 360 cm2 with two different clamping strategies (strategy 1 and strategy 2) under the standard operating conditions. The conditions contain: coolant temperature of 75°C, anode/cathode back pressure of 100/95 kPa, anode/cathode stoichiometry of 1.5/2.2, anode/cathode RH of 80%/80%. The segmented cell technology is from Yuanzhu technology (www.szfcet.com), which is applied to test the current density distribution of the fuel cell and the stress distribution is obtained by pressure measurement film.

### **3 Results and Discussion**

#### **3.1 Influence of Different Structures on Stress Distribution**

**Fig. 1.** The stress distribution (**a**) and the flow field structure (**c**) of the contact interface of the cell with strategy 1; the stress distribution (**b**) of the cell with strategy 2; cross section of the parallel areas (**d**) and the vertical areas (**e**, **f**) of the flow field structure.

Figure 1 analyzes the influence of different structures on stress distribution. It can be found that the stress distribution of the cell with strategy 1 (Fig. 1(a)) is less uniform than the stress distribution of the cell with strategy 2 (Fig. 1(b)). As shown in Fig. 1(c), the contact interface of the cell with strategy 1 can be divided into two parts. One is the area where the ribs on the contact interface are vertical to each other (the triangle areas of Fig. 1(c)), and the other is the area where the ribs on the contact interface are parallel to each other (the oval areas of Fig. 1 (c)). Comparing the stress distribution and the flow field structure on the contact interface, we can find that the contact where ribs are vertical to each other is better than that where ribs are parallel to each other. In vertical areas, the ribs and channels on both sides of the interface are interlaced. The contact forms of the interface can be divided into three types: rib to rib, rib to channel, channel to channel (Fig. 1(e, f)). From Fig. 1(a) we can find that both the rib-to-rib form and the rib-to-channel form can achieve good contact stress, while the channel-to-channel form gains poor contact. As most of the vertical areas are rib-to-rib and rib-to-channel forms, these areas generally result in good contact. On the other hand, the parallel areas only have rib-to-rib and channel-to-channel forms. And as the proportion of the channelto-channel form which has poor contact is considerably large (near 50%), the parallel areas generally result in poor contact. As shown in Fig. 1(b), the interface of the cell with strategy 2 only has rib-to-channel and rib-to-rib forms (Fig. 1(e)). Thus, the stress distribution on the entire contact interface is more uniform and generally better.

#### **3.2 Relationship Between Stress Distribution and Current Density Distribution**

**Fig. 2.** The current density distribution of the cell with strategy 1 (**a**) and the cell with strategy 2 (**b**) at 1000 mA/cm2.

Figure 2 shows the current density distribution of the cell with two clamping strategies at 1000 mA/cm2 respectively, which investigates the relationship between stress distribution and current density distribution. As shown in Fig. 2(a), it is found that the local current density in some areas (the oval areas of Fig. 2(a)) is lower than the normal value. It can be found from the stress distribution and flow field structure in Fig. 1(c) that these low current density areas are consistent with the poor contact areas (parallel areas). Therefore, the local poor contact may affect the reaction in these areas and lead to low local current density. On the other hand, by comparing the triangle areas of the stress distribution (Fig. 1 (c)) and the current density distribution (Fig. 2(a)), the results show that local good contact contributes to the reaction in these areas, which leads to the increase of local current density. Moreover, the current density distribution of the cell with strategy 2 (Fig. 2(b)) is much more uniform than that of the cell with strategy 1 (Fig. 2(a)). By combining the phenomenon with the stress distribution of the cells with two clamping strategies respectively (Fig. 1(a) and (b)), we can find that a more uniform stress distribution can lead to a more uniform reaction current density distribution.

#### **3.3 Relationship Between Current Density Distribution and the Cell Overall Performance**

As shown in Fig. 3, the red line is the polarization curve of strategy 2, and the black line is the polarization curve of strategy 1. It can be found that the overall performance of the cell with strategy 2 is better than that of the cell with strategy 1, especially at medium and

**Fig. 3.** Polarization curves of the cell with both clamping strategies.

low current density region. Since the stress distribution and current density distribution of the cell with strategy 2 is more uniform than that of strategy 1 (Figs. 1 and 2), it can be deduced that the well uniformity of the stress distribution and current density distribution has a positive impact on the improvement of the cell overall performance.

## **4 Conclusions**

In this study, two different clamping strategies are applied to a commercial PEMFC with the active area of 360 cm2. The stress distribution and current density distribution as well as the overall performance of the cell in both cases are tested under the same operating conditions. The results show that the vertical areas where rib-to-rib and rib-tochannel forms occupy the majority have good contact while the parallel areas where the proportion of the channel-to-channel form which has poor contact is considerably large have poor contact, which leads to the uneven stress distribution of the cell with strategy 1. The cell with strategy 2 of which the contact interface only has rib-to-channel and rib-torib forms has more uniform and generally better stress distribution on the entire contact interface. By comparing the relationship between the stress distribution and the current density distribution of both strategies, it is found that a more uniform stress distribution can lead to a more uniform reaction current density distribution. From the results of the overall performance of the cells in both strategies, it can be considered that the good uniformity of the stress distribution and current density distribution has a positive impact on the improvement of the cell overall performance. Thus, it is significant to improve the clamping strategy and the contact form of the flow field in order to improve the uniformity of the stress distribution and reaction current density distribution, and ultimately improve the cell overall performance.

**Acknowledgement.** The authors gratefully acknowledge the financial support by Key projects of intergovernmental international scientific and technological innovation cooperation (No. 2022YFE0102900) and the International Cooperation Program of Science and Technology Commission of Shanghai Municipality (No.<PHONE_NUMBER>0).

## **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Hydrogen Leakage Detection in Confined Spaces with Drift Suppression Based on Subspace Alignment**

Haifeng Se, Jinhai Jiang, Chuanyu Sun, and Kai Song(B)

School of Electrical Engineering and Automation, Harbin Institute of Technology, Harbin 150001, China<EMAIL_ADDRESS>
**Abstract.** Gas sensor has been widely used in flammable gas detection. In order to solve the problem that gas sensors are susceptible to drift in hydrogen leakage detection in confined spaces, a drift compensation framework based on subspace alignment is proposed. First, global domains and subdomains are aligned simultaneously in the new subspace, thus reducing the distribution difference between the source domain and the target domain. Then, the proposed method utilizes extreme learning machine (ELM) to iteratively refine the prediction results of the target domain, and continuously optimizes the subspace and classifier. In this way, the proposed method realizes drift compensation at the feature level. Compared with the existing methods, the proposed method achieves the highest accuracy of 79.83% in the long-term drift scenario. Therefore, the experimental results show that the proposed method is competent for hydrogen leakage detection with drift, and can provide a reference for the design of drift compensation method based on gas sensors.

**Keywords:** Gas sensor · Hydrogen leakage detection · Drift compensation · ELM

## **1 Introduction**

Compared with oil and natural gas, hydrogen energy has a higher conversion rate and lower carbon emissions, which is regarded as one of the feasible solutions to solve the problem of fossil fuel environmental pollution and energy crisis [1]. However, hydrogen is a flammable gas and hydrogen leakage in confined spaces may occur in the process of hydrogen transportation and hydrogen storage. In this case, it has potential risks of combustion and explosion. Therefore, it is necessary to design hydrogen leakage methods.

With the development of material science and artificial intelligence, gas sensors have been widely used in the field of explosive detection [2]. Most gas detection methods utilize gas sensors + machine learning to identify flammable gases. They generally assume that the training and test sets are independent and equally distributed. The prediction model is built on the training set and applied to the test set. However, due to aging, temperature and humidity interference, even the same concentration of the target gas, the response of the gas sensor is still different, i.e., sensor drift. In this context, the methods which assume that the training set and the test are identically distributed are not competent for gas detection in the drift scene in practical application.

For gas detection with sensor drift, Zhang et al. [3] proposed a machine learning method based on subspace projection to suppress gas sensor drift. Artursson et al. [4] introduced a drift counteraction method based on principal component analysis and partial least squares, and thus suppress sensor drift by removing drift direction component. Vergara et al. [5] designed a machine learning approach, namely the integration of classifiers, to solve the problem of gas identification over long periods of time. However, the accuracy of existing methods is still low. Therefore, a gas sensor compensation method based on subspace alignment is proposed. This method is suitable for gas identification in drift scene with high accuracy by aligning global domains and subdomains in the new subspace.

### **2 Methodology**

The source domain with *D*-dimension are denoted as *XS* = [*x*1 *s* , *x*2 *s* , ··· , *x NS s* ]∈*D*×*NS* , and the target domain with *D*-dimension are denoted as *XT* = [*x*1 *T* , *x*2 *T* , ··· , *x NT T* ] ∈ *D*×*NT* . *P* ∈ *D*×*d* represents a base matrix. ·2 is the *l*2-norm. *tr*(·) denotes the trace operator.

#### **2.1 Subspace Alignment Method (SAM)**

To suppress sensor drift, the proposed SAM utilizes *P* to achieve domain alignment, and thus reduce distribution divergence between source domain and target domain. First, we focus on the alignment of the global distribution based on marginal distribution distance in the subspace.

$$\min\_{P} \|\mu\_S - \mu\_T\|\_2^2 = \min\_P \left\|\frac{1}{N\_S} \sum\_{l=1}^{N\_S} \mathbf{y}\_S^l - \frac{1}{N\_T} \sum\_{j=1}^{N\_T} \mathbf{y}\_T^j \right\|\_2^2 \tag{1}$$

where μ denotes the centre of the domain in the subspace. *yi S* = *PT xi S* (*i* = 1, 2, ...,*NS* ) and *y j T* = *PT x j T* (*j* = 1, 2, ...,*NT* ) are the projected samples in the subspace. Then we rewrite (3) as:

$$\begin{aligned} &\min\_{P} \left\| \frac{1}{N\_S} \sum\_{i=1}^{N\_S} \mathbf{y}\_S^i - \frac{1}{N\_T} \sum\_{j=1}^{N\_T} \mathbf{y}\_T^j \right\|\_2^2 \\ &= \min\_{P} \left( \mathbf{P}^T (\boldsymbol{\mu}\_S - \boldsymbol{\mu}\_T)(\boldsymbol{\mu}\_S - \boldsymbol{\mu}\_T)^T P \right) \end{aligned} \tag{2}$$

where *u* is the centre of the domain in the original space. For example, *uS* = 1 *NS NS i*=1 *xi S* and *uT* = 1 *NT NT j*=1 *xj T* .

We then consider aligning subdomains with the same category between domains, which is formulated as follows

$$\begin{aligned} &\min\_{P} \sum\_{c=1}^{C} \left\| \frac{1}{N\_{\mathcal{S}}^{c}} \sum\_{\mathbf{y}\_{\mathcal{S}}^{i} \in class} \mathbf{y}\_{\mathcal{S}}^{i} - \frac{1}{\hat{N}\_{T}^{c}} \sum\_{\hat{\mathbf{y}}\_{\mathcal{S}}^{i} \in class} \hat{\mathbf{y}}\_{T}^{\ell} \right\|\_{2}^{2} \\ &= \min\_{P} tr \left( \sum\_{c=1}^{C} P^{T} (\boldsymbol{u}\_{\mathcal{S}}^{c} - \hat{\boldsymbol{u}}\_{r}^{c}) (\boldsymbol{u}\_{\mathcal{S}}^{c} - \hat{\boldsymbol{u}}\_{r}^{c})^{T} P \right) \end{aligned} \tag{3}$$

where *y*ˆ *j T* represents the *j*-th sample whose classifier prediction result belong to class *c* from the target domain. *N*ˆ *c T* is the number of samples whose classifier prediction result belong to class *c*. *u*ˆ *c T* is the centre of c-class samples from the target domain.

To prevent data distortion from source and target domains, it is necessary to maximize the following term

$$\begin{aligned} &\max\_{P} ((1 - \alpha\_1) \times tr(P^T X\_S X\_S^T P) + \alpha\_1 \times tr(P^T X\_T X\_T^T P)) \\ &= \max\_{P} (tr((1 - \alpha\_1) \times P^T X\_S X\_S^T P + \alpha\_1 \times P^T X\_T X\_T^T P)) \end{aligned} \tag{4}$$

where α1 ∈ [0, 1][0, 1].

Finally, the proposed SAM simultaneously optimizes (2), (3) and (4), which is formulated as:

$$\min\_{P} \frac{\alpha\_{2} \times tr(P^{T} (\boldsymbol{u}\_{\mathcal{S}} - \boldsymbol{u}\_{T})(\boldsymbol{u}\_{\mathcal{S}} - \boldsymbol{u}\_{T})^{T}P) + \alpha\_{3} \times tr(\sum\_{c=1}^{C} P^{T} (\boldsymbol{u}\_{\mathcal{S}}^{c} - \hat{\boldsymbol{u}}\_{T}^{c})(\boldsymbol{u}\_{\mathcal{S}}^{c} - \hat{\boldsymbol{u}}\_{T}^{c})^{T}P)}{tr((1 - \alpha\_{1}) \times P^{T} \boldsymbol{X}\_{\mathcal{S}} \boldsymbol{X}\_{\mathcal{S}}^{T} P + \alpha\_{1} \times P^{T} \boldsymbol{X}\_{T} \boldsymbol{X}\_{T}^{T} P)} \tag{5}$$

To achieve a unique solution, (5) is rewritten as

$$\begin{cases} \min\_{P} (\alpha\_2 \times tr(P^T (u\_S - u\_T)(u\_S - u\_T)^T P) + \alpha\_3 \times tr(\sum\_{c=1}^C P^T (u\_S^c - \hat{u}\_T^c)(u\_S^c - \hat{u}\_T^c)^T P) \\\\ s.t.(1 - \alpha\_1) \times P^T X\_S X\_S^T P + \alpha\_1 \times P^T X\_T X\_T^T P = I \end{cases} \tag{6}$$

where *I* ∈ *d*×*d* is the identity matrix.

After that, the Lagrange function *L* is introduced:

$$L = \alpha\_2 \times tr(\boldsymbol{P}^T (\boldsymbol{u}\_S - \boldsymbol{u}\_T)(\boldsymbol{u}\_S - \boldsymbol{u}\_T)^T \boldsymbol{P}) + \alpha\_3 \times tr(\sum\_{c=1}^C \boldsymbol{P}^T (\boldsymbol{u}\_S^c - \hat{\boldsymbol{u}}\_T^c)(\boldsymbol{u}\_S^c - \hat{\boldsymbol{u}}\_T^c)^T \boldsymbol{P})$$

$$-tr((\boldsymbol{I} - ((\boldsymbol{I} - \boldsymbol{\alpha}\_1) \times \boldsymbol{P}^T \boldsymbol{X}\_S \boldsymbol{X}\_S^T \boldsymbol{P} + \boldsymbol{\alpha}\_1 \times \boldsymbol{P}^T \boldsymbol{X}\_T \boldsymbol{X}\_T^T \boldsymbol{P})) \boldsymbol{\Phi}) \tag{7}$$

where = diag(φ1, φ2,...,φ*d* ) ∈ *d*×*d* denote the Lagrange multipliers.

According to ∂*L* ∂*P* = 0, we can get

$$((1 - \alpha\_1) \times X\_S X\_S^T + \alpha\_1 \times X\_T X\_T^T))^{-1} (\alpha\_2 (\mu\_S - \mu\_T)(\mu\_S - \mu\_T)^T$$

$$+ \alpha\_3 \sum\_{c=1}^C (\boldsymbol{u}\_s^c - \hat{\boldsymbol{u}}\_r^c)(\boldsymbol{u}\_s^c - \hat{\boldsymbol{u}}\_r^c)^T)P = P\Phi \tag{8}$$

*P* can be solved by eigenvalue decomposition based on *A*. i.e.,

$$AP = P\Phi$$

where

$$A = ((1 - \alpha\_2) \times X\_S X\_S^T + \alpha\_2 \times X\_T X\_T^T)^{-1} \left( (1 - \alpha\_1) \times (\mu\_S - \mu\_T)(\mu\_S - \mu\_T)^T \right)$$

$$+ \alpha\_1 \times \sum\_{c=1}^{C} (\mu\_S^c - \hat{\mu}\_T^c)(\mu\_S^c - \hat{\mu}\_T^c)^T \right).$$

Therefor, *P*∗ is formulated as

$$P^\* = [P\_1, P\_2, \dots, P\_d] \tag{10}$$

where *Pk* (*k* = 1, 2, ..., *d*) denotes the eigenvector with respect to the *k*-th minimum eigenvalue.

#### **2.2 The Proposed Framework**

The proposed framework first calculates *P* according to (10), and the source domain and target domain are projected into the subspace, thus reducing the distribution differences between domains. Then, the data of the subspace is input into ELM [6] for training, and thus obtain the prediction result of the target domain. Finally, the target domain is iteratively optimized to predict the results and *P* (see Fig. 1).

### **3 Experiments**

#### **3.1 Dataset**

The benchmark sensor drift dataset was collected over 36-month from UCSD [5], including acetone, acetaldehyde, ethanol, ethylene, ammonia, and toluene. There are 13,910 samples in the data set and each sample consists of 128-dimensional feature vector. The benchmark sensor drift dataset is divided into ten batches, as shown in Table 1.

#### 184 H. Se et al.

**Fig. 1.** Average accuracy with different Iterations.

**Table 1.** The sensor drift dataset.

#### **3.2 Experimental Results**

The experimental environment is an Intel(R) core (TM) i7-10750H<EMAIL_ADDRESS>GHz, 16 GB memory, Windows 11 and Python 3.6. We set batch 1 as the training set (source domain) and batch *t* (*t* = 2, 3, …, 10) as the test set (target domain). The proposed method is compared with several drift compensation methods, including SVM, ELM and DRCA [3].

The experimental results are shown in Table 2. Compared with SVM, ELM and DRCA, the proposed method achieves the highest classification accuracy by reducing the distribution difference. Therefore, these results show that the proposed drift compensation strategy is competent for gas detection in long-term drift scenarios.

**Table 2.** Experiment results.

During training we hope to use as few iterations as possible to achieve the highest accuracy. Figure 2 shows the average accuracy of the proposed method with iterations. The results show that the proposed method is effective and the highest average accuracy is achieved in fewer iterations.

**Fig. 2.** Average accuracy with different Iterations.

## **4 Conclusion**

Hydrogen leakage detection in confined spaces is the guarantee of successful hydrogen transport and hydrogen storage technology. In this paper, focusing on the practical application scenario, a drift compensation method for gas sensor is proposed. In this method, subspace alignment is realized by aligning subdomains and global domains, and thus reducing the distribution difference due to sensor drift. Experimental results show that the proposed method achieves the highest classification accuracy in the long-term drift with an average accuracy of 79.83%.

**Acknowledgments.** This work was supported by the Central Science and Technology Commission of China under grant JZKJW20200036.

## **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Design of 10 L/s Liquid Hydrogen Pump with Integrated Inducer and Centrifugal Impeller**

ZiWei Li1,2, Cui Lv1(B) , Jihao Wu1,2, Linghui Gong1, Jinzhen Wang1, Jin Shang1, and Ming He1

1 Key Laboratory of Cryogenic Science and Technology, State Key Laboratory of Technologies in Space Cryogenic Propellants, Technical Institute of Physics and Chemistry, Chinese Academy of Sciences, Beijing, China

<EMAIL_ADDRESS>
2 Key Laboratory of Cryogenic Science and Technology, University of Chinese Academy of Sciences, Beijing, China

**Abstract.** With the global focus on energy and environmental issues, hydrogen energy development has become more and more important. The liquid hydrogen pump is an important equipment for hydrogen transportation, whose performance is significantly decided by the impeller. In this paper, we designed an integrated impeller with an inducer and an centrifugal impeller and optimized the axial length and outlet angle of the centrifugal impeller. Grid independence and cross validation prove that the simulation model and results are reliable. It is worth noting that a high roughness is given to the impeller surfaces in our simulation to obtain results closer to the real situation. By analyzing the pressure and relative liquid flow angle distribution in the fluid domain, the structural design of the impeller is optimized, and then the performance curve is obtained. The simulation results will provide a guidance for the manufacture of our liquid hydrogen pump, and also provide a reference for the design and manufacture of other liquid hydrogen pumps.

**Keywords:** Liquid hydrogen pump · Integrated impeller · High speed

The development of hydrogen energy has attracted global focus, which is greatly concerned with energy and environmental issues. Recently, many countries around the world have deployed hydrogen energy-related technologies. Due to the advantages of low energy consumption and high efficiency, the liquid hydrogen pump is widely valued for liquid hydrogen transferring. The impeller is the most important part of the pump. Considering the cavitation issue, we used an integrated impeller with an inducer and centrifugal impeller to improve the performance.

At present, there is less literature to be investigated. In 2015, Li et al. [1] designed a semi-open integrated impeller for aviation fuel pumps. The performance requirements are 8000 rpm, 77,000 L/h, and the inlet and outlet pressure difference is 0.956 MPa. According to calculations, the pump does not produce cavitation in the operating range, and the efficiency is 65%. And then, they optimized the profile line to improve the performance [2]. Liu et al. [3] optimized the blade wrapping angle. In 2018, Li et al. [4] designed a fuel centrifugal pump with an integrated inducer and impeller influenced by an inlet flow ejector. The pump performance requirements are 1500 L/h, 10678 rpm, 23.07 m. And the simulation efficiency is 69%. The experiment efficiency is 65%. Tian et al. [5] designed an integrated impeller for an aero fuel pump and optimized the centrifugal impeller blade wrapping angle. The pump performance requirements are 3.9 kg/s, 9000 rpm, -P ≥ 80 kPa. And the simulation efficiency is 81.5%.

In this paper, we designed an integrated impeller with an inducer and centrifugal impeller and optimized axial length and outlet angle. By analyzing the flow field pressure and relative liquid flow angle distribution, the structural design of the impeller is determined, and the performance curve is obtained.

## **1 Preliminary Design of Integrated Impeller**

## **1.1 Performance Parameters**

**Table 1.** Performance parameters

According to Table 1, the liquid hydrogen integrated impeller is designed as shown in Fig. 1. The impeller has 3 main blades and 3 auxiliary blades.

**Fig. 1.** Closed integral impeller

### **1.2 Calculation and Boundary Conditions**

Adopt the SST k − ω turbulence model and secondary windward space discrete format. The convergence residual is set to 10–5. Set the inlet as pressure boundary condition as 1.3 bar and outlet as volume flow rate boundary as 10 L/s. Set wall roughness as 1.2 μm.

### **1.3 Results of Preliminary Impeller Design**

The results of the preliminary design are shown in Table 2. The efficiency is 89.65%. The flow field pressure and velocity distribution are shown in Fig. 2.

**Table 2.** The results of preliminary impeller design

a. Total pressure b. Absolute velocity

## **2 Optimization of the Integrated Impeller Structure**

#### **2.1 Optimization of the Axial Length**

After changing the axial length of the centrifugal impeller from 15 to 40 mm with the interval of 5 mm, we get the calculation results as shown in Fig. 3. When the axial length of the centrifugal impeller is 25 mm, the efficiency is optimal.

**Fig. 3.** Efficiency and head change with the axial length of the centrifugal impeller

**Fig. 4.** Relative flow angle of each plane when axial length of the centrifugal impeller is 25 mm.

The relative flow angle contribution is shown in Fig. 4. 20, 35, 45, and 55 mm correspond to the inducer plane, the inlet, the middle, and the outlet plane of the centrifugal impeller. On each plane, with the increase of radius, the circumferential direction velocity increases, and the relative flow angle decreases. At the 20 mm plane, the relative flow angle is about 35° at the roof of the blade and about 18° at the top of the blade. The relative flow angle is consistent with the β angle of the blade design. At the 35 mm plane, the relative flow angle is about 34° at the roof of the blade and about 26° at the top of the blade. The design β angle is 43° and 23°. The low relative flow angle area at the inner diameter is affected by the leading edge of the auxiliary blade. At the 45 mm plane, the relative flow angle is about 39.5° at the roof of the blade and the design β angle is 41.5°. The relative flow angle is about 42.4° at the top of the blade and the design β angle is 24.5°. The β angle of the blade limits the development of fluid flow and the counter-current zone appears. At the 55 mm plane, the relative flow angle is about 45° at the roof of the blade and the design β angle is 32°. At top of the blade, the flow angle is 25°, which is close to the design of β angle 27°. In general, the impeller structure needs to be improved.

**2.2 Optimization of the Outlet Angle**

The outlet angle affects the pump performance including hydraulic efficiency, and characteristic curve. Six groups of outlet angles are shown in Table 3. And the calculation results are shown as follows.

**Table 3.** Impeller performance at different outlet angles

**Fig. 5.** Relative flow angle distribution at different outlet angles

As shown in Table 3, when selecting a larger outlet angle such as ➅, the efficiency is highest under the design working conditions. Due to 2.1, the area where the blade angle is not reasonable is between the length of 45 and 55 mm. And the distribution of liquid flow angle at the shaft length section of 45 mm and 55 mm under the condition of six groups of blade shapes was analyzed.

As showed in Fig. 5, in group ➀, at 45 mm plane, the liquid flow angles at the root and top of the leaf at 45 mm were about 35° and 24°, and the design β angles of the leaves were 38° and 23°, respectively. The design angle is appropriate. At 55 mm plane, the liquid flow angle at the root and top of the leaf is about 48° and 32°, and the design β angles of the blade are 28° and 23°. It can be seen that the outlet angle setting is too small. In group ➁, at 45 mm plane, the design angle is appropriate. At 55 mm plane, the liquid flow angle at the root and top of the leaf is about 40° and 23°, and the design β angles of the blade are 30° and 25°. It can be seen that the outlet angle setting is small but has improved compared to ➀. ➂, ➃, ➄ is also gradually improving compared to the previous one. In group ➅, the liquid flow angles at the root and top of the leaf were about 43° and 29° at 45 mm plane, and the design β angles of the leaves were 43° and 28.5°. At 55 mm plane, the liquid flow angle at the root and top of the leaf is about 38° and 33°, and the design β angles of the blade are 38° and 34°. The design outlet angle is appropriate. So we adopt the group ➅.

### **3 Results**

**Fig. 6.** Impeller performance curve

The Impeller performance curve is shown in Fig. 6. Under the design flow rate condition of 10 L/s, the impeller efficiency reaches the optimal value of 92.63%, and the head is 167.24 m (Fig. 7).

**Fig. 7.** Total pressure and absolute velocity of the final design

Finally, the complete structure of the impeller and volute under the design condition was calculated, and the efficiency was 85.69%. The efficiency loss in volutes was 6.9%.

## **4 Conclusion**

The liquid hydrogen pump impeller is designed with an integrated inducer and centrifugal impeller in this paper. After the axial length of the centrifugal impeller and are outlet angle are optimized, simulation is performed on the impeller with the design conditions. Even with the high surface toughness, the impeller efficiency and overall efficiency with the volute shell reach 92.63% and 85.69%, respectively. With the reliable results shown in this paper, our liquid hydrogen pump will build up and test soon, and the test performance will be published. This paper will benefit for the development of the low energy consumption and high efficiency liquid hydrogen pump.

**Acknowledgments.** This work was supported by the National Science Foundation of China (Grant No. 52106034), the Youth Innovation Promotion Association Innovative, Chinese Academy of Sciences, and the Talent Cultivation Project of Technical Institute of Physics and Chemistry, Chinese Academy of Sciences (Grant No. 2021-LC).

## **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Experimental Research on High-Pressure Hydrogen Leakage and Diffusion of Hydrogen Refueling Station**

Bin Jiang(B) , Tongshen Zhen, and Fangfang Fang

State Power Investment Corporation R&D Institute, Changping Beijing 102209, China<EMAIL_ADDRESS>
**Abstract.** As the infrastructure to provide hydrogen for hydrogen fuel cell vehicles, hydrogen refueling station is a very important part of hydrogen energy utilization. However, due to the characteristics of hydrogen, such as flammability and explosion, low density, wide range of explosive limit concentration, hydrogen refueling station accidents often occur. The existing research on hydrogen refueling stations often uses the method of numerical simulation, and mainly considers the leakage of hydrogen storage tank. There are few relevant experimental studies and little consideration is given to the case of hydrogen pipeline leakage. In order to explore the phenomenon and rule of high-pressure hydrogen leakage and diffusion in the pipeline of hydrogen refueling stations, a full-size high-pressure hydrogen leakage test facility is built based on a real hydrogen refueling station. The vehicle-mounted high-pressure hydrogen storage tank is used as the high-pressure hydrogen gas source to provide constant hydrogen pressure to the test section through the combination of different valves in the pipeline and the instrument control system. By changing different leakage sizes and pressures, the concentration distribution and influence factors after hydrogen leakage are analyzed, which provides an important basis for the optimal layout and operation and maintenance of the safety facilities of the existing hydrogen refueling station.

**Keywords:** Hydrogen leakage · Hydrogen refueling station · Test facility

## **1 Introduction**

Hydrogen refueling station as the infrastructure to provide hydrogen for hydrogen fuel cell vehicles, hydrogen refueling station is a very important part of hydrogen energy utilization. However, due to the characteristics of hydrogen, such as flammability and explosion, low density, wide range of explosive limit concentration, hydrogen refueling station accidents often occur. The existing research on hydrogen refueling stations often uses the method of numerical simulation, and mainly considers the leakage of hydrogen storage tank. There are few relevant experimental studies and little consideration is given to the case of hydrogen pipeline leakage. In order to explore the phenomenon and rule of high-pressure hydrogen leakage and diffusion in the pipeline of hydrogen refueling stations, a full-size high-pressure hydrogen leakage test facility is built based on a real hydrogen refueling station. The vehicle-mounted high-pressure hydrogen storage tank is used as the high-pressure hydrogen gas source to provide constant hydrogen pressure to the test section through the combination of different valves in the pipeline and the instrument control system. By changing different leakage sizes and pressures, the concentration distribution and influence factors after hydrogen leakage are analyzed, which provides an important basis for the optimal layout and operation and maintenance of the safety facilities of the existing hydrogen refueling station.

## **2 Experimental System**

#### **2.1 Experimental System Design**

Based on the design of pipelines and valves for typical hydrogen energy facilities and equipment, a high-pressure hydrogen storage tank (with a pressure more than 35 MPa) is used as the high-pressure hydrogen gas source. through the combination of different valves in the pipeline and the instrument control system, a constant hydrogen pressure is provided to the pipeline in the test section, and a high-pressure hydrogen gas leakage diffusion test facility is built.

Using hydrogen concentration sensors, sound sensors, fiber optic sensors, and infrared sensors to monitor the diffusion behavior of hydrogen after high-pressure leakage, analyze the flow field structure and concentration spatiotemporal distribution of hydrogen under different leakage pressures and leakage ports, analyze the diffusion behavior, acceleration, and changes in hydrogen leakage vibration measured by sound sensors, and comprehensively analyze the impact of different fracture sizes on hydrogen diffusion behavior, Provide data support for the hydrogen leakage diffusion dynamics model and the prediction method of combustion and explosion risk areas.

#### **2.2 Test Facility**

The high pressure hydrogen leakage test facility is shown in Fig. 1.

The system mainly consists of a gas supply system, a pressurization system, a test section, and a relief system. Vibration sensors are arranged on the test section, sound sensors are arranged near the leakage port, and several hydrogen concentration sensors are arranged above the leakage port to monitor the leakage behavior of hydrogen. The actual test facility is shown in Fig. 2.

The experimental system is equipped with two leakage nozzles with diameters of 0.05 and 1 mm. The nozzle can be horizontal and vertical. The hydrogen booster can pressurize the hydrogen gas to a maximum of 35 MPa. The input parameters of this experimental system mainly refer to the relevant parameters of the actual demonstration hydrogen refueling station, as is shown in Table 1.

The material of the test section pipeline is Q345 or 316 L seamless steel pipe, and the maximum size of simulated leakage in the test section is the inner diameter of the pipeline. According to the actual operation of a hydrogen refueling station, the pipeline size is 3/4 inch, with an outer diameter of 19.05 mm, which corresponds to an inner diameter of 10.97 mm (0.432 in) at a working pressure of 35 MPa.

**Fig. 1.** Schematic diagram of the test facility

**Fig. 2.** High pressure hydrogen leakage and diffusion test facility

**Table 1.** Overall parameters of the test facility

## **3 Test Results**

## **3.1 Horizontal Hydrogen Ejection Condition**

Horizontal hydrogen ejection condition is shown in Fig. 3. Three sensors were arrange at a distance of 1m downstream of the nozzle, the lower sensor is located directly in front of the nozzle, and the distance between the upper two sensors is 50 cm. The pressure at the nozzle is 52 bars. Within 2 s after the start of hydrogen injection, the hydrogen sensor measures the hydrogen concentration value. The hydrogen sensor facing the nozzle measures the maximum value, and the middle hydrogen sensor delays for 1 s to measure the hydrogen concentration data. The top hydrogen sensor hardly measures the hydrogen gas data. When the hydrogen injection ends, the wake concentration is measured to be about 1%.

## **3.2 Vertical Hydrogen Ejection Condition**

Vertical hydrogen ejection condition is shown in Fig. 4. Three sensors are arranged. The middle sensor is located directly above the nozzle, and the distance between the sensors on both sides is 70 cm from the middle sensor. The pressure at the nozzle is 42 bars. Within 2.5 s after the hydrogen gas is sprayed, the hydrogen sensor measures the hydrogen concentration value, while the middle hydrogen sensor measures the maximum value. The two symmetrically arranged side sensors measure a smaller value, and the maximum hydrogen concentration values of the two are almost the same.

## **4 Conclusions**

Further research need to be done such as obstacles, different nozzle shapes and higher pressures (700 bars).

**Fig. 3.** High pressure hydrogen leakage—horizontal injection

**Fig. 4.** High pressure hydrogen leakage—vertical ejection

## **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Coupled Heat Transfer Characteristics of SiC High Temperature Heat Exchanger in Solid Oxide Fuel Cell**

Jia-Liang Song1(B) , Hua Chen2, Yong-Dong Chen1, Gai-Ge Yu1, Hong-Wei Zou1, and Bing-Chuan Han1

1 Hefei General Machinery Research Institute Co., Ltd. (HGMRI), Hefei 230031, China<EMAIL_ADDRESS><EMAIL_ADDRESS>
2 University of Science and Technology of China, Hefei 230027, China

**Abstract.** High temperature heat exchanger is a crucial equipment in the BOP of SOFC. Replacing the commonly used metal materials with high-temperature resistant SiC ceramic materials for the manufacturing of SOFC high-temperature heat exchanger is a revolutionary technology with great application potential. This paper focused on SiC-based cathodic air preheater which is a novel SOFC high temperature heat exchanger, and firstly investigated the coupled radiationconduction-convection heat transfer characteristics between flue gas and air at extremely high temperature conditions. The DO model in ANSYS Fluent was utilized to analyze the radiation heat transfer characteristics of high-temperature flue gas and the effect of gas absorption coefficient, and the simulation results were compared with the S2S model and non-radiation model. The results showed that radiation heat transfer cannot be ignored at high flue gas inlet temperature. With flue gas inlet temperature in the range of 800–1100 °C and low air/flue gas flow rate ratio, the gas radiation heat transfer and the effect of flue gas absorption coefficient can be ignored.

**Keywords:** SOFC · Cathode air preheater · SiC · Heat radiation · Coupled heat transfer characteristics

## **1 Introduction**

In China, Solid oxide fuel cell (SOFC) is favorable in heat power cogeneration and is of great significance in supporting the "dual carbon" target [1]. At present, the fledgling SOFC industry in China mainly focused on the research of stack, however the focus on the BOP (Balance of Plant), which is another major component of SOFC, is not enough. BOP is the crucial technology in SOFC including functions of inlet/outlet gas heat exchange, fuel/air flow distribution, waste heat utilization and so on. It is the key to ensure the efficient and stable operation of the stack.

The operating temperature of SOFC is in the range of 700–1000 °C. It is required that the cathode inlet air temperature of the stack should not be less than 700 °C after stable operation. The BOP recovers the unreacted anode fuel and cathode air of the stack, and sends them to the catalytic combustion chamber for full combustion. The combustion exhaust preheats the cathode air in the high temperature heat exchanger to an appropriate temperature range and then enters the stack to participate in the reaction.

For the high temperature heat exchanger, on the one hand, the increasing stack reaction temperature has the risk of exceeding the temperature tolerance limit of metal heat exchanger represented by superalloy materials. On the other hand, the long-term use of typical commercial metal alloy heat exchanger suffers from the risk of chromium volatilization, which is a known cathodic degradation mechanism, reducing the performance and lifetime of SOFC [2]. Furthermore, the strength of metal alloy heat exchanger decreases rapidly above temperature of 550 °C [3], thus the long-term operation under adverse conditions such as high temperature and oxidation will affect the dynamic service performance of alloy heat exchangers. SiC ceramics and their composite materials are considered as the most promising ceramic materials for high temperature heat exchangers due to their high intensity at high temperature, low density, low thermal expansion coefficient, high thermal conductivity, anti-oxidation, low creep at high temperature and excellent thermal shock resistance [4, 5], thus the application prospects of SiC in SOFC power generation system are promising.

For the high temperature nonmetal heat exchanger, relevant research in China has not been carried out yet, while the international research is just getting started. José Luis Córdova et al. [6, 7] designed and developed a novel cathode air preheater made of alumina ceramics for SOFC. The heat exchanger comprised a group of overlapping quasi-spiral flow channels with rectangular cross section, which was tested to achieve high heat transfer efficiency of 92% and low pressure drop. M. Dev Anand [8] and Umayorupagam P. et al. [9] compared the performance of compact heat exchangers with rectangular channels using SiC and AlN materials. Moreover, researchers have also investigated the application of SiC based high temperature heat exchangers in other fields [10–12]. The research on heat transfer and flow process is relatively simple, and the effect of heat radiation is also ignored in order to simplify heat transfer calculation. However, the temperature of flue gas in SOFC cathode air preheater can be as high as 1100 °C, and the flue gas side contain polar molecules such as CO2 and H2O, thus the flus gas has a certain radiation capacity. Furthermore, SiC ceramic heat transfer element has a dark color and rough surface, thus exhibits high blackness at high surface temperature, which can be approximately considered as blackbody. Therefore, the effect of thermal radiation may not be negligible.

To sum up, replacing the commonly used metal materials with high-temperature resistant SiC ceramic materials for the manufacturing of SOFC high-temperature heat exchanger is a revolutionary technology with great application potential. However, the research on coupled radiation-conduction-convection heat transfer characteristics of SiC high temperature heat exchanger in SOFC needs to be carried out urgently. This paper focused on SiC-based cathodic air preheater which is a novel SOFC high temperature heat exchanger. Based on ANSYS Fluent, the coupled heat transfer characteristics between flue gas and air at extremely high temperature conditions were firstly investigated and the effect of radiation heat transfer were analyzed.

## **2 Structure and Operating Parameters of SOFC Air Preheater**

SOFC cathode air preheater comprises a series of high-temperature flue gas channels and air channels. In order to study the coupled radiation-conduction-convection heat transfer characteristics between flus gas and air, the heat exchanger is simplified as a tubular heat exchanger composed of a hot channel and a cold channel. Figure 1 shows the structure and operating conditions of tubular heat exchanger simulated in this paper. Table 1 shows the simulation parameters of the structure and operating conditions of a natural gas SOFC air preheater. The hydraulic diameter of the internal air channel is 5–10 mm, the internal diameter of the external flue gas channel is 10–20 mm, and the length along the channel is within 500 mm. The heat exchanger material is SiC. The flue gas inlet temperature is in the range of 800–1100 °C, with pressure of 1 bar and flow rate of 1–5 m/s. The flue gas components are H2O (20–30%), CO2 (10–20%), O2 (0–10%) and N2. The inlet air temperature is 25 °C, with pressure of 1 bar and flow rate of 1–5 m/s. The outlet air temperature is required no less than 700 °C.

**Fig. 1.** Structure and operating conditions of simulated heat exchanger

## **3 Numerical Model**

#### **3.1 Radiative Transfer Equation**

Considering the effect of heat radiation, the radiation heat flux must be introduced into the energy equation. In order to calculate the radiation heat flux, the Radiative Transfer Equation (RTE) must be introduced, and the temperature field should be solved by coupling the energy equation and RTE.

The RTE for radiation participating medium at position *r*, in the direction *s* is as follows:

$$\frac{dI(\vec{r},\vec{s})}{ds} + (\alpha + \sigma\_s)I(\vec{r},\vec{s}) = \alpha n^2 \frac{\sigma T^4}{\pi} + \frac{\sigma\_s}{4\pi} \int\_0^{4\pi} I(\vec{r},\vec{s}') \Phi(\vec{s},\vec{s}') d\Omega' \tag{1}$$

where, *s* is the position vector. *r* is the direction vector. *s* is the scattering direction vector. *s* is the path length. α is the absorption coefficient. *n* is the refractive index. σs is the scattering coefficient. σ is the Stefan-Boltzmann constant. *I* is the radiation intensity. *T* is the local temperature. is the phase function. is the solid angle.

**Table 1.** Simulation parameters of the structure and operating conditions.

#### **3.2 Optical Thickness**

Optical thickness (τλ) represents the heat radiation attenuation ability of radiation participating gas at a specific wavelength along given path length. High optical thickness value corresponds to high attenuation effect of the radiation participating gas on the incident radiation energy. τλ is the line integral of the attenuation coefficient along any path in a scattering and absorbing medium which can be expressed as follows:

$$\mathbf{r}\_{\lambda} = \int\_{0}^{L} \kappa\_{\lambda}(l) \mathbf{d}l \tag{2}$$

where, κλ is the spectral attenuation coefficient which depends on the gas components, pressure, temperature and the wavelength of incident radiation energy. κλ is the sum of the spectral absorption coefficient (αλ) and spectral scattering coefficient (σs,λ) at a specific wavelength:

$$
\kappa\_{\lambda} = \alpha\_{\lambda} + \sigma\_{\text{s},\lambda} \tag{3}
$$

For uniform medium, spectral attenuation coefficient is constant along the path length *L*, thus optical thickness has a simple physical interpretation as the length of a path in units of mean free path which can be expressed as follows:

$$
\pi\_{\lambda} = \kappa\_{\lambda} L \tag{4}
$$

In this simulation, the diameter of flue gas channel is in the range of 5–10 mm, and the flue gas contains two polar molecules, CO2 and H2O. The absorption coefficient of CO2 is about 0.5 m−1, and the concentration of CO2 is in the range of 10–20%. The absorption coefficient of H2O is about 0.4 m−1, and the concentration of H2O is in the range of 20–30%. Thus, the optical thickness of flue gas is less than 0.005, which belongs to optical thin gas. In engineering estimation, the absorption effect of CO2 and H2O on incident radiation energy cannot be ignored, but their scattering effect can be ignored. Therefore, to preliminarily analyze the effect of flue gas radiation on heat exchanger performance, the pure absorption and non-scattering condition of gray gas is simulated in this study.

### **3.3 Radiation Heat Transfer Model**

In this paper, DO model was selected to simulate the radiation heat transfer between high-temperature flue gas and air. In addition, S2S model calculation results can be used as comparison to analyze the effect of gas radiation.

#### (1) **DO model**

DO model is applicable to radiation problems in entire range of optical thicknesses. It allows the solution of enclosure radiative transfer without participating media and participating radiation. The DO model solves the RTE of a finite number of discrete solid angles, and the solution precision is controlled by the precision of the discrete solid angles. DO model considers the RTE in the position *s* as a field equation as follows:

$$\nabla(I(\vec{r},\vec{s})\vec{s}) + (\alpha + \sigma\_8)I(\vec{r},\vec{s}) = \alpha n^2 \frac{\sigma T^4}{\pi} + \frac{\sigma\_8}{4\pi} \int\_0^{4\pi} I(\vec{r},\vec{s}') \Phi(\vec{s},\vec{s}') d\Omega' \tag{5}$$

#### (2) **S2S model**

The S2S model is suitable for the solution of enclosure radiative transfer without participating media, however it assumes that all surfaces are diffuse and gray radiation. Thus, S2S model cannot be used to model participating radiation problems. The emissivity (ε) equals the absorptivity (α) and the transmissivity could be neglected. When radiant energy (*E*) is incident on the surface, part (ρ*E*) is reflected, part (α*E*) is absorbed, and part (τ*E*) is transmitted. Considering the conservation of energy, that α + ρ = 1. Since α = ε, thus ρ = 1 − ε.

The radiation heat flux leaving a given surface consists of directly emitted heat flux and reflected heat flux. The reflected heat flux depends on the incident heat flux from the surroundings, which can be expressed by the radiant heat flux leaving all other surfaces. The radiation heat flux leaving from surface *k* can be expressed as follows:

$$
\varepsilon q\_{\rm out,k} = \varepsilon\_k \sigma T^4 + \rho\_k q\_{\rm in,k} \tag{6}
$$

where, *q*out,*k* is the radiation heat flux leaving the surface. ε*k* is the emissivity. ρ*k* is the reflectivity. *q*in,*k* is the radiation heat flux incident on the surface from the surroundings.

The radiant heat flux incident from one surface to another can be calculated by the surface to-surface view factor (*F*jk). Therefore, incident radiant heat flux can be expressed as radiant heat flux leaving all other surfaces as:

$$q\_{\rm in,k} = \sum\_{j=1}^{N} F\_{kj} q\_{\rm out,j} \tag{7}$$

Then, the total radiant heat flux leaving surface *k* can be rewritten as:

$$q\_{\rm out,k} = \varepsilon\_k \sigma T^4 + \rho\_k \sum\_{j=1}^N F\_{kj} q\_{\rm out,j} \tag{8}$$

Equation (8) can be written as follows:

$$J\_k = E\_k + \rho\_k \sum\_{j=1}^{N} F\_{kj} J\_j \tag{9}$$

where, *Jk* represents the heat flux given off from surface *k*, and *Ek* represents the emissive heat flux of surface *k*.

#### **3.4 Grid Independence Study**

A two-dimensional axisymmetric model of coupled radiation-conduction-convection heat transfer between flue gas and water was built according Fig. 1a. The flue gas domain, air domain and channel domain were divided into grids as shown in Fig. 2. The RTE and energy equations are solved based on ANSYS Fluent. A grid independence study was conducted to assess the appropriate number of grid elements. Different grid elements of 2 × 103, 5 × 103, and 1 × 104 were evaluated. Figure 3 showed the comparison of air temperature variation along the channel length calculated with different grid elements. It is found that grid elements of 5 × 103 is adequate for the simulation.

**Fig. 2.** Meshing of simulated heat exchanger

### **4 Results and Discussion**

#### **4.1 Radiation Model Comparison**

In this section, DO model and S2S model in ANSYS Fluent were utilized to analyze the coupled radiation-conduction-convection heat transfer characteristics between hightemperature flue gas and air, and the simulation results were compared with the nonradiation model to analyze the effect of radiation. The detailed simulation parameters were listed in Table 2.

**Fig. 3.** Air temperature variation at different grid elements

**Table 2.** Simulation parameters of different cases.

Figure 4 showed the temperature distributions of coupled heat transfer calculated by the DO radiation model, S2S radiation model and non-radiation model at case 1. It was found that the air outlet temperature was low without considering the radiation while the air outlet temperature is high considering the radiation. This means that the effect of radiation heat transfer cannot be ignored at high flue gas inlet temperature. Comparing the results of different radiation models, it can be seen that the calculation results of DO radiation model and S2S radiation model are very close. Thus, it was preliminarily considered that the effect of gas radiation on coupled heat transfer is tiny, and the effect of gas radiation can be ignored in the coupled heat transfer process, which indicates that the gray gas hypothesis is reasonable.

**Fig. 4.** Comparison of temperature distribution of different radiation models at case 1

Figures 5 and 6 showed the temperature variations of flue gas and air calculated by the DO radiation model, S2S radiation model and non-radiation model at case 1. It was found that the air temperature along the channel increases slightly when considering the radiation. When the flue gas inlet temperature is 800 °C, the air outlet temperature can only be heated to 600 °C. The air temperature along the channel increases greatly when considering the radiation, and the air outlet temperature can be heated to 700 °C. With high flue gas inlet temperature, the contribution of radiation heat transfer to air heating can be as high as 100 °C, thus the influence of radiation heat transfer is very important. Comparing the results of different radiation models, it can be seen that the calculation results of DO radiation model and S2S radiation model are very close, so the influence of gas radiation on heat transfer can be ignored at case 1.

**Fig. 5.** Air temperature variation at case 1

Figures 7 and 8 showed the temperature variations of flue gas and air calculated by the DO radiation model, S2S radiation model and non-radiation model at case 2. The diameter of air channel in case 2 is larger than that in case 1, thus the ratio of air flow rate to flue gas flow rate in case 2 is larger than that in case 1 under the same inlet velocity, which leads to lower air outlet temperature in case 2. Regardless of radiation heat transfer, the air outlet temperature can only be heated to 320 °C. Considering the radiation heat transfer, the air outlet temperature calculated by S2S radiation model is 420 °C, and air outlet temperature calculated by DO radiation model is 400 °C. On one hand, the contribution of radiation heat transfer to air heating at case 2 is up to 100 °C,

**Fig. 6.** Flue gas temperature variation at case 1.

which further proves the importance of radiation heat transfer to coupled heat transfer. On the other hand, comparing the results of different radiation models, it is found that there is a deviation between the DO radiation model and the S2S radiation model results, that is, the gas radiation will affect the proportion of radiation heat transfer at case 2, which indicates that the gas radiation has a certain effect on the coupled heat transfer when the air/gas flow ratio is large. In addition, the air outlet temperature calculated by the S2S radiation model is higher than that calculated by the DO radiation model, which indicates the existence of radiation participating gases slightly reduce the coupled heat transfer characteristics.

**Fig. 7.** Air temperature variation at case 2

Figures 9 and 10 showed the temperature variations of flue gas and air calculated by the DO radiation model, S2S radiation model and non-radiation model at case 3. The diameter of flue gas channel in case 3 is smaller than that in case 2, thus the ratio of air flow rate to flue gas flow rate in case 3 is larger than that in case 1 and case 2 under the same inlet velocity, which leads to lower air outlet temperature in case 3. The air outlet temperature can only be heated to 250 °C without considering radiation heat transfer. Considering the radiation heat radiation heat transfer, air outlet temperature calculated by the S2S radiation model is 290 °C and the DO radiation is 270 °C. On one hand, the contribution of radiation heat transfer to air heating is 20–40 °C in case 3, which means the effect of radiation heat transfer on coupled heat transfer is reduced in case 3. On the

**Fig. 8.** Flue gas temperature variation at case 2.

other hand, comparing the results of different radiation models, it can be seen that the air outlet temperature calculated by S2S radiation model is higher than that calculated by DO radiation model, which further proves that gas radiation has certain effect on the coupled heat transfer when the air/gas flow ratio is large.

**Fig. 9.** Air temperature variation at case 3.

**Fig. 10.** Flue gas temperature variation at case 3.

#### **4.2 Effect of Absorption Coefficient**

According to above analysis, the gas radiation has a certain impact on the coupled heat transfer characteristics only at large air/gas flow ratio. Since the SOFC cathode air preheater in this study is designed to heat the air to above 700 °C using high-temperature flue gas, the air/ gas flow ratio is generally small, thus the effect of gas radiation is unapparent. In addition, the optical thickness of flue gas is very tiny (less than 0.005), thus the gas radiation effect is very tiny. Based on the above analysis, the flue gas radiation has little effect on the coupled heat transfer characteristics.

To further verify the effect of flue gas radiation on the coupled heat transfer characteristics, the gray gas assumption is adopted in this section considering the flue gas absorption coefficient to be independent of wavelength, thus the average absorption coefficient is used in calculation. The DO radiation model is used to analyze the effect of the absorption coefficient on the coupled heat transfer characteristics of flue gas and air with the flue gas absorption coefficient in the range of 0–1 m−1. Figures 11, 12 and 13 showed the air and flue gas temperature variations along the channel with different absorption coefficients at case 1–3. It is found that the absorption coefficient of flue gas in the DO radiation model has little effect on the air and flue gas temperature variation with inlet flue gas temperature of 800 °C at each case. Thus, when inlet flue gas temperature is not very high, the effect of gas radiation can be ignored.

**Fig. 11.** Air and flue gas temperature variation along the channel at case 1 with different absorption coefficients.

Figures 14 and 15 showed the effect of absorption coefficient on air and flue gas temperature variation at case 2 with inlet flue gas temperature of 1100 °C and 1500 °C, respectively. It is found that the absorption coefficient of flue gas has little effect on the air temperature variation with high inlet flue gas temperature, while has small effect on the flue gas temperature variation. Moreover, the higher the inlet flue gas temperature, the greater the effect on the flue gas temperature variation. Since the flue gas temperature in the SOFC cathode air preheater is in the range of 800–1100 °C, the effect of flue gas absorption coefficient and gas radiation can be ignored.

**Fig. 12.** Air and flue gas temperature variation along the channel at case 2 with different absorption coefficients.

**Fig. 13.** Air and flue gas temperature variation along the channel at case 3 with different absorption coefficients

**Fig. 14.** Effect of absorption coefficient on flue gas temperature variation at inlet flue gas temperature of 1100 °C.

### **5 Conclusions**

(1) The air outlet temperature calculated by the non-radiation model is apparently lower than that of DO radiation model and S2S radiation model, thus the effect of radiation heat transfer is very important and cannot be ignored at high flue gas inlet temperature.

**Fig. 15.** Effect of absorption coefficient on flue gas temperature variation at inlet flue gas temperature of 1500 °C.

**Acknowledgements.** This research is sponsored by the Doctor's Science and Technology Fund of HGMRI<PHONE_NUMBER>) and Youth's science and technology Fund of SINOMACH (QNJJ-ZD-2022-06).

## **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Deep Neural Network for Prediction of Adsorbent Selectivity on Hydrogen Purification**

Chenglong Li1,2, Chengsi Xie1, Yi Zong4, Richard Chahine3, Tianqi Yang2(B) , Feng Ye5(B) , and Jinsheng Xiao2,3

1 Automobile Technology and Service College, Wuhan City Polytechnic, Wuhan 430064, Hubei, China

<EMAIL_ADDRESS>

5 School of Artificial Intelligence, Jianghan University, Wuhan 430056, Hubei, China

<EMAIL_ADDRESS>
**Abstract.** With emergence of new materials, more and more materials are available for adsorption and separation processes. The adsorption selectivity of adsorbent to adsorbate is one of the important indicators in choosing materials. Because the adsorption experiment of the mixture is time-consuming and difficult, the selectivity of the adsorbent is generally calculated by the ideal adsorbed solution theory (IAST). Taking the CO2/H2 gas mixture as an example, this paper proposes a new adsorption selectivity calculation method based on a deep neural network (DNN) with 5 hidden layers, which takes the molar fraction of CO2, adsorption pressure and Langmuir adsorption isotherm parameters as the inputs of DNN. Combining the DNN and the NIST/ARPA-E database to quickly and accurately calculate the adsorption selectivity, the hydrogen purification and carbon dioxide storage materials can be quickly screened.

**Keywords:** Hydrogen purification · Ideal adsorbed solution theory · Langmuir isotherm · Selectivity · Deep neural network

## **1 Introduction**

Integrated gasification combined cycle (IGCC) is an advanced power system that combines coal gasification technology and an efficient combined cycle. IGCC will produce a lot of CO2 and H2. It has a good development prospect to capture CO2 efficiently and obtain high-purity H2 by pressure swing adsorption (PSA) technology [1]. The commonly used commercial adsorbents are activated carbon (AC) and zeolite. However, with the emergence of new materials, metal-organic frameworks (MOFs) are increasingly proven to be effective as adsorbents to separate CO2 and H2. Therefore, quickly and accurately selecting the appropriate adsorbent as the separation of CO2 and H2 becomes particularly important.

At present, the selectivity of the adsorbent is one of the important indicators to evaluate the performance of the adsorbent. Adsorption selectivity is the thermodynamic characteristic of the equilibrium-based separations system, which determines the separation efficiency [2]. The method commonly used to solve selectivity is the ideal adsorbed solution theory (IAST). The main practical advantage of IAST is that the adsorption selectivity of adsorbent can be estimated simply from the adsorption isotherm of a single gas, so no special equipment is required for the separation measurement of mixed gas [3]. However, the solution of IAST requires integral calculation from different adsorption isotherm models to solve the spreading pressure. Different adsorption isotherm models will obtain different analytical models, and the solution of selectivity is also relatively complex [4].

Deep neural network (DNN) has good advantages in processing non-linear characteristics and has been increasingly used in the simulation, prediction and optimization of gas adsorption and separation. In this paper, a DNN model is used to predict the separation selectivity. The inputs of DNN are CO2 molar fraction, adsorption pressure and Langmuir adsorption isotherm parameters, and the output of DNN is separation selectivity. The training data of DNN comes from the selectivity by IAST calculations.

### **2 Mathematical Model**

#### **2.1 Langmuir Model for Calculating Adsorption Amount**

Langmuir adsorption model is the most widely used isotherm model. The basic Langmuir adsorption model has simple from:

$$n = \frac{n\_s bp}{1 + bp} \tag{1}$$

where *n* represents the equilibrium adsorption amount, *p* represents the equilibrium adsorption pressure, and *ns* and *b* are saturation capacity and affinity, which are obtained by fitting the experimental values.

The extended Langmuir is used for the adsorption of multi-component gas:

$$m\_l = \frac{n\_{sl} b\_l p\_l}{1 + \sum\_{l} b\_l p\_l} \tag{2}$$

where *i* is the gas component.

#### **2.2 Ideal Adsorbed Solution Theory (IAST) for Calculating Adsorption Amount**

The ideal adsorbed solution theory (IAST) assumes that the adsorbed mixture is an ideal solution at constant spreading pressure and temperature. From the IAST, the spreading pressure π is given by [5],

$$
\pi\_i^0 \left( p\_i^0 \right) = \frac{RT}{A} \int\_0^{p\_i^0} n d \ln p \tag{3}
$$

where a is the adsorbent specific surface area, *R* is the molar gas constant. *p*0 *i* is the gas pressure of component *i* that corresponds to the spreading pressure π*.*

The spreading pressure for a Langmuir isotherm via Eq. (3) is

$$\frac{A}{RT}\pi\_i^0 \left(p\_l^0\right) = n\_s \log\left(1 + bp\_l^0\right) \tag{4}$$

At a constant temperature *T*, the spreading pressure of a single component is the same.

For binary adsorption of component 1 and 2, combining Eqs. (3) and (4),we can get,

$$n\_{s1} \log \left( 1 + b\_1 p\_1^0 \right) = n\_{s2} \log \left( 1 + b\_2 p\_2^0 \right) \tag{5}$$

The IAST requires

$$\mathbf{x}\_{\mathbf{l}}\mathbf{y}\_{\mathbf{l}}p\_{\mathbf{l}} = \mathbf{x}\_{\mathbf{l}}p\_{\mathbf{l}}^{0}, (\mathbf{l} - \mathbf{y}\_{\mathbf{l}})p\_{\mathbf{l}} = (\mathbf{l} - \mathbf{x}\_{\mathbf{l}})p\_{\mathbf{2}}^{0} \tag{6}$$

where *yi* and *xi* are the molar fractions of component *i* in gas phase and adsorbed phase, respectively, and *pt* is the total gas pressure.

#### **2.3 Selectivity Calculated by IAST and Extended Langmuir Equation**

Adsorption selectivity in a binary mixture of components 1 and 2 is defined as

$$S\_{12} = \left(\frac{\chi\_1}{\chi\_1}\right) \Big/ \left(\frac{\chi\_2}{\chi\_2}\right) = \left(\frac{\chi\_1}{\chi\_2}\right) \Big/ \left(\frac{\chi\_2}{\chi\_1}\right) = \left(\frac{n\_1}{n\_2}\right) \Big/ \left(\frac{\chi\_2}{\chi\_1}\right) \tag{7}$$

The adsorption selectivity can be solved by combining Eqs. (5) and (6). If the extended Langmuir equation is used to solve the adsorption selectivity, combining Eqs. (2), (7) becomes

$$LS\_{12} = \left(\frac{n\_{s1}b\_1p\_1}{1+\sum\_{l}b\_lp\_l} / \frac{n\_{s2}b\_2p\_2}{1+\sum\_{l}b\_lp\_l}\right)\left(\frac{y\_2}{y\_1}\right) = \frac{n\_{s1}b\_1}{n\_{s2}b\_2}\tag{8}$$

where *pi* is calculated from the total pressure *pt* by multiplying the mole fraction *yi*. In this paper, component 1 represents CO2, and component 2 represents H2. So here *S*12 represents the selectivity of an adsorbent for CO2 compared to H2. Higher *S*12 is good for hydrogen purification.

#### **2.4 Deep Neural Network (DNN) for Selectivity**

Artificial neural network (ANN) has strong nonlinear processing ability, which includes an input layer, hidden layers and an output layer. A deep neural network (DNN) is a neural network with multiple hidden layers. The more the number of hidden layers, the better the nonlinear fitting characteristics of the ANN, and the better it is to deal with complex problems. In this paper, the number of hidden layers is set to 5, and the number of neurons in each hidden layer is 10. The structure of the DNN in this paper is shown in Fig. 1. The inputs of DNN are CO2 molar fraction, adsorption pressure and Langmuir adsorption isotherm parameters. And the output of DNN is separation selectivity. The chosen method of training data generation is via Latin hypercube sampling (LHS). The significant advantage of LHS is that a small number of samples can represent the entire sample space. The upper and lower boundaries of the LHS are shown in Table 1. The total number of sample points is set to 6000, and the adsorption selectivity values of the selected sample points are calculated by IAST. Due to the large difference between the input parameter values, the calculated selectivity difference is also large. The adsorption selectivity values are set between 3 and 1500. Finally, 3602 sample points are selected to train DNN.

**Fig. 1.** DNN structure for predicting separation selectivity.

**Table 1.** Lower (lb) and upper (ub) boundaries of Latin hypercube sampling (LHS).

### **3 Results and Discussion**

The NIST/ARPA-E database of Novel and Emerging Adsorbent Materials is a free, web-based catalog of adsorbent materials and measured adsorption properties of numerous materials, thus making it a treasure trove for data-driven analysis [6, 7]. As shown in Fig. 2, the Isotherm Visualization Tool of the NIST/ARPA-E database can quickly obtain the adsorption isotherm and Langmuir fitting parameters by searching the adsorbent material and gas species. The adsorption of CO2 and H2 on activated carbon and CuBTC is randomly selected in the NIST/ARPA-E database. Therefore, by combining NIST/ARPA-E Database and DNN model, the selectivity of different adsorbents to CO2 and H2 can be quickly calculated so as to achieve the purpose of fast screening materials.

The algorithm of DNN is completed by Matlab software. The correlation coefficients of the training set, validation set, test set and whole data set of the DNN are 0.99998, 0.9998, 0.99688 and 0.99948, respectively, indicating that the DNN model can well predict the selectivity. In order to further verify the accuracy of the DNN model, the selectivity of activated carbon and CuBTC under different pressures and CO2 molar fractions is calculated by extended Langmuir, IAST and DNN. As shown in Fig. 3, by comparison, the selectivity values calculated by the IAST agree well with the values calculated by the DNN model, which further verifies the robustness of the DNN model. In Fig. 3 (a), for an 80:20 H2/CO2 mixture, CuBTC selectivity under different adsorption pressures is greater than activated carbon. In Fig. 3 (b), at 10 atm, the selectivity of CuBTC is also greater than that of activated carbon at different mole fractions of carbon dioxide. Therefore, for ideal H2 purification and CO2 capture materials, CuBTC is more suitable than activated carbon. Compared with activated carbon, CuBTC has a larger surface area and higher pore volume. There are also related articles about the application of CuBTC in PSA hydrogen purification. For the 81:19 H2/CO2 mixture, the purity 99.99+ of hydrogen can be obtained by a 4-step PSA cycle with CuBTC as an adsorbent [8].