Author details

Marvin E. Goldstein1 \* and Pierre Ricco2

1 National Aeronautics and Space Administration, Glenn Research Centre, Cleveland, OH, USA

2 Department of Mechanical Engineering, The University of Sheffield, Sheffield, UK

\*Address all correspondence to<EMAIL_ADDRESS>
© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Leading Edge Receptivity at Subsonic and Moderately Supersonic Mach Numbers DOI: http://dx.doi.org/10.5772/intechopen.83672


*2.1.1 Dynamics of laser map*

Most biological systems exhibit enormous diversity and structurally multicomponent resulting in ecological imbalance and disorder/disharmony in environment. Inspired by articles of Lotka, Volterra, and Allee, numerous articles appeared with diversity in assumptions depending of species and their living envi-

A dynamical system be chaotic then it must be (i) sensitive to initial conditions, (ii) topologically mixing and (iii) its periodic orbits must be dense. In chaotic systems, there exists a strange attractor, a chaotic set, which has fractal structure. Complex systems are also sensitive to their initial conditions and two complex systems that are initially very close together in terms of their various elements and dimensions can end up in distinctly different places. Wide discussions on complex system may found in some pioneer literatures,

similarity that chaotic sets have fractal structure, [60, 68–73].

Chaos measured by Lyapunov exponents, (also called Lyapunov characteristic components or LCEs); LCE > 0 indicates existence of chaos and LCE < 0 indicates regularity, [52–62]. A complex system can better understood by measuring (i) chaos, (ii) Topological entropies and (iii) correlation dimension. Topological entropy, a non-negative number, provides a perfect way to measure complexity of a system. More topological entropy in any system signifies more complexity in it. Actually, it measures the evolution of distinguishable orbits over time, thereby providing an idea of how complex the orbit structure of a system is, [48–50, 61–69]. A system may be chaotic with zero topological entropy. In addition, a significant increase in topological entropy does not justify that it is chaotic. The book by Nagashima and Baba, [62], gives a very clear definition of topological entropy. The correlation dimension provides the dimensionality of the chaotic attractor. Correlation dimensions are non-integers and this is one reasons besides self-

It emerges from a good number of recent researches that chaos appearing in dynamical system be controlled and suggested number techniques to control chaos, [74–88]. These techniques have some limitations depending on the models and

Objective of this article is to investigate the emergence of chaos and complexity

in nonlinear dynamical systems through examples of nonlinear models. Numerical simulations carried out for bifurcation analysis, plotting of LCEs and topological entropies for different systems. Numerical calculations extended to obtain correlation dimensions for certain chaotic attractors emerging in different systems. The study further extended to explain different types of chaos controlling technique. Studies confined to one, two and three-dimensional

Real systems are mostly nonlinear and many of them are with multicomponent structure. Their individual elements possess individual properties. Such systems are

During evolution, a complex system exhibits chaos in some parameter space but also some other phenomena called complexity. This complexity is due to the interaction among multiple agents within the system displayed in the form of coexistence of multiple attractors, bistability, intermittency, cascading effects, exhibit of hysteresis properties etc. Thus, complexity can viewed as its systematic nonlinear properties and it is due to the interaction among multiple agents within the system. Foundation work and elaborate descriptions on complexity can viewed from some pioneer articles on complexity in nonlinear dynamics presented in [45–51]. Study of complexity means to know the results that emerging from a collection of

ronmental conditions in predator-prey models, [20–44].

*A Collection of Papers on Chaos Theory and Its Applications*

termed as the complex system.

interacting parts.

[14, 18, 45, 46, 48, 50, 51].

nature of nonlinearity.

systems only.

**186**

A highly simplified type discrete nonlinear model for laser system, arising from Laser Physics, described in articles, [12, 50, 89–91]. The model describes evolution of certain Fabry-Perot cavity containing a saturable absorber and driven by an external laser represented by

$$\mathbf{x}\_{\mathbf{n}+1} = \mathbf{Q} - \frac{\mathbf{A} \cdot \mathbf{x}\_{\mathbf{n}}}{\mathbf{1} + \mathbf{x}\_{\mathbf{n}}^2}, \forall \mathbf{e} \in \mathcal{R}, \mathbf{n} \in \mathbf{N} \tag{1}$$

Here Q is the normalized input field and A is a parameter depends on the specifics of the parameters and A > 0. The fixed points of the map are the real root of equation

$$\mathbf{x}^3 - \mathbf{Q}\mathbf{x}^2 + (\mathbf{1} + \mathbf{A})\mathbf{x} - \mathbf{Q} = \mathbf{0} \tag{2}$$

This equation has either three real roots or one real and a pair of complex conjugate roots depending on parameter space ð Þ *A*, *Q* . Stability occur in the form of stability and bistability, [89].

#### **Fixed Points and Bifurcations**:

For *Q* fixed, *Q* ¼ 2*:*76, and *A* <4*:*3793, only one stable steady state solution exits and stable two cycle starts when *A* exceeds this value. Thus, approximately, *A* ¼ 4*:*3793, is the bifurcation point. At value *A* ¼ 4*:*3, the stable steady state solution is x\* = 0.720533.

Keeping *Q* ¼ 2*:*76 and varying parameter *A*, bifurcation diagrams are drawn, **Figure 1**, for four different ranges of values of *A*. Similarly, keeping *A* fixed,

**Figure 1.** *Bifurcation diagrams of map (1) for four cases: when Q = 2.76 and parameter A varies.*

*A* ¼ 5*:*4 and varying *Q* in four different ranges, bifurcation diagrams are drawn, **Figure 2**. One observe clearly the appearance of periodic windows within chaotic region of bifurcations as an indication of intermittency and other complex phenomena. Periodic windows become gradually shorter and appearance become more frequent while moving forward in parameter space.

Both time series plots shown in **Figure 3** are for chaotic evolution of system (1) and correspond to parameters (a) ð Þ *A*, *Q* = (5.3, 2.76), due to which an unstable fixed point obtained as *x*� ¼ 0*:*58531, and parameters (b) ) ð Þ *A*, *Q* = (5.4, 2.9), due to which an unstable fixed point obtained as *x*� ¼ 0*:*572218. For both cases, initial point taken is *x*0 ¼ 0*:*5 which lies nearby these points and so, also, unstable.


**1. Introduction**

The European Union aims to reach climate neutrality by the 2050, which means that Europe needs to be carbon neutral years before. Clean energy technology deployment must accelerate rapidly to meet climate goals [1]. Today, hydrogen is mainly produced from fossil fuels, resulting in close to 900 Mt of CO2 emissions per year. Approximately 48% of H2 comes from natural gas through the steam reforming process, 30% comes from naphtha/oil reforming in the chemical industry and 18% comes from coal gasification. Only the remaining 4% is generated by water electrolysis, allowing a CO2-free process, only if the electricity comes from renewable sources [2,3].

Therefore, in order to significantly contribute to the clean energy transition, it is critical to develop low-carbon hydrogen production routes. Natural gas reforming with carbon capture and storage (CCS) (Blue H2) and electrolysis (Green H2) are today the two main alternatives [1]. However, a new path based on methane pyrolysis (Turquoise H2) has garnered increasing interest in recent years and may represent an interesting option during the transition to a long-term sustainable society [4].

Methane cracking is based on the splitting of methane into hydrogen and carbon, without any associated CO2 emissions [5]. The main reaction is endothermic and characterized by strong kinetic limitations due to a high activation energy, which is between 356 and 452 kJ/mol [2].

$$\text{CH}\_4\text{ (g)} = 2\text{H}\_2\text{ (g)} + \text{C}\_{\text{(s)}} \quad \Delta H\_{298K}^0 = -74.8 \text{ kJ/mol} \tag{1}$$

**Citation:** Palo, E.; Cosentino, V.; Iaquaniello, G.; Piemonte, V.; Busillo, E. Thermal Methane Cracking on Molten Metal: Kinetics Modeling for Pilot Reactor Design. *Processes* **2023**, *11*, 1537. https://doi.org/10.3390/ pr11051537

Academic Editors: Elio Santacesaria, Riccardo Tesser and Vincenzo Russo

Received: 30 March 2023 Revised: 8 May 2023 Accepted: 15 May 2023 Published: 17 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The reaction can occur in a thermal, catalytic and combined process. Thermal cracking occurs theoretically above 300 ◦C; however, only around 1000 ◦C can a reasonable high conversion be reached [2]. Figure 1 shows the CH4 conversion calculated at equilibrium at different pressures and temperatures using Aspen Plus software v11 by Aspen Tech, which presents the same trend reported in the literature [6].

**Figure 1.** CH4 equilibrium conversion at different temperatures [◦C] and pressures [bar] calculated on Aspen Plus v11.

Catalyst addition enables a reduction in the activation energy to the range of 96.1–236 kJ/mol, according to the catalyst type. Different kinds of catalysts are studied in the literature for this application. Solid catalysts have been used to lower the activation energy of the reaction and to increase the reaction rate at moderate temperatures (below 1000 ◦C). Such catalysts are either metallic catalysts (Ni, Pd, Pt, Fe, Ni-Pd, etc.) usually supported on metal oxides (Al2O3, MgO, etc.), or carbonaceous catalysts (activated char, carbon black, etc.) [2]. The active metals used in conventional methane cracking can be based on transition metals (Ni > Co > Fe) and noble metals (Pt, Pd), with catalytic activity, however, impacted by coking in time on stream tests, due to the blockage of active sites [7]. Ni-based catalysts are the most active and stable at temperatures between 500 and 700 ◦C, but they rapidly deactivate with increasing temperatures. Moreover, Ni is toxic, therefore any downstream use of a carbon product is enabled if purification is carried out. Fe-based catalysts are cheaper and use a nontoxic metal, but require higher temperatures (700–900 ◦C) [8]. Carbonaceous catalysts, such as activated char/biochar/coal char and carbon black, have also been reported for methane cracking. Their main benefits are a lower cost, no contamination of by-products carbon and more resistance to sulfur [9].

One of the main concerns regarding catalytic cracking is the clogging of the catalyst formed by the carbon, which usually therefore needs to be regenerated in steam or air. Regeneration is not always simple and for some metals, such as Ni, it is almost impossible to recover the initial activity of the catalyst. Furthermore, even if carbonaceous catalysts do not theoretically necessitate regeneration, their catalytic activity also drops after a few hours of operation [2].

To conclude, regardless of whether methane cracking is catalytic or not, the main problem is the cumulative deposition of coke, which can plug not only the catalyst, but the reactor itself after a few hours of operation [10].

A possible solution is to carry out the reaction in a molten media, such as in liquid bubble column reactors [11,12]. The feed gas is then injected at the bottom of the reactor, which rises as bubbles through the molten metal. In this way, the carbon can be easily separated in continuous operation since the low-density solid carbon floats on the surface of the liquid, preventing carbon accumulation and the blockage of the reactor [2,12,13]. This solution still represents a challenge regarding its applicability on an industrial scale [14].

Molten media present other advantages, such as their ability to improve heat transfer and thermal inertia due to their heat capacity, to increase the residence time due to the liquid viscosity, to enable efficient heating using renewable electricity and thus replace fossil fuel combustion, and the liquid phase may catalyze the reaction [2,12].

Molten media can be both molten salts and molten metals. Molten metals show a higher performance in terms of catalytic activity than molten salts. In particular, Ni, Pd, Pt, Co and Fe have a high catalytic behavior, while In, Ga, Sn, Pb and Bi are considered inert metals with low catalytic activity, even if their combination could highly modify the alloy activity and bring unexpected results, exceeding the performance of active metal alloys [2,13]. Furthermore, molten metals should offer isothermal conditions during cracking, and generally have high thermal conductivities. Therefore, they can homogenize the temperature, enhancing methane decomposition [2,12].

Starting from all the above considerations, it was decided that the thermal methane cracking using molten tin in a liquid bubble reactor would be analyzed. The first patent on molten methane cracking was proposed by Tyrer in 1931, although this solution was not considered until the 1990s [15]. Among the studies of the last twenty years, Steinberg was one of the first to propose a bubble reactor configuration with a tin/copper molten metal bath to promote heat transfer and carbon separation [16]. Serban et al. conducted experiments using a stainless-steel feed tube or a porous metal sparger, injecting the methane gas from the top into liquid tin or lead, with or without a packed bed of SiC; they achieved a conversion of 57% at 750 ◦C with a packed bed/tin combination and a porous metal sparger [17]. Geißler et al. studied the thermo-chemical modeling of hydrogen production in a liquid metal bubble column reactor filled with tin, using a packed bed and a methane volume flow rate of 50 mL/min, obtaining a maximum hydrogen yield of 30% at 1000 ◦C [18]. Catalan et al. proposed a model that takes into account the coupling between kinetics and hydrodynamics for thermal and catalytic cracking. The model combines a drift flux correlation for gas hold-up, one of the most critical parameters in this type of application, and a thermodynamically consistent kinetic model of methane decomposition [13]. Subsequently, this first model was improved so that the interfacial area and gas hold-up can be calculated as a function of the operating conditions in the reactor, with no assumptions about the bubble size [19].

The purpose of this work is to analyze in particular the kinetic mechanisms involved in the reaction in the molten media. The main factors to be considered are the high temperature, high residence time in the liquid phase and the diameter of the bubbles along the reactor. A simplified mathematical model to predict methane conversion is developed and solved for both lab-scale and industrial-scale reactors. The results of the first case are compared with experimental data produced by the Department of Chemical Engineering University of La Sapienza in Rome, in order to carry out a first validation of the approach used.

This work is connected to the evaluations of the process scheme and reactor design described in two applied patent applications, one of which will be published by 2023 [12].

#### **2. Materials and Methods**

#### *2.1. Experimental Set-Up*

Experimental tests were carried out at the Industrial Chemistry laboratory, Department of Chemical Engineering University of La Sapienza in Rome.

Three tubular fixed-bed quartz reactors with different diameters were tested: 1.5 cm, 3.3 cm, 6.0 cm ID; in all three cases, a height of up to 20 cm was filled.

All the reactors were heated by an electric furnace, with the temperature varying between 950 ◦C and 1070 ◦C. The reactors were placed vertically inside the furnace. The temperature was controlled using a K-type thermocouple. The selected temperature range complies with the literature: looking at thermodynamics, the reaction is favored at high temperatures; however, almost zero methane conversions are recorded below 900 ◦C [13].

The flow rate of the methane could be varied in a range of 30–140 mL/min. The methane flow rate was controlled by a flowmeter, which was injected into the reactor from a capillary centered at the bottom of the reactor and bubbling inside the molten tin.

A scheme of the experimental set up and the reactor configuration schemes are reported in Figures 2 and 3.

**Figure 2.** Experimental set up scheme.

**Figure 3.** Experimental lab reactor: reactor scheme.

The possibility of cooling the tin surface with Argon was evaluated for the 6 mm ID reactor. In this way, it was possible to slow down the cracking kinetics and avoid a reaction in the gas phase.

The carbon produced was separated by a trap at the reactor exit and the gas was sent to a mass spectrometer analyzer QGA-Hiden.

The experimental CH4 conversion for each reactor diameter is reported in Table 1.

**Table 1.** Experimental data of methane conversion.

#### *2.2. Simplified Kinetics Model*

The methane cracking mechanism in a liquid bubble reactor is determined by both the reaction kinetics and the hydrodynamics of gas rising in the bath. In this study, only the kinetics aspects are analyzed inside the molten media phase by assuming that there are bubbles with a constant diameter along the entire reactor height. Therefore, the mass transfer phenomena effects are not currently considered. This simplified hypothesis will be removed in an updated model based on a dedicated experimental campaign.

In general, in a molten metal, as the methane bubbles rise, two types of reactions occur: one at the center of the bubbles where the methane is not in contact with the liquid, and one at the bubble interface where the gas is in direct contact with the molten media [2]. In this study, it is assumed that the prevalent contribution to methane conversion is due to the heterogeneous surface reaction at the gas/liquid–metal interface; the reaction rate is therefore proportional to the specific contact surface area, which depends on the size distribution of the bubbles [13].

The metal bath, due to its high thermal conductivity, is assumed to have a uniform temperature, independent of the heat required by the reaction [6,20].

The mathematic model developed is based on the methane mass balance inside the molten tin. The methane mass balance is defined in space, along the axial axis of the reactor, and in time. The equations reported from here on were solved using the Gproms software.

The methane mass balance (A = CH4) can be written as follows:

$$V\mathcal{g} \* \frac{d\mathcal{C}a(z,t)}{dt} = \mathcal{Q}\left(z,t\right) \* \mathcal{C}a(z,t)\Big|\_{z} - \mathcal{Q}(z,t) \* \mathcal{C}a(z,t)\Big|\_{z+\Delta Z} - (-ra) \* \mathcal{S} \tag{2}$$

A first-order linear reaction kinetics with respect to methane is considered [21]:

$$(-ra) = k(T) \, \text{Ca}(z, t) \tag{3}$$

where *k*(*T*) is the rate coefficient of the catalytic reaction (Equation (17)) and *Ca*(*z*, *t*) is the concentration of methane as a function of height and time.

The decomposition of methane is a reversible reaction. In the model developed, only the forward reaction is considered, since it can be assumed that the contribution of the backward reaction tends to be zero at the high temperatures analyzed [22,23].

The methane flow rate, as a function of height and time, can be written as follows:

$$Q(z,t) = Q^0(1 + \varepsilon A \* Xa) \tag{4}$$

*Xa* represents methane conversion, defined as follows:

$$Xa = \frac{Q^0 \* \mathbb{C}a^0 - Q(z, t)\mathbb{C}a(z, t)}{Q^0 \* \mathbb{C}a^0} \tag{5}$$

where *Q*0 represents the methane inlet flow rate, *Ca*0 represents the methane initial concentration and *εA* is derived from stoichiometry and defined as follows:

$$
\varepsilon A = \sum \beta i \* y i^0 \tag{6}
$$

where ∑ *βi* represents the sum of the stoichiometric coefficients and *yi*0 represents the molar fraction of methane.

The final balance obtained by combining all the equations above is as follows:

$$\frac{d\mathbb{C}a}{dt} = -\frac{1}{\pi \ast R^2} \left[ \mathbb{Q}\left(z, t\right) \frac{d\mathbb{C}a(z, t)}{dz} + \mathbb{C}a(z, t) \frac{d\mathbb{Q}(z, t)}{dz} \right] - k(T) \ast \mathbb{C}a(z, t) \ast a \tag{7}$$

$$Q(z,t) = Q^0 \left[ 1 + \varepsilon A \* \left( \frac{Q^0 \* \mathbb{C}a^0 - Q\mathbb{C}a}{Q^0 \* \mathbb{C}a^0} \right) \right] \tag{8}$$

with the following initial conditions:

$$\begin{cases} Z = 0, & \mathbb{C}a = \mathbb{C}a^{\circ} \\ Z = 0, & \mathbb{Q} = \mathbb{Q}^{\circ} \end{cases} \tag{9}$$

The specific contact surface area between the methane bubble and the liquid metal, *a*, is an important process parameter, since by increasing its value, a higher methane conversion can be reached. It is calculated as follows [24]:

$$a = \frac{6\varepsilon}{d\_b \* (1 - \varepsilon)}\tag{10}$$

where *db* is the bubble diameter and *ε* is the gas hold up.

As reported in Equation (10), the specific contact surface is also inversely proportional to the diameter of the bubble. This confirms that small bubbles lead to high methane conversions, as also reported in the literature [2].

As explained at the beginning of this section, the bubble diameter is assumed to be constant along the reactor height. The bubble diameter is estimated using Tate's law, assuming that there are no viscous effects, no interactions between the bubbles, that there is a vertical rising pathway, and that there are no variations in the molten surface tension [25]:

$$d\_b = \left(\frac{12R\_0\sigma}{(\rho\_{\rm Tin} - \rho\_\S)g}\right)^{\frac{1}{g}},\tag{11}$$

where

Gas hold-up can be defined as the volume fraction of gas in the total volume of the gas– liquid phase in the bubble column. A smaller gas hold-up requires a larger reactor volume given the same hydrogen production rate, methane conversion and reaction temperature. It can be estimated as the ratio of the gas volume over the total volume of the molten bed.

$$
\varepsilon = \frac{V\_{\text{S}}}{V\_{\text{Sn}}} \tag{12}
$$

where

$$V\_{\mathcal{S}} = Q \ast \mathfrak{r}\_{\mathfrak{m}} \tag{13}$$

$$\pi\_{\rm H} = \frac{H\_{\rm Sn}}{v\_{\rm b}} = \text{rising time} \tag{14}$$

$$V\_{Sn} = \pi \ast R^2 \ast H\_{Sn} \tag{15}$$

$$w\_b = 29.69 \ast d\_b^{0.316} = \text{rising velocity} \left[ \frac{\text{cm}}{\text{s}} \right] \tag{16}$$

with

As the bubble rising velocity is a function of the bubble diameter only, it is assumed constant along the entire reactor height.

The kinetic constant is defined by Arrhenius' law:

$$k(T) = k\_0 \* \exp^{\left(-\frac{\mathbb{k}a}{\mathbb{k}\cdot\mathbb{T}}\right)}\tag{17}$$

Various studies in the literature were analyzed in order to find the values of the preexponential factor *k*0 and the activation energy *Ea* for thermal cracking in molten tin. The values from Rodat et al. [26] were selected since they better fitted the experimental data.

In Table 2, all the parameters fixed to solve the methane balance are reported.

**Table 2.** Parameters fixed to solve methane material balance.

#### **3. Results**

The mathematical model described above is solved using Gprom by entering different tin heights into the software. For each tin height, it is possible to calculate the rising time *τm* (Equation (14)), the volume of the molten bed *VSn* (Equation (15)), the gas volume *Vg* (Equation (13)) and the gas hold up *ε* (Equation (12)).

By combining these parameters to solve the mass balance, it is possible to obtain the methane conversion as a function of the tin height.

To validate this approach, the model results should be compared with the experimental data reported in Table 1. Therefore, the model was solved for 1.5, 3.3 and 6.0 cm ID reactors, considering an initial methane flowrate of 30 mL/min. For all these three cases analyzed, the bubble diameter (Equation (11)) was 2 mm and, accordingly, the rising velocity (Equation (16)) was 178.5 mm/s.

### *3.1. Case A: 1.5 cm and 3.3 cm ID Reactors*

The first two cases analyzed use a 1.5 and 3.3 cm ID reactor. In Figure 4, the conversion profile vs. the tin height is reported. It is evident that the experimental data for the first two reactors are not consistent with the mathematical prediction since the experimental conversions results are much higher than the modelling ones, which tend towards an asymptote for relatively low values.

**Figure 4.** Model and experimental CH4 conversion [%] vs. tin height [m]. (**a**) Case for 1.5 cm ID reactor. (**b**) Case for 3.3 cm ID reactor.

This behavior can be attributed to the small size of the reactors. The wettability of the quartz tube with tin is very low and therefore the bubbles tend to flow along the walls of the tube where the pressure drops are lower. This effect is called the wall effect. Bubble velocity increases on the quartz wall and therefore bubbles leave the tin bed quickly.

This phenomenon is inversely proportional to the diameter of the reactor and is accentuated when the tin height, and therefore the pressure drop, increases. The main effect is that the residence time of methane in the tin bath becomes very limited, leading to a methane conversion that occurs predominantly in the gas phase above the molten metal. The latter must be avoided since the solid carbon, formed outside the bath in the headspace above the molten metal, can deposit in an uncontrolled manner on the reactor surfaces, causing the consequent problems of fouling and/or clogging.

The experimental data in Figure 4 are therefore much higher than the prediction of the model since they consider both the conversion contribution in the liquid phase and in the gaseous phase. As the height of the tin increases, the values diverge more and more, confirming that the conversion occurs mostly in the headspace above the bath.

#### *3.2. Case B: 6.0 cm ID Reactor*

By following the modeling approach already described, it is possible to calculate the methane conversion at different tin heights also for the 6.0 cm ID reactor. The model profile is then compared with the experimental data from Table 1, and the results are shown in Figure 5.

**Figure 5.** Model and experimental CH4 conversion [%] vs. tin height [m] for 6.0 cm ID reactor.

The model conversion still presents an asymptotic behavior. However, unlike the previous cases, the experimental data are more congruent with the model prediction. This result can be attributed on the one hand to an increase in the diameter, which discourages the wall effect with the consequent flow of the bubbles. On the other hand, in the experimental tests of the 6.0 ID reactor, the conversion was measured only in the molten tin, removing the contribution in the gaseous phase. This was possible thanks to a flow of cold Argon, with the same flow rate as the methane feed, which is sent to the surface of the tin bath; in this way, the temperature in the headspace decreases with a slowdown of the kinetics, preventing conversion into the gaseous phase.

The described procedure is only applicable for the 6.0 cm ID reactor due to geometrical constraints.

The main result is that, despite the fact that the mass of the tin is greater than in the previous cases (the volume is greater since at the same relative height, the diameter is bigger), the experimental points in Figure 5 have shifted downwards with respect to those in Figure 4. The values of conversion for the 1.5 cm and 3.3 cm ID reactors should then be corrected by subtracting the conversion in the gaseous phase:

$$X\_{a\\_tim} = X\_{a\\_exp} - X\_{a\\_gas} \tag{18}$$

The conversion into the gaseous phase can be estimated using experimental tests on the empty reactor, which will be the object of a future experimentation campaign.

#### **4. Discussion**

By analyzing the results obtained, it may be concluded that the simplified model describes the experimental data quite accurately. Figure 5 also reflects the methane conversion trend reported in the literature for a similar reactor configuration [2,13]. The curve flattening is due to the small size of the reactor, which leads to a maximum conversion of 20%. This value corresponds to a tin height of about 20 cm, above which the molten metal seems to have no more influence.

It can be concluded that the approach developed is validated.

After this validation, the model is used to predict the methane conversion for an industrial liquid bubble reactor, which is schematized in Figure 6.

**Figure 6.** Liquid bubble reactor scheme for CH4 cracking.

The reactor details are described in the patent applied, which will be published in 2023 [12]. It consists of a metallic reactor internally coated with refractory material and filled with molten tin. The metal is maintained at a constant temperature of 1000 ◦C by electric elements immersed in the molten bath. The methane enters inside the molten tin as bubbles through a gas distributor placed at the bottom of the reactor. The geometry and the dimensions of the reactor have been optimized in the same patents mentioned above.

The main parameters fixed are reported in Table 3.

**Table 3.** Parameters fixed for an industrial reactor.

As already stated, the conversion is strongly influenced by the diameter of the gas bubble since it increases as the bubble diameter decreases. For this reason, three different bubble diameters are analyzed: 0.2 cm, 2.0 cm and 5.0 cm. For each of them, rising velocities (Equation (16)) of 17.85 cm/s, 36.96 cm/s and 49.37 cm/s were calculated, respectively.

Like the lab cases, the model is solved on Gproms, by entering different tin heights into the software: for each of them, the rising time *τm* (Equation (14)), the volume of the molten bed *VSn* (Equation (15)), the gas volume *Vg* (Equation (13)), the gas hold up ε (Equation (12)) and the corresponding methane conversion are calculated.

In Figure 7, the methane conversion profile vs. the tin height is reported for each bubble diameter. Figure 7 shows the equilibrium conversion as well, evaluated in the same operating conditions. In addition, in this case, the conversion presents the same trend as shown in Figures 4 and 5; what changes is the value at which the curve flattens out and the corresponding tin height. The latter is independent of the diameter of the bubble and is equal to 1.5 m in all cases. Instead, the methane conversion, as expected and also reported in the literature [6,26], is strongly influenced by the bubble diameter: with big bubbles, a maximum methane conversion of about 37% is obtained; this value can become double only by decreasing the diameter of the bubbles and keeping all other operating conditions

equal. This demonstrates that the bubble diameter is one of the key factors for this process architecture.

The bubble diameter depends both on the characteristics of the fluid and on the geometry of the sparger. Conversions close to equilibrium can be obtained by suitably designing the sparger to guarantee a bubble diameter in the order of mm or even less.

#### **5. Conclusions**

Methane cracking on molten media represents a promising process for low-emission hydrogen production. It becomes important to understand the thermodynamic and kinetic mechanisms on which the cracking reaction occurs, in order to predict the methane conversion and optimize the reactor design.

For this purpose, a mathematical model was developed in this work. The model is based on the methane mass balance inside the molten tin. The differential equations are solved using Gproms software.

The model results were compared with the experimental data, showing an almost equal trend in the methane conversion during the molten phase.

Afterwards, the mathematical model developed is used to estimate the methane conversion in an industrial reactor at different bubble diameters. The results obtained are promising since quite high conversions can be achieved with small bubbles.

This proves that the bubble diameter is one of the main factors that plays a role in methane cracking. Finding the right compromise between the size of the sparger orifice and the tin bath height is a key parameter in the development of the technology.

**Author Contributions:** Conceptualization, V.C. and E.P.; methodology, V.P., G.I.; validation, E.P., G.I. and V.P.; data curation, E.B.; writing—original draft preparation, V.C.; writing—review and editing, V.C.; supervision, E.P., G.I. and V.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.


On the basis of the three-dimensional model of the pump, a geometric model of the volume of operating fluid filling its interior was prepared (**Figure 1b**). As expected, the result is a very complex structure in terms of geometry. Due to the particular interest in the phenomena occurring in the suction channel of the pump, the calculations used a fragment of the geometric model of the operating fluid volume filling the interior of the pump, which is the volume of oil filling the pump from the inlet to the suction kidneys supplying the fluid to the inter-vane spaces (**Figure 2a**). The separated volume is contained in the immovable elements of the structure, which further simplifies the formulation of the flow problem and the

*(a) 3D model of the vane pump with integrated mechatronic electric drive and (b) 3D model of the operating*

In this chapter, selected examples of numerical calculations will be described, showing the possibility of using CFD methods to solve machine and equipment problems with a rotating operating element, often found in indus-

**146**

**Figure 1.**

*fluid volume filling the pump.*

*(a) Suction channel geometry and (b) discrete model.*

#### **Figure 3.**

*Generic geometry—suction channel pressure distribution for different rotational velocities: (a) 500 rpm, (b) 1000 rpm, (c) 1500 rpm, (d) 2000 rpm, (e) 2500 rpm and (f) 3000 rpm.*

*Generic geometry—suction channel velocity distribution for different velocities: (a) 500 rpm, (b) 1000 rpm, (c) 1500 rpm, (d) 2000 rpm, (e) 2500 rpm and (f) 3000 rpm.*

In order to obtain the most accurate results of the simulation, the "pressure inlet" condition at the inlet and the "mass flow rate" at the outlet were assumed. The mass flow rate was determined using the formula:

$$q\_1 = 2xb\left[\frac{\pi}{2}\left(R\_2^2 - R\_1^2\right) - w\left(R\_2 - R\_1\right)\right] \tag{1}$$

**149**

**Figure 5.**

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans…*

The results of the calculations, apart from pressure distributions, were presented in the form of velocity distributions in the considered area, which are presented in **Figure 4**. From the obtained velocity distributions, it appears that the rotational velocity of the pump significantly influences the velocity of fluid flow in one of the supply channels for both the generic and the modified geometry. It is worth noting that the area where the highest velocities were identified corresponds to the area of the lowest pressures observed in the suction channel. The velocity of the fluid decreases with the lowering of

*Generic geometry—streamlines in the investigated suction channel area for different velocities: (a) 500 rpm,*

*(b) 1000 rpm, (c) 1500 rpm, (d) 2000 rpm, (e) 2500 rpm and (f) 3000 rpm.*

**Figure 5** shows the fluid flow in the form of streamlines, for which the inflow plane to the domain is assigned as the beginning. The results obtained confirm the previous assumptions that the fluid flows evenly and without major turbulences through both inlet channels. Uneven velocity distribution and different pressure values due to asymmetrical layout of channels did not affect the fluid flow. The results obtained on the basis of numerical calculations are the basis for evaluation

Another object under consideration with rotating operating elements was a radial fan. The aim of the numerical simulation was to improve its efficiency. The flow of real gas through a fan with a finite amount of blades is carried out by the cost of loss of energy, called hydraulic losses. Those losses are a consequence of the friction of air molecules occurring on the blade walls and fan housing, vortexes developed in the gas stream, etc. The influence of hydraulic losses on the working characteristic of the radial fan is described by a hydraulic efficiency coefficient, which is defined as the ratio of the useful power to the power delivered by the impeller. This coefficient also defines the real delivery height to the theoretical delivery height—obtained for the finite amount of impeller blades. The impeller geometry considered in possible options, i.e. with eight (factory option) and nine

the rotational velocity, but in the case of simplified geometry, it is slightly lower.

of the structure of channels supplying fluid to the inter-vane volumes.

**3. Radial fan: characteristics and performance improvement**

(suggested option) vanes, are shown in **Figures 6** and **7**.

*DOI: http://dx.doi.org/10.5772/intechopen.82266*

where *q* is the specific mass flow rate; z is the number of vanes; b is the width of a vane; w is the thickness of a vane; *R1* is the small race radius; and R2 is the large race radius.

On this basis, the numerical values entered into the simulation for each impeller velocity were obtained. Within the framework of the study, the analysis of the operating medium flow through the suction channel of the vane pump was performed for various rotational velocity values—changed within the range of 500–3000 rpm.

**Figure 3** shows the pressure distributions in the suction channel of the tested pump for the generic geometry. For each of the cases considered, the lowest pressure occurs in one of the channels supplying fluid to the suction kidneys directly at the inlet to the channel. It was found that the negative pressures for the whole range of rotational velocities are higher than the pressure of oil evaporation, which prevents the occurrence of cavitation phenomena.

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans… DOI: http://dx.doi.org/10.5772/intechopen.82266*

#### **Figure 5.**

*Advanced Computational Fluid Dynamics for Emerging Engineering Processes...*

*Generic geometry—suction channel velocity distribution for different velocities: (a) 500 rpm, (b) 1000 rpm,*

In order to obtain the most accurate results of the simulation, the "pressure inlet" condition at the inlet and the "mass flow rate" at the outlet were assumed.

where *q* is the specific mass flow rate; z is the number of vanes; b is the width of a vane; w is the thickness of a vane; *R1* is the small race radius; and R2 is the large

On this basis, the numerical values entered into the simulation for each impeller velocity were obtained. Within the framework of the study, the analysis of the operating medium flow through the suction channel of the vane pump was performed for various rotational velocity values—changed within the range of 500–3000 rpm. **Figure 3** shows the pressure distributions in the suction channel of the tested pump for the generic geometry. For each of the cases considered, the lowest pressure occurs in one of the channels supplying fluid to the suction kidneys directly at the inlet to the channel. It was found that the negative pressures for the whole range of rotational velocities are higher than the pressure of oil evaporation, which

) − *w*(*R*2 − *R*1)] (1)

\_\_ *π z* (*R*2 2 − *R*1 2

*(c) 1500 rpm, (d) 2000 rpm, (e) 2500 rpm and (f) 3000 rpm.*

*q* = 2*zb*[

The mass flow rate was determined using the formula:

prevents the occurrence of cavitation phenomena.

**148**

**Figure 4.**

race radius.

*Generic geometry—streamlines in the investigated suction channel area for different velocities: (a) 500 rpm, (b) 1000 rpm, (c) 1500 rpm, (d) 2000 rpm, (e) 2500 rpm and (f) 3000 rpm.*

The results of the calculations, apart from pressure distributions, were presented in the form of velocity distributions in the considered area, which are presented in **Figure 4**. From the obtained velocity distributions, it appears that the rotational velocity of the pump significantly influences the velocity of fluid flow in one of the supply channels for both the generic and the modified geometry. It is worth noting that the area where the highest velocities were identified corresponds to the area of the lowest pressures observed in the suction channel. The velocity of the fluid decreases with the lowering of the rotational velocity, but in the case of simplified geometry, it is slightly lower.

**Figure 5** shows the fluid flow in the form of streamlines, for which the inflow plane to the domain is assigned as the beginning. The results obtained confirm the previous assumptions that the fluid flows evenly and without major turbulences through both inlet channels. Uneven velocity distribution and different pressure values due to asymmetrical layout of channels did not affect the fluid flow. The results obtained on the basis of numerical calculations are the basis for evaluation of the structure of channels supplying fluid to the inter-vane volumes.

#### **3. Radial fan: characteristics and performance improvement**

Another object under consideration with rotating operating elements was a radial fan. The aim of the numerical simulation was to improve its efficiency. The flow of real gas through a fan with a finite amount of blades is carried out by the cost of loss of energy, called hydraulic losses. Those losses are a consequence of the friction of air molecules occurring on the blade walls and fan housing, vortexes developed in the gas stream, etc. The influence of hydraulic losses on the working characteristic of the radial fan is described by a hydraulic efficiency coefficient, which is defined as the ratio of the useful power to the power delivered by the impeller. This coefficient also defines the real delivery height to the theoretical delivery height—obtained for the finite amount of impeller blades. The impeller geometry considered in possible options, i.e. with eight (factory option) and nine (suggested option) vanes, are shown in **Figures 6** and **7**.

#### *Advanced Computational Fluid Dynamics for Emerging Engineering Processes...*

**Figure 6.** *Impeller shape: eight vanes.*

#### **Figure 8.**

*Discrete model with division into tetrahedral elements on the impeller and vanes.*

#### **Figure 9.**

*Comparison of the calculation results and the results of the technical documentation for the impeller with the eight and nine vanes.*

**151**

**Figure 10.**

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans…*

For the calculations, the model of impeller according to the enclosed documentation was used as the output model. Calculations have been made for both impeller variants. For both of the cases, the discrete model was based on tetrahedral elements (as exemplary shown in **Figure 8**). Elements near walls were compacted. The flow was modelled as turbulent, using the RANS method and the two-equation turbulence model k-ε. In the first stage of the study, the analysis of the impeller with eight (**Figure 6**)

and nine (**Figure 7**) vanes was carried out. For the eight vanes, the results of the simulation were also compared with the available results in the technical

*Total velocity [m/s] distribution for calculated impeller operating points with eight and nine vanes.*

*DOI: http://dx.doi.org/10.5772/intechopen.82266*

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans… DOI: http://dx.doi.org/10.5772/intechopen.82266*

**Figure 10.** *Total velocity [m/s] distribution for calculated impeller operating points with eight and nine vanes.*

For the calculations, the model of impeller according to the enclosed documentation was used as the output model. Calculations have been made for both impeller variants. For both of the cases, the discrete model was based on tetrahedral elements (as exemplary shown in **Figure 8**). Elements near walls were compacted. The flow was modelled as turbulent, using the RANS method and the two-equation turbulence model k-ε.

In the first stage of the study, the analysis of the impeller with eight (**Figure 6**) and nine (**Figure 7**) vanes was carried out. For the eight vanes, the results of the simulation were also compared with the available results in the technical

*Advanced Computational Fluid Dynamics for Emerging Engineering Processes...*

*Discrete model with division into tetrahedral elements on the impeller and vanes.*

*Comparison of the calculation results and the results of the technical documentation for the impeller with the*

**150**

**Figure 9.**

*eight and nine vanes.*

**Figure 8.**

**Figure 7.**

**Figure 6.**

*Impeller shape: nine vanes.*

*Impeller shape: eight vanes.*

documentation and found to be similar (**Figure 9**). Furthermore, the overall performance of the two types of impellers found with aid of CFD calculation maintains in similar level.

In order to verify the correctness of the calculation of the main dimensions of the impeller, a theoretical design process was carried out. On the basis of known designs, the influence of impeller parameters on its performance, compression and efficiency was simulated. It was necessary to maintain the existing parameters of the impeller, improving only its efficiency. The modifications were limited by the external dimensions of the impeller in order to be able to work with the existing collecting volute.

After a number of variant combinations, the outlet angle of the vane was changed to 23° and the vane profile modified to improve efficiency. The results show that by changing the outlet angle, the average efficiency for the eight-vane impeller was increased by 2.3% and for the nine-vane impeller by 2.9% in relation to the basic eight-vane impeller.

**Figure 10** shows a comparison of the flow images for the impellers with eight and nine vanes with a 23° outlet angle.

The best results were obtained for the nine-vane impeller and the changed outlet angle. An average efficiency increase of 2.9% was achieved in relation to the impeller from the technical documentation. The flow images are correct. There are no particularly dangerous phenomena, such as interruption of flow or turbulence.

#### **4. Centrifugal pump with collecting channel: undetermined flow with cavitation**

Another object of the study was a single-stage centrifugal pump with a spiral volute cooperating with two similar types of impellers, commonly used in such a device. Those impellers are denoted as W13 and W17. The W17 impeller differs from the W13 impeller only by the shape of a vane. Both impellers had eight vanes each. The analysis of the impellers with the two-dimensional peculiarity method for non-viscous medium suggested higher cavitation resistance of the W13 impeller.

In the first stage, calculations were made of the undetermined flow through the pump without cavitation in order to determine the most favourable boundary conditions to be applied when analysing the flow through the pump and determining the calculation characteristics of the pump and the impeller.

The calculations reflect the full three-dimensional geometry of the pump (**Figure 11**) consisting of a straight section of the pipeline before the inlet to the impeller, a centrifugal impeller, a spiral collecting volute, a diffuser, and a short section of pipeline after the pump.

Separate discreet models have been built in the inlet and outlet impeller areas. On the cylindrical surface between the impeller and the volute, these models were not connected by common nodes and remained unfit. Thus, during the calculation it was possible to use the "sliding mesh" technique, which is used to model the rotation of the impeller in relation to the stationary casing. The discrete model is built with approximately 1.3 million tetrahedral elements in total. The elements were also compacted near the vane surface and in the area between the impeller and the collecting channel (**Figure 12**).

The mathematical model of the flow is described by the Reynolds-averaged Navier-Stokes equations (RANS). For the description of the turbulence, a twoequation k-ε model was used. The following control surfaces were used, where static pressure was monitored during the calculation:

**153**

ANSYS Fluent.

**Figure 11.**

**Figure 12.**

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans…*

• The inlet section at the beginning of the suction channel (A-A)

*Discrete model by type of tetrahedral element of the impeller surface, on the hub and rear disc side.*

*The calculation area under consideration and its characteristic cross sections (W13 impeller pump).*

• The cross section at the end of the diffuser (C-C)

Calculations were made according to the scheme:

• The cylindrical surface at the outlet from the impeller inter-vane channel (B-B)

• The outlet section at the end of the cylindrical section of the pipeline (D-D)

In the inlet section (A-A), a homogeneous velocity field was set with the value resulting from the flow rate and the channel section area c = Q/A and the direction corresponding to the connector axis ("velocity inlet" boundary condition). In cross section (D-D), a high static pressure of 1000 kPa was set so that the pressure in the impeller would not drop below the saturation vapour pressure ("pressure outlet" boundary condition). A two-phase flow "mixture" model was selected for the calculations. During the calculations, equations describing the formation of the gaseous phase (cavitation) were excluded. This approach is suggested by

*DOI: http://dx.doi.org/10.5772/intechopen.82266*

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans… DOI: http://dx.doi.org/10.5772/intechopen.82266*

#### **Figure 12.**

*Advanced Computational Fluid Dynamics for Emerging Engineering Processes...*

in similar level.

collecting volute.

**cavitation**

impeller.

the basic eight-vane impeller.

and nine vanes with a 23° outlet angle.

section of pipeline after the pump.

collecting channel (**Figure 12**).

static pressure was monitored during the calculation:

documentation and found to be similar (**Figure 9**). Furthermore, the overall performance of the two types of impellers found with aid of CFD calculation maintains

After a number of variant combinations, the outlet angle of the vane was changed to 23° and the vane profile modified to improve efficiency. The results show that by changing the outlet angle, the average efficiency for the eight-vane impeller was increased by 2.3% and for the nine-vane impeller by 2.9% in relation to

**Figure 10** shows a comparison of the flow images for the impellers with eight

**4. Centrifugal pump with collecting channel: undetermined flow with**

Another object of the study was a single-stage centrifugal pump with a spiral volute cooperating with two similar types of impellers, commonly used in such a device. Those impellers are denoted as W13 and W17. The W17 impeller differs from the W13 impeller only by the shape of a vane. Both impellers had eight vanes each. The analysis of the impellers with the two-dimensional peculiarity method for non-viscous medium suggested higher cavitation resistance of the W13

In the first stage, calculations were made of the undetermined flow through the pump without cavitation in order to determine the most favourable boundary conditions to be applied when analysing the flow through the pump and determin-

The calculations reflect the full three-dimensional geometry of the pump (**Figure 11**) consisting of a straight section of the pipeline before the inlet to the impeller, a centrifugal impeller, a spiral collecting volute, a diffuser, and a short

Separate discreet models have been built in the inlet and outlet impeller areas. On the cylindrical surface between the impeller and the volute, these models were not connected by common nodes and remained unfit. Thus, during the calculation it was possible to use the "sliding mesh" technique, which is used to model the rotation of the impeller in relation to the stationary casing. The discrete model is built with approximately 1.3 million tetrahedral elements in total. The elements were also compacted near the vane surface and in the area between the impeller and the

The mathematical model of the flow is described by the Reynolds-averaged Navier-Stokes equations (RANS). For the description of the turbulence, a twoequation k-ε model was used. The following control surfaces were used, where

ing the calculation characteristics of the pump and the impeller.

The best results were obtained for the nine-vane impeller and the changed outlet angle. An average efficiency increase of 2.9% was achieved in relation to the impeller from the technical documentation. The flow images are correct. There are no particularly dangerous phenomena, such as interruption of flow or turbulence.

In order to verify the correctness of the calculation of the main dimensions of the impeller, a theoretical design process was carried out. On the basis of known designs, the influence of impeller parameters on its performance, compression and efficiency was simulated. It was necessary to maintain the existing parameters of the impeller, improving only its efficiency. The modifications were limited by the external dimensions of the impeller in order to be able to work with the existing

**152**

*Discrete model by type of tetrahedral element of the impeller surface, on the hub and rear disc side.*

Calculations were made according to the scheme:

In the inlet section (A-A), a homogeneous velocity field was set with the value resulting from the flow rate and the channel section area c = Q/A and the direction corresponding to the connector axis ("velocity inlet" boundary condition). In cross section (D-D), a high static pressure of 1000 kPa was set so that the pressure in the impeller would not drop below the saturation vapour pressure ("pressure outlet" boundary condition). A two-phase flow "mixture" model was selected for the calculations. During the calculations, equations describing the formation of the gaseous phase (cavitation) were excluded. This approach is suggested by ANSYS Fluent.

**Figure 13.** *Example pressure pulsation diagram as pressure difference between vane outlet (interface\_2) and inlet (inlet), depending on iteration (time).*

On the internal walls of the flow channel, the condition of zero velocity of the fluid in relation to the wall was set. The increase of static pressure (increase of hydrostatic height) between inlet and outlet cross sections of the pump was the expected value and allowed to reproduce flow characteristics. During the calculations, the average static pressure was monitored on the four control surfaces mentioned above. The calculations were interrupted after repeated oscillations of the static pressure on these surfaces were obtained, which took place after 6–8 rotations of the impeller. An example of a pressure pulsation diagram is shown in **Figure 13**. A fixed time step of Δt = 5,75E–5 s, corresponding to an impeller rotation by 1°, was used for the calculations.

The calculated flow characteristics of the entire pump and the W13 impeller are presented in **Figure 14**. The course of the relevant experimental characteristics is also presented.

The pump characteristics indicate a pressure increase between the cross sections A-A and C-C, characteristics of the impeller—between sections A-A and B-B. The pressure drop in the suction channel is insignificant compared to the pressure drop in the impeller.

#### **Figure 14.**

*Pump and impeller flow characteristics W13 determined by calculation of the transient flow (spiral collecting channel model): comparison with experimental data.*

**155**

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans…*

Cavitation in the pump is associated with a pressure drop in the suction area of the first degree [11]. This causes the fluid-vapour biphasic flow to occur and the continuity of the flow through the pump to be interrupted. In centrifugal pumps, cavitation shall be characterised by a clearly visible disturbance in the following characteristics: flow H = f(Q ), power consumption P = f(Q ) and efficiency η = f(Q ). If the suction height increases at a given velocity and flow rate (or the intake height decreases), then the boundary value of the suction height at which the pump enters the cavitation state is obtained. In this way, taking into account a certain safety margin, it is possible to obtain a curve of the required excess of the energy of a fluid at the pump inlet section over the energy of evaporation of this fluid in the form of NPSH = f(Q ) (net positive suction head). The NPSH parameter

> *ps* − *pv* γ + *cs* 2 \_\_

where *ps* is the absolute pressure at the inlet cross section of the pump and cs is the fluid velocity at the pump inlet cross section (average). Typically, this surplus is

Due to the symmetry of geometry, the flow through the impeller is determined. The elimination of pressure pulsations has significantly accelerated the iterative calculation process. The flow field in the impeller still remained a periodic-symmetric field, but it was the same in all the vane channels. This allowed the calculation area to be limited to one inter-vane channel of the impeller. As a result, the calculation

A discrete model consisting of about 300,000 hexahedral cells was used. Since the discrete model remains stationary during the calculation, a *moving reference*

For the calculations, the "velocity inlet" and "pressure outlet" boundary condi-

Cavitation test in the impeller was performed for several selected values of the flow rate. Calculations were carried out in which equations describing cavitation and two-phase flow were included. The static pressure at the outlet was gradually

It was found that the lowest pressure in the impeller was initially higher than the saturated vapour pressure pmin > pv; then it was already limited by the pv value. For each set outlet pressure, the static inlet pressure was recorded. In the W13 impeller, cavitation occurs on the impeller vanes, close to the incidence edge on the concave side of the vane. In the W17 impeller, cavitation appears on the convex side of

When the outlet pressure is further reduced, it reaches a constant boundary value, depending on the flow rate—fully developed cavitation. Further lowering of the outlet pressure leads to a loss of convergence and interruption of the calculation.

/h). Selected images of the development of cavitation are

tions were used on the outer surface of the annular collecting channel and the two-phase flow "mixture" model. During the calculation, the average static pressure value at the cross sections A-A (inlet) and B-B (outlet from the impeller inter-vane

Determination of the cavitation state in the impeller for a given flow rate requires many calculations of the pressure distribution in the inter-vane space at the decreasing inlet pressure. The simulation assumes that a simplified geometric model of a collective channel can be used to determine the flow characteristics of the impeller itself. Instead of a spiral, an axial-symmetrical guide was used as a

related to a state where the first-stage total head drops by 3% (NPSH3).

time corresponding to one characteristic point has been reduced.

*frame* was used which rotates at the impeller velocity.

2*g* (2)

*DOI: http://dx.doi.org/10.5772/intechopen.82266*

expresses the "suction power" of the pump:

drainage element for the medium.

channel) was monitored.

the vane (for Q = 70 m3

presented in **Figures 15** and **16**.

reduced from 800 to 580 kPa.

*NPSH* = \_\_\_\_\_

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans… DOI: http://dx.doi.org/10.5772/intechopen.82266*

Cavitation in the pump is associated with a pressure drop in the suction area of the first degree [11]. This causes the fluid-vapour biphasic flow to occur and the continuity of the flow through the pump to be interrupted. In centrifugal pumps, cavitation shall be characterised by a clearly visible disturbance in the following characteristics: flow H = f(Q ), power consumption P = f(Q ) and efficiency η = f(Q ). If the suction height increases at a given velocity and flow rate (or the intake height decreases), then the boundary value of the suction height at which the pump enters the cavitation state is obtained. In this way, taking into account a certain safety margin, it is possible to obtain a curve of the required excess of the energy of a fluid at the pump inlet section over the energy of evaporation of this fluid in the form of NPSH = f(Q ) (net positive suction head). The NPSH parameter expresses the "suction power" of the pump:

$$\text{NPSH} \quad = \frac{p\_s - p\_v}{\text{\textdegree\uparrow}} + \frac{c\_s^2}{\text{\textdegree\downarrow}} \tag{2}$$

where *ps* is the absolute pressure at the inlet cross section of the pump and cs is the fluid velocity at the pump inlet cross section (average). Typically, this surplus is related to a state where the first-stage total head drops by 3% (NPSH3).

Determination of the cavitation state in the impeller for a given flow rate requires many calculations of the pressure distribution in the inter-vane space at the decreasing inlet pressure. The simulation assumes that a simplified geometric model of a collective channel can be used to determine the flow characteristics of the impeller itself. Instead of a spiral, an axial-symmetrical guide was used as a drainage element for the medium.

Due to the symmetry of geometry, the flow through the impeller is determined. The elimination of pressure pulsations has significantly accelerated the iterative calculation process. The flow field in the impeller still remained a periodic-symmetric field, but it was the same in all the vane channels. This allowed the calculation area to be limited to one inter-vane channel of the impeller. As a result, the calculation time corresponding to one characteristic point has been reduced.

A discrete model consisting of about 300,000 hexahedral cells was used. Since the discrete model remains stationary during the calculation, a *moving reference frame* was used which rotates at the impeller velocity.

For the calculations, the "velocity inlet" and "pressure outlet" boundary conditions were used on the outer surface of the annular collecting channel and the two-phase flow "mixture" model. During the calculation, the average static pressure value at the cross sections A-A (inlet) and B-B (outlet from the impeller inter-vane channel) was monitored.

Cavitation test in the impeller was performed for several selected values of the flow rate. Calculations were carried out in which equations describing cavitation and two-phase flow were included. The static pressure at the outlet was gradually reduced from 800 to 580 kPa.

It was found that the lowest pressure in the impeller was initially higher than the saturated vapour pressure pmin > pv; then it was already limited by the pv value. For each set outlet pressure, the static inlet pressure was recorded. In the W13 impeller, cavitation occurs on the impeller vanes, close to the incidence edge on the concave side of the vane. In the W17 impeller, cavitation appears on the convex side of the vane (for Q = 70 m3 /h). Selected images of the development of cavitation are presented in **Figures 15** and **16**.

When the outlet pressure is further reduced, it reaches a constant boundary value, depending on the flow rate—fully developed cavitation. Further lowering of the outlet pressure leads to a loss of convergence and interruption of the calculation.

*Advanced Computational Fluid Dynamics for Emerging Engineering Processes...*

*Example pressure pulsation diagram as pressure difference between vane outlet (interface\_2) and inlet (inlet),*

On the internal walls of the flow channel, the condition of zero velocity of the fluid in relation to the wall was set. The increase of static pressure (increase of hydrostatic height) between inlet and outlet cross sections of the pump was the expected value and allowed to reproduce flow characteristics. During the calculations, the average static pressure was monitored on the four control surfaces mentioned above. The calculations were interrupted after repeated oscillations of the static pressure on these surfaces were obtained, which took place after 6–8 rotations of the impeller. An example of a pressure pulsation diagram is shown in **Figure 13**. A fixed time step of Δt = 5,75E–5 s,

The calculated flow characteristics of the entire pump and the W13 impeller are presented in **Figure 14**. The course of the relevant experimental characteristics is

The pump characteristics indicate a pressure increase between the cross sections A-A and C-C, characteristics of the impeller—between sections A-A and B-B. The pressure drop in the suction channel is insignificant compared to the pressure drop in the impeller.

corresponding to an impeller rotation by 1°, was used for the calculations.

*Pump and impeller flow characteristics W13 determined by calculation of the transient flow (spiral collecting*

**154**

**Figure 14.**

*channel model): comparison with experimental data.*

**Figure 13.**

*depending on iteration (time).*

also presented.

*Advanced Computational Fluid Dynamics for Emerging Engineering Processes...*

#### **Figure 15.**

*Cavitation development image on the impeller W13 disc surface at Q = 70 m3 /h and decreasing static pressure at the inlet (percentage of gas phase is given).*

Cavitation image – model with spiral collection channel and immovable impeller (*Moving Reference Frame*) – **Figure 17**. Cavitation image – model with spiral collecting channel and rotating impeller (*Moving Mesh*) – **Figure 18**.

The cavitation fields for the axial-symmetric model are correctly symmetrical. However, the behaviour of the tested impellers is different:

The calculations converge quickly. However, the cavitation fields in the *moving reference frame* model are non-physical, and the cavitation area expands very

#### **Figure 16.**

*Cavitation development image on the impeller W17 disc surface at Q = 70 m3 /h and decreasing static inlet pressure (percentage of gas phase is given).*

**157**

duration.

**Figure 18.**

**Figure 17.**

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans…*

*Cavitation area for parameters: (a) inlet = 189 kPa and outlet = 800 kPa, (b) inlet = 95 kPa and outlet = 700 kPa, (c) inlet = 58 kPa and outlet = 650 kPa, (d) inlet = 27 kPa and outlet = 622 kPa and (e)*

*inlet = 27 kPa and outlet = 600 kPa. The percentage of gas phase is given.*

quickly. Cavitation starts in the direction of the smallest radius of the collecting spiral. The *moving mesh* model produces the best results (mainly physical). However, the problem is the slow convergence of calculations and their long

*Cavitation area for parameters: (a) inlet = 28.3 kPa and outlet = 675 kPa, (b) inlet = 27.8 kPa and outlet = 650 kPa and (c) inlet = 27.8 kPa and outlet = 622 kPa. The percentage of gas phase is given.*

*DOI: http://dx.doi.org/10.5772/intechopen.82266*

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans… DOI: http://dx.doi.org/10.5772/intechopen.82266*

#### **Figure 17.**

*Advanced Computational Fluid Dynamics for Emerging Engineering Processes...*

Cavitation image – model with spiral collection channel and immovable impeller (*Moving Reference Frame*) – **Figure 17**. Cavitation image – model with spiral collect-

The cavitation fields for the axial-symmetric model are correctly symmetrical.

The calculations converge quickly. However, the cavitation fields in the *moving reference frame* model are non-physical, and the cavitation area expands very

• Impeller W13: cavitation is formed on the concave side of the vane.

• Impeller W17: cavitation is formed on the convex side of the vane.

ing channel and rotating impeller (*Moving Mesh*) – **Figure 18**.

However, the behaviour of the tested impellers is different:

*Cavitation development image on the impeller W13 disc surface at Q = 70 m3*

**156**

**Figure 16.**

**Figure 15.**

*the inlet (percentage of gas phase is given).*

*pressure (percentage of gas phase is given).*

*Cavitation development image on the impeller W17 disc surface at Q = 70 m3*

*/h and decreasing static inlet*

*/h and decreasing static pressure at*

*Cavitation area for parameters: (a) inlet = 189 kPa and outlet = 800 kPa, (b) inlet = 95 kPa and outlet = 700 kPa, (c) inlet = 58 kPa and outlet = 650 kPa, (d) inlet = 27 kPa and outlet = 622 kPa and (e) inlet = 27 kPa and outlet = 600 kPa. The percentage of gas phase is given.*

#### **Figure 18.**

*Cavitation area for parameters: (a) inlet = 28.3 kPa and outlet = 675 kPa, (b) inlet = 27.8 kPa and outlet = 650 kPa and (c) inlet = 27.8 kPa and outlet = 622 kPa. The percentage of gas phase is given.*

quickly. Cavitation starts in the direction of the smallest radius of the collecting spiral. The *moving mesh* model produces the best results (mainly physical). However, the problem is the slow convergence of calculations and their long duration.

#### **5. Cavitation resistance of the pump**

The different cavitation properties of the two impellers can be explained by the significantly different inlet angle of the β1 vane—30°40′ (W13) and 21° (W17)—as with the same other geometric data, resulted in a very different position of the ideal inflow point. This is confirmed by the experimental characteristics of the pumps H = f(Q ) and η = f(Q ) from operation.

The analyses indicated the possibility of obtaining information on cavitation resistance of the designed structure through the rational use of CFD programs. The alternative solution of designing a prototype pump and carrying out a series of experiments may be challenging.

#### **6. Conclusion**

The CFD analysis made it possible to identify areas where cavitation is more likely to occur and to assess its intensity in relation to the rotational velocity. The results showed that one of the inlet channels has both negative pressure and increased fluid flow velocity. Calculations made for different pump rotational velocities and different suction channel geometries have shown that the intensity of these phenomena increases with the rotational velocity. However, these phenomena are not strong enough to contribute to the development of the phenomenon of cavitation. A series of simulations for different suction channel geometries have confirmed that no modification of the suction channel geometry is required. Considering the designs presented here the cavitation occurred either on the convex or concave side of the vane. The main difference between the vanes was the angle of its inclination. Hence there is a specific angle between 30 and 21° at which the transition occurs. The volumetric flow rate was unchanged in both of the impeller designs, although the inlet pressures were found to be different. For blades inclined at 30°, the inlet pressures were almost twice lower than in the case of 21°. Hence lower inclination of blades is more immune to cavitation development.

The CFD calculations were made to check the selection of the main dimensions of the radial fan. After performing many variant calculations, it was found that by changing the number of blades and the outlet angle of the blades, it is possible to increase the efficiency of the fan. It appeared that the efficiency is greater for impellers with greater amount of vanes. Furthermore the efficiency increased when the vanes were inclined to 23°, and as stated above, at such angle the cavitation occurs at higher inlet pressures and represent higher immunity to cavitation. Therefore the increase of efficiency may be partially a consequence of lack of cavitation.

It was found that the characteristics of the centrifugal pump from CFD calculations are consistent with the characteristics obtained experimentally. Based on the CFD analysis, cavitation resistance of the designed centrifugal pump was determined.

**159**

**Author details**

Wroclaw, Poland

provided the original work is properly cited.

Wieslaw Fiebig\*, Paulina Szwemin and Maciej Zawislak

\*Address all correspondence to<EMAIL_ADDRESS>
© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

Faculty of Mechanical Engineering, Wroclaw University of Science and Technology,

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans…*

*DOI: http://dx.doi.org/10.5772/intechopen.82266*

*CFD Simulation of Flow Phenomena in Selected Centrifugal Pumps, Industrial Fans… DOI: http://dx.doi.org/10.5772/intechopen.82266*

### **Author details**

*Advanced Computational Fluid Dynamics for Emerging Engineering Processes...*

The different cavitation properties of the two impellers can be explained by the significantly different inlet angle of the β1 vane—30°40′ (W13) and 21° (W17)—as with the same other geometric data, resulted in a very different position of the ideal inflow point. This is confirmed by the experimental characteristics of the pumps

The analyses indicated the possibility of obtaining information on cavitation resistance of the designed structure through the rational use of CFD programs. The alternative solution of designing a prototype pump and carrying out a series of

The CFD analysis made it possible to identify areas where cavitation is more likely to occur and to assess its intensity in relation to the rotational velocity. The results showed that one of the inlet channels has both negative pressure and increased fluid flow velocity. Calculations made for different pump rotational velocities and different suction channel geometries have shown that the intensity of these phenomena increases with the rotational velocity. However, these phenomena are not strong enough to contribute to the development of the phenomenon of cavitation. A series of simulations for different suction channel geometries have confirmed that no modification of the suction channel geometry is required. Considering the designs presented here the cavitation occurred either on the convex or concave side of the vane. The main difference between the vanes was the angle of its inclination. Hence there is a specific angle between 30 and 21° at which the transition occurs. The volumetric flow rate was unchanged in both of the impeller designs, although the inlet pressures were found to be different. For blades inclined at 30°, the inlet pressures were almost twice lower than in the case of 21°. Hence lower inclination of blades is more immune to cavitation

The CFD calculations were made to check the selection of the main dimensions of the radial fan. After performing many variant calculations, it was found that by changing the number of blades and the outlet angle of the blades, it is possible to increase the efficiency of the fan. It appeared that the efficiency is greater for impellers with greater amount of vanes. Furthermore the efficiency increased when the vanes were inclined to 23°, and as stated above, at such angle the cavitation occurs at higher inlet pressures and represent higher immunity to cavitation. Therefore the increase of efficiency may be partially a consequence of lack of

It was found that the characteristics of the centrifugal pump from CFD calculations are consistent with the characteristics obtained experimentally. Based on the CFD analysis, cavitation resistance of the designed centrifugal pump was

**5. Cavitation resistance of the pump**

H = f(Q ) and η = f(Q ) from operation.

experiments may be challenging.

**6. Conclusion**

development.

cavitation.

determined.

**158**

Wieslaw Fiebig\*, Paulina Szwemin and Maciej Zawislak Faculty of Mechanical Engineering, Wroclaw University of Science and Technology, Wroclaw, Poland

\*Address all correspondence to<EMAIL_ADDRESS>
© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **References**

[1] Fiebig W, Cependa P, Jedraszczyk P, Kuczwara H. Innovative solution of an integrated motor pump assembly. In: ASME/BATH 2017 Symposium on Fluid Power and Motion Control. Sarasota, USA: ASME; 2017

[2] Frosina E. A three dimensional cfd modeling methodology applied to improve hydraulic components performance. Energy Procedia. 2015;**82**:950-956

[3] Mancò S, Nervegna N, Rundo M, Armenio G. Modelling and Simulation of Variable Displacement Vane Pumps for IC Engine Lubrication, SAE Technical Paper 2004-01-16012004

[4] Inaguma Y. Theoretical analysis of mechanical efficiency in vane pump. JTEKT Engineering Journal, English Edition No. 1007E Technical paper. 2010:28-35

[5] Frosina E, Stelson KA, et al. Vane pump power split transmission: Three dimensional computational fluid dynamics modeling. In: ASME/BATH 2015 Symposium on Fluid Power & Motion Control. Chicago, USA; 2015

[6] Houzeaux G, Codina R. A finite element method for the solution of rotary pumps. Computers and Fluids. 2007;**36**:667-679

[7] Hyun K, Hazel M, Suresh P. Two-Dimensional CFD Analysis of a Hydraulic Gear Pump. Washington: American Society for Engineeguide Education; 2007

[8] Mochala M. Intermittent CFD simulation of interlocked hydraulic pumps industrial use, basic conditions and prospect. In: FLUIDON Konferenz "Simulation in mechanischen Umfeld", Aachen; 2009

[9] Stryczek J et al. Visualisation research of the flow processes in the outlet chamber—outlet bridge—inlet chamber zone of the gear pumps. Archives of Civil and Mechanical Engineering. 2014;**15**(1):95-108. DOI: 10.1016/j.ACME 2014.02.010

[10] Jedraszczyk P, Fiebig W. CFD Model of an External Gear Pump Proceedings of 13th International Conference Computer Aided. Springer International Publishing; 2017. Available from: https://www.spguideerprofessional. de/en/cfd-model-of-an-external-gearpump/12181546

[11] Tabaczek T, Zawiślak M, Zieliński AK. Flow and cavitation analysis of a centrifugal pump. Systems: Journal of Transdisciplinary Systems Science. 2012;**16**(2):385-394

**160**

*Advanced Computational Fluid Dynamics for Emerging Engineering Processes...*

[9] Stryczek J et al. Visualisation research of the flow processes in the outlet chamber—outlet bridge—inlet chamber zone of the gear pumps. Archives of Civil and Mechanical Engineering. 2014;**15**(1):95-108. DOI:

10.1016/j.ACME 2014.02.010

[10] Jedraszczyk P, Fiebig W. CFD Model of an External Gear Pump Proceedings of 13th International Conference

Computer Aided. Springer International

[11] Tabaczek T, Zawiślak M, Zieliński AK. Flow and cavitation analysis of a centrifugal pump. Systems: Journal of Transdisciplinary Systems Science.

Publishing; 2017. Available from: https://www.spguideerprofessional. de/en/cfd-model-of-an-external-gear-

pump/12181546

2012;**16**(2):385-394

**References**

USA: ASME; 2017

2015;**82**:950-956

2010:28-35

2007;**36**:667-679

Education; 2007

Aachen; 2009

[1] Fiebig W, Cependa P, Jedraszczyk P, Kuczwara H. Innovative solution of an integrated motor pump assembly. In: ASME/BATH 2017 Symposium on Fluid Power and Motion Control. Sarasota,

[2] Frosina E. A three dimensional cfd modeling methodology applied to improve hydraulic components performance. Energy Procedia.

[3] Mancò S, Nervegna N, Rundo M, Armenio G. Modelling and Simulation of Variable Displacement Vane Pumps for IC Engine Lubrication, SAE Technical Paper 2004-01-16012004

[4] Inaguma Y. Theoretical analysis of mechanical efficiency in vane pump. JTEKT Engineering Journal, English Edition No. 1007E Technical paper.

[5] Frosina E, Stelson KA, et al. Vane pump power split transmission: Three dimensional computational fluid dynamics modeling. In: ASME/BATH 2015 Symposium on Fluid Power & Motion Control. Chicago, USA; 2015

[6] Houzeaux G, Codina R. A finite element method for the solution of rotary pumps. Computers and Fluids.

[7] Hyun K, Hazel M, Suresh P. Two-Dimensional CFD Analysis of a Hydraulic Gear Pump. Washington: American Society for Engineeguide

[8] Mochala M. Intermittent CFD simulation of interlocked hydraulic pumps industrial use, basic conditions and prospect. In: FLUIDON Konferenz "Simulation in mechanischen Umfeld",

### *Edited by Albert S. Kim*

As researchers deal with processes and phenomena that are geometrically complex and phenomenologically coupled the demand for high-performance computational fluid dynamics (CFD) increases continuously. The intrinsic nature of coupled irreversibility requires computational tools that can provide physically meaningful results within a reasonable time. This book collects the state-of-the-art CFD research activities and future R&D directions of advanced fluid dynamics. Topics covered include in-depth fundamentals of the Navier–Stokes equation, advanced multi-phase fluid flow, and coupling algorithms of computational fluid and particle dynamics. In the near future, true multi-physics and multi-scale simulation tools must be developed by combining micro-hydrodynamics, fluid dynamics, and chemical reactions within an umbrella of irreversible statistical physics.

Published in London, UK © 2019 IntechOpen © Vanessa Ives / unsplash

Advanced Computational Fluid Dynamics for Emerging Engineering Processes -

Eulerian vs. Lagrangian

Advanced Computational

Fluid Dynamics for Emerging

Engineering Processes

Eulerian vs. Lagrangian

*Edited by Albert S. Kim*


**7.4 Sialic acids**

Most of Sialic acid-specific lectins was found in invertebrates such as those from the Indian horseshoe crab (Mohan et al., 1982), marine crab *Scylla serrata* (Mercy and Ravindranath, 1992), lobster, tunicalase, fungus *Hericium arinaceum* (Kawagishi et al., 1994) and leaves of mulberry (Ratanapo et al., 1998). A lectin from the white shrimp *Litopenaeus setiferus* (LsL) hemolymph is a heterotetramer of two 80 kDa and two 52 kDa subunits, *N*-acetylated sugars, such as GlcNAc, GaINAc, and NeuAc, were the most effective inhibitors of the LsL hemagglutinating activity. Desialylation of erythrocytes or inhibitory glycoproteins abolished their capacity to bind LsL, confirming the relevance of sialic acid in LsL-ligand interactions (Alpuche et al., 2005). In 2009 the *Phaseolus coccineus* lectin (PCL) specificity towards sialic acid showed the molecular mass of 30 kDa consisting of homodimer subunits. Moreover the purified PCL was devoid of antifungal activity against *Candida albicans* and Penicillium *italicum*, but markedly inhibited the growth of *Hericium maydis*, *Rhizoctonia solani*, Gibberella *sanbinetti,* and Sclerotinia *sclerotiorum* while the same concentration of PCL decrease the 50% hemagglutinating activity was inhibited by sialic acid it suggesting a significant correlation between sialic acid-specific site and its bi-functional bioactivities (Chen et al., 2009).


*3.2. CoCVaR Calculation Results*

In addition to numerical results presented in the following sections for new measures, we calculated systemic risk contributions based on the CoCVaR method. We reproduced the case study described in Section 3.3.4 of Huang and Uryasev (2017) with corrected input data (corrected the wrong sign in return data of financial instruments). We considered negative returns (losses) for each bank and the index. We averaged each bank's contributions across time and ranked them accordingly, where larger values correspond to stronger contributions to system's CoCVaR (the units of all reported values are 100%):