Carbon dioxide free production of hydrogen

The present report summarizes the theoretical modelling and experimental investigation results of the study on the direct thermal methane cracking. This work is a part of the LIMTECH-Project (Liquid Metal Technologies) funded of Helmholtz Alliance and was carried out from 2012 to 2017. The Project-part B5 “CO2-free production of hydrogen” focused on experimental testing and particularly on modelling the novel methane cracking method based on liquid metal technology. The new method uses a bubble column reactor, filled with liquid metal, where both the chemical reaction of methane decomposition and the separation of gas fraction from solid carbon occur. Such reactor system was designed and built in the liquid metal laboratory (KALLA) at KIT. The influences of liquid metal temperature distribution in reactor and feed gas flow rate on methane conversion ratio were investigated experimentally at the temperature range from 930°C to 1175 °C and methane flow rate at the reactor inlet from 50 to 200 mLn/min. In parallel with experimental investigations, a thermochemical model, giving insight in the influence of the above mentioned parameters has been developed at KIT and a CFD model was developed at LUH to get an overview about the bubble dynamics in the reaction system. The influence of different bubble sizes and shapes, multi-inlet coalescence effects as well as the potential of electromagnetic stirring have been investigated.


Introduction
The carbon dioxide emission (CO 2 ) from human activities influences the earth's ecosystem significantly. The worldwide economic growth and development plays a key role in the upward trend in CO 2 emission. As a result, there are many ongoing research works on development of efficient technologies that will be able to reduce the CO 2 emission. One of the technologies, whose usage led to producing a considerable amount of carbon dioxide, is the steam methane reforming (SMR). Depending on the feedstock, this technology provides production of 9 to 12 tons of CO 2 as a byproduct per one ton of produced hydrogen. In the future, current SMR process might be replaced by other environment-friendly hydrogen production technologies based on carbon free sources. However, from the economic point of view, it is not to expect that in the near future hydrogen will be completely produced from renewable and sustainable resources. The proposed direct methane cracking method can in relatively short time provide an attractive industrial solution that leads to the reducing of global carbon dioxide emission while using fossil feedstock, i.e. natural gas.
The chemical reaction of the direct methane cracking can be described by an overall endothermic reaction: The reaction products are gaseous hydrogen and solid elemental carbon. The energy requirement for this chemical reaction is moderate, but a considerable shift of the thermodynamic equilibrium in direction of the reaction products is only achievable beyond 800°C (1 atm). At 1200°C, the theoretical methane conversion efficiency is about 95%.
The direct thermal decomposition of methane is not a new idea, see [2][3][4] for further reference. In the past, a lot of research works were dedicated to the investigation and realizing this process experimentally. The main problem in the practical implementation was the very short operating time of the developed systems. Indeed, the produced carbon is deposited in the reaction zone and this causes the blockage of the reactor. In comparison to previous approaches, the main advantage of the suggested method is the continuous separation and transport of the produced solid carbon from the reaction zone. It occurs owing to a difference of density between liquid metal and carbon and it is supported by the rising gas bubbles. After reaching the liquid metal surface, the carbon powder can be removed from the system with available conventional methods used in the industry (skimming, gas floatation, etc.).
Due to the opaqueness of the liquid metal, visualization of bubble dynamics in the reactor is hardly possible. Numerical simulations have therefore been done in order to investigate dependencies of bubble residence times, bubble shapes, potential carbon sediment locations etc. on operation parameters like gas flow rate. To identify potential optimizations, reactor geometry, inlet types and the heating system have been varied. Especially in this article the aspect of multi-inlets as well as the effect of electromagnetic stirring with regards to increase of hydrogen yield is put in the focus.

Basic process description and experimental set-up
The simplified process flow diagram is shown in Figure 1. The methane cracking reaction occurs in a vertical bubble column reactor filled with liquid metal. Natural gas or pure methane can be used as the feed gas. It is injected at the bottom of the column reactor. The gas bubbles formed at the injector rise to the top of the reactor. Thereby, methane is heated up to liquid metal temperature level and split into gaseous hydrogen and solid carbon. The separation of produced hydrogen from the methane-hydrogen mixture takes place in a separator. The unreacted methane and other possible gaseous by-products can be mixed with the feed gas via mixing unit and moved again to the reactor inlet. Most  cracking reaction collects on the liquid metal free surface and can be removed from the reactor relatively easy. The filter installed within the outlet gas line holds back residual, gas-borne carbon particles.
The construction of experimental test facility HELiS (Hydrogen production Experiment in Liquid Sn) has been part of this project and has been operated at the KArlsruhe Liquid metal LAboratory (KALLA) of KIT. In particular, the facility consists of a gas supply system, an experimental port for the liquid metal reactor, an outlet gas analysis system and an instrumentation and control system. The gas supply system is designed for feeding the working gas in required mixing proportion (methane/argon or nitrogen) into the chemical reactor. It consists of a manifold of standard cylinders with commercial high purity argon at pressure of 200 atm, pipelines, pressure gauges, and gas flow sensors. The use of liquid tin (Sn) as the working fluid in the bubble column reactor has been based upon a careful selection process along several criteria like inert behaviour towards the reaction gases. Major advantages are of course the high thermal conductivity, high density compared to carbon, nontoxicity and long-term chemical stability. The strong corrosion attack of liquid tin at high temperatures towards metallic, particularly steel container materials poses special challenges, which will be detailed in the following.
During the study, two types of the reactor design have been developed and tested: stainless steel reactor and combined quartz glass-steel reactor designs. In general, stainless steel is the favoured construction material due to the great design flexibility and the availability of connectors and flanges even for high temperature applications. However, in contact with molten tin, stainless steel alloying elements are dissolved increasingly with increasing temperature. In the experimental setup, 316 and 1.4549 steels were used as reactor materials. In both cases, experiments showed that materials are corroded quite fast. The lifetime of these reactors was insufficient. In comparison with stainless steel, the quartz glass is chemically stable and showed no corrosion impact towards tin within the planned temperature range. The main idea of this reactor design concept is to couple the advantages of the chemically stable glass material and the mechanical stability of stainless steel. The non-corrosive quartz glass tube getting in contact with liquid metal is inserted into the steel tube. This design proved reliable during the experimental campaigns [2,3,4,5].
In detail, the chemical reactor used in the experiments is a vertical quartz glass tube with inner diameter of 40.6 mm and length of 1268 mm, which is inserted into a stainless steel tube with inner diameter of 49.25 mm and length of 1150 mm (Fig. 2). A single orifice (Ø 0.5 mm), at the bottom of the quartz glass tube was used as gas injector. The combined reactor was placed inside of an 8 kW electric furnace with a maximum operating temperature of 1200 °C. The upper part of the reactor, so called gas zone, is ending above the furnace and was insulated separately. The temperature inside of the reactor was continuously measured during the operating time by applying an alumina tube with type K thermocouples at 10 positions along the whole reactor length. The temperature measured at a level of 600 mm in the middle of the reaction zone was used as reference temperature. Glass fragments or cylindrical rings made of quartz glass have been used to form packed beds. The space porosity of the packed beds was varied from 76% to 84%. The separate supply tank placed above the rector was used for melting of tin prior to filling up the reactor. The gaseous reaction products were filtered via sinter metal elements made of stainless steel with pore size 0.5 µm. The analysis of product gas was realized by using a gas chromatograph (GC). During the experiments, the mole-fraction was continuously measured for methane, hydrogen, nitrogen, ethane, ethylene and acetylene at various operating conditions.

Thermo-chemical modelling and experimental results
The following section describes the thermo-chemical modelling and the experimental results and extends our previous work [3].

Thermo-chemical modelling
To estimate the influence of different parameters on the process as well as for future process description and hydrogen yield predictions, a thermo-chemical (TC) model was developed. The model is based on three coupled one dimensional partial differential equations with spherical coordinates, solved using MATLAB ® parabolic PDE solver. The energy equation (3), the species equation (4) and a pressure equation (5) are solved inside of a bubble traveling through the liquid metal reactor. The bubble is assumed as a spherical body with a constant radius and a constant residence time in the reactor depending on the inventory material and height. As the bubble radius is kept constant, volume expansion due to the pressure drop and the chemical reaction along the liquid metal reactor length is not included in the model. The initial bubble radius in equation (2) is calculated by Tate´s law, which bases on a force balance between buoyancy and surface tension [6] using temperature dependent surface tension and density from Alchagirov et al. [7], [8].
where d orf is the diameter of the orifice,  o and  o are the surface tension and density respectively. In the thermo-chemical model, all gas properties are functions of the present local temperature, pressure and gas composition. The source terms in equation (3) and (4) imply temperature dependent reaction rates, following Arrhenius law and assuming a first order chemical reaction. Energy equation where is the gas density, , is the specific heat capacity of the gas, T is the temperature, r is the bubble radius,  G is the heat conductivity of the gas, 0 is the pre-exponential factor, ∆ is the reaction enthalpy, p is the absolute pressure, is the universal gas constant, is the mole fraction, is the activation energy. Species equation where D is the diffusion coefficient.
where is the residence time.
Assuming linear pressure drop from the bottom to the top of the reactor, whereby the local liquid metal hydrostatic pressure in the column equals the pressure inside of the bubble. The hydrostatic pressure is calculated using the temperature dependent tin density equation from Alchagirov et al. [8].
The feed gas conditions at the bottom of the reactor are used as the initial conditions in equation (3) Symmetry boundary conditions are set at the center of the bubble for the bubble temperature, the species fraction and the pressure (equation (10)). At the bubble interface, symmetry boundary conditions are set for the species fraction and the pressure (equation (11)). As the bubble is rising in the liquid metal due to density differences, and assuming the reference coordinate system fixed to the travelling bubble, the liquid metal is flowing around the bubble. Assuming further, that gas convection inside of the bubble can be neglected and the bubble behaves like a solid sphere, the following Nusselt correlation for forced convection over a sphere from Melissari et. al [9] can be used to model heat transfer from the tin to the bubble. The correlation is valid in a wide range of Pr and Re numbers, including the present Pr number for tin (Pr = 0.0031) at the lowest investigated temperature.
with BB LM ud Re where is the bubble velocity, is the bubble diameter, is the kinematic viscosity of the liquid metal. Even at low Re numbers, resulting from small bubble rise velocities and diameters, the Nu number always leads to a value of 2 or higher. Evaluating the heat transfer limitation, the Bi number in equation (9) reveals a minimum value of 520 for the highest investigated temperature where is the heat conductivity of tin and 4 is the heat conductivity of methane. Consequently the heat transfer limitation is on the gas phase side and the interfacial bubble temperature is set to the local liquid metal temperature and changes with time, and thus with the liquid metal height.
Inside of the bubble ( = 0), the conditions are defined as follows On the gas liquid interface of the bubble ( = ), the conditions are defined as follows To estimate the performance of the process using the measured gas mole fractions from the experiments, the hydrogen yield, the methane conversion and the ratio between the produced hydrogen and the reacted methane are investigated. As the product gas analysis in the GC only considers the components in the gas phase, without taking into account the produced solid carbon, basic calculations to clarify the difference between the measured mole fractions in the GC and the mole fractions of the components, considering carbon production, are necessary. For pure methane as feed gas, assuming the behaviour of an ideal gas, the initial total molar flow rate is given by where ̇4 ,0 is an initial molar flow rate of methane, ̇4 ,0 is an initial volume flow rate of methane. Usually, given by mass conservation, the sum of all components present in the chemical reaction leave the reactor. In this process, the carbon stays inside of reactor, which reduces the total mass and consequently the total molar flow of the product gas. This leads to a new total molar flow entering the GC, changing the mole fraction of each component in the product gas, without taking the carbon mole fraction into account. With the new GC molar flow rate, the calculation of the methane conversion, using GC mole fraction data and Equations (14) Besides the hydrogen yield, the ratio between the produced hydrogen and the reacted methane is mandatory to estimate possible deviations due to the formation of intermediates This ratio leads to a value of 2, in case of all the methane in the feed gas converts to hydrogen and carbon without the formation of intermediates.

Comparison of calculation and measurements
As mentioned above, to quantify the chemical process of the direct methane cracking in liquid metal, several experimental campaigns were performed at the liquid metal laboratory (KALLA). The hydrogen yield measured by the gas chromatograph was used as an indicative parameter to estimate the potential of the process. At the experiments, special focus was given on the influence of the feed gas volume flow rate at different liquid metal temperatures and on the reactor clogging due to carbon formation.
The figure 3 shows the experimental results of methane conversion compared to calculated data as a function of the average liquid metal temperature in reactor at different feed methane volume flow rates. The model prediction showed very good agreement with the experiments.
Especially the liquid metal temperature has a significant influence on the resulting methane conversion, whereby the dependency of the applied methane volume flow rate is quite moderate. The influence of the gas residence time in the liquid metal reactor is only slightly visible by changing the methane feed volume flow rate as the packed bed itself enhances the residence time significantly compared to a reactor design without a packed bed. Nevertheless, applying a pure methane volume flow rate of 200 mL n /min at 1175 °C results in an almost 4 times higher overall hydrogen standard volume flow rate at the outlet compared to a hydrogen output applying a 50 mL n /min initial methane feed volume flow rate.
During the experiments, almost no intermediate reaction products occurred, apart from a maximum mole fraction of 0.5 mol-% ethane and 1.0 mol-% ethylene in the product gas could be detected, whereby no acetylene was measured. As a consequence, the hydrogen selectivity in this process is almost one, leading to a potential product portfolio of only hydrogen and solid carbon.
While disassembling the quartz glass reactor, a thin carbon layer with a thickness in the range of 10 µm was found between the inventory and the quartz glass reactor wall. In the upper part of the reactor, above the metal surface, a mixture of the packed bed and the produced carbon, appearing as a powder, was found (Figure 4). By means of the observed amount, most of the produced carbon accumulated in this area. A smaller amount, transported by the product gas, was collected in the filter elements between the off-gas tube and the GC.

Numerical investigation of fluid dynamic effects within the reactor
To get an idea about the bubble dynamics inside the reactor system of the HELiS facility and to evaluate potential optimization opportunities, additionally to the experiments numerical simulations have been done. According to an efficient production, the dynamics and residence times for different bubble diameters as well as optimization methods like electromagnetic stirring or multi-inlet usages were discovered. Simulations were done using the open source software package OpenFOAM.

Phase description models for numerical bubble simulation
To describe the physics of a two phase bubble flow and to get a numerical description of the dynamics, generally two different methods are available. The discrete phase models for particle tracking based on Euler-Lagrange description and Euler-Euler approaches with a Eulerian description based on an inertial observer formulation for both phases. The Euler-Euler-approaches can be also subdivided into Volume of Fluid, mixture and Euler methods. Prior simulations, especially magnetohydrodynamic modeling have been done using Volume of Fluid method. In the Volume of Fluid method, a volume fraction equation is solved to calculate volume-fraction-averaged parameter that is used to solve the momentum equation. A detailed explanation can be found in [10].
In actual simulations the Eulerian model, described in [11], is used. A solver for coupling this method also with electromagnetic calculations in the software getDP has been developed. One advantage is a more exact calculation for greater bubble volume fractions and greater Stokes numbers, especially at the surface and for different velocities v  of the phases. In contrast to the Volume of Fluid model the Eulerian model solves a separated set of equations for every phase. The coupling between the phases is reached by pressure-and phase-exchange coefficients. The continuity equation, given by is used to solve the volume fraction q  in the momentum equation for every fluid phase: Here ρq describes the density, ∇p the pressure gradient and ṁpq the mass transfer between phase p and q. The q th phase stress tensor τ ̿ q is a function of phase fraction and velocity. External forces ⃗ , lift forces ⃗ , and virtual mass forces ⃗ , can be calculated by correlations described in [12]. The phase exchange coefficient K pq between the dispersed phase q (methane) and continuous phase p (tin) is a function of drag and will be calculated by a Schiller-Naumann-Correlation [13] in the following studies. For the description of energy transfer between the two phases, a separate enthalpy equation has to be solved for each phase: The parameter h q describes the enthalpy, Sq includes possible source terms. In addition, it includes the heat flux Q pq that is a function of heat transfer coefficient between the phases. It is dependent on the Nusselt number. In the following studies the correlation from Ranz and Marshall [14] is used for its description. The fluid dynamics in the system methane and tin are turbulent. To describe the turbulent flow structures in the system a mixture-k-epsilon-model from Behzadi et. al [15] is used. The main idea of this model is based on the assumption that turbulence is dictated by the continuous phase, so turbulent kinetic energy and dissipation for each phase can be linked with constants.

Bubble size studies and slug flow dynamics
In the following studies simulation results of bubble flow in reactor using an Eulerian model will be described. In this section bubble residence times will be analyzed as a function of bubble diameters. As well as the system temperature the residence time is also coupled with the hydrogen yield and offers a potential for optimization. In these models, a tin temperature of 900°C is assumed. Varying the temperature would also influence the rates of ascents. Generally, it has been shown that carbon sediments occur also at walls near the bottom of the reactor, so a greater range from 1 mm bubbles up to bubble sizes near the reactor diameter were discovered (see figure 5.b-d.). The maximum possible bubble size due to reactor diameter is limited to 40 mm. The results of the diameter studies are shownin figure 6. According to experimental results with the system air/water, a spherical bubble shape occurs at smaller bubble diameters, which results in lower rates of ascent and respectively higher residence times. For the spherical shape, rates of ascents increase with bubble diameter until the point of 7 mm. For greater diameters, oscillations occur that results in lower bubble velocities. Further increasing the bubble diameter leads to bubble shape change due to greater Eötvös-number that increases the velocity again. For bubble diameters greater than 20 mm another effect develops. As it can be seen in figure 5.b, stronger wall effects occur, especially at the lower part of the reactor, these effects also lead to bubble dispersion. This is also a possible explanation for carbon sediments occurrence at the reactor walls near the lower part of the reactor and the described velocity reduction. After dispersion effects at the lower part of the reactor, coalescence occur after a short time, so a single bubble leaves the melt after the described residence time (figure 5.c).  For the simulation of constant volume flows used in experiment, the effect of generating greater bubbles and splitting can be seen, as well as the occurrence of wobbling with increasing volume flow, which also fits to other experiments. Another opportunity of comparison, the change of inclination angle of the reactor has been discovered. Increasing the angle, the residence time decreases down to a minimum at 45° and then goes up for greater angles, as it can be seen in experiments for air/water.

Magnetohydrodynamic flow control
To influence bubble residence times, one opportunity is given by the electromagnetic stirring and heating. This can be implemented by an inductor around the bubble column. To calculate the influence of an electromagnetic field on the bubble flow, additional relations, in detail described in [16], are necessary. The Lorentz force that is caused by the induced electromagnetic field of an inductor, has an influence on the velocity field of the melt. This influence and the coupling between hydrodynamic and electromagnetic effects can be described by the following equations: Here ⃗⃗ describes the magnetic field induction, ⃗ is the current density, the magnetic constant, the electrical conductivity of the liquid, the electrical potential. The vector ⃗ describes the volume forces that include gravity g and Lorentz force ⃗ =⃗ x ⃗⃗ . In this case an additional force vector ⃗ appears in the momentum equation described in section 4.1.
The Lorentz forces generate whirl structures in the melt with velocity maxima near the reactor wall (see figure 4.4). As in many melting processes, electromagnetic stirring can be achieved, depending on the inductor frequency and the corresponding skin depth. To optimize the bubble dynamics and the resulting hydrogen yield, the idea is to use the turbulent shear force of the induced melt velocities to split bubbles into blebs. According to result in section 4.2., especially for bubble diameters around 7 mm, a diameter reduction results in higher residence times, besides the surface for heat transfer and reaction increases.  In former studies a volume of fluid method and a k-ω-SST-turbulence model has been implemented in the commercial tool Ansys. To discover the influence of an electromagnetic field on the bubbles, a coil with three windings and a width of 10 mm each is implemented in the lower part of the reactor using In actual studies inductor current and frequencies have been varied. A higher skin depth respectively a decreased the frequency, will increase the stirring effect due to modified distribution of the Lorentz forces. Bubble splitting can be better achieved when decreasing the skin depth by increasing the inductor frequency. But if the current is not changed, a higher power lever would also be needed. Depending on the relation of reactor diameter and skin depth, energy losses reduces with higher frequency which can be used for splitting the bubbles. Generally, it has to be noticed that the change of coupling between the Lorentz forces and the multiphase flow cannot be done with every time step. To identify such phenomena clearly also additional validations using neutron radiography are planned.

Multi-inlet coalescence studies
For a possible scaling up to industrial scale it is interesting to use more than one inlet in the reactor to increase the efficiency. In previous studies, a porous sparger has been used as an alternative inlet. The problem, that is numerical described in [17], is a strong coalescence effects of the resulting bubbles above the porous media. To avoid this, the following passage will include different studies for the relation between the characteristic lengths of the flow problem. This has also been done with regards to a possible scaling up. In these studies, three inlets were discovered. For this case the initial distances between rising bubbles as well as reactor diameter were varied. Decreasing the reactor diameter and the initial distances of the 3 bubbles down to a critical point leads to coalescence (figure 9). Here the reactor volume (dark blue area in figure 9) is reduced. Detailed results are shown in figure 10. Here dR describes the characteristic length of the continuous phase (reactor diameter), while dB can be described as the characteristic length of the dispersed phase (initial bubble diameter). In the diagram the points without bubble coalescence are marked with blue rings. A range for the rates of ascents for the three bubbles is marked. It can be seen that especially for 10  between the bubbles is necessary and for 20 mm also smaller distances avoid coalescence effects. Especially for 5 mm bubbles, an avoidance of coalescence leads to higher residence times. For 10 mm no strong influence can be seen while for 20 mm bubbles residence times even decrease without coalescence. Generally, especially the middle of the three bubbles is influenced by the other ones, and also wobbling effects appear. As seen in section 4.2, an increase of bubble size on smaller bubbles, e.g. by coalescence, has stronger influence on the rates of ascents.
In these studies, the distance between the bubbles are equally separated to the reactor, so the relations dR/dB can be recalculated to initial bubble distances to bubble diameter (aB/dR). So dR/dB = 12 equals aB/dR = 2, dR/dB = 14 equals aB/dR = 2,5 and dR/dB=16 equals aB/dR = 3. For the case dR/dB = 16 (aB/dR = 3) in every discovered cases no coalescence appeared. In other words, the initial bubble distance has to be 3 times higher than initial bubble size to avoid bubble coalescence in the discovered range. It has to be noticed, that the zigzag behavior of the bubbles is increased with greater volume flows, which can be seen in simulations and also experiments using neutron radiography. If it is not possible to generate single bubbles with approximately same diameter, especially at higher constant volume flow, greater distances between the inlets would necessary to avoid the described coalescence effects.

Conclusion
The main objective of the presented work was to better understand the influence of some major influencing parameters on the direct thermal methane cracking within a liquid metal bubble column reactor. The feasibility of the method based on liquid metal technology was successfully demonstrated and presented in [2][3][4].
The results of a mechanistic thermo-chemical 1D-Model have confirmed the strong influence of liquid metal temperature on the methane cracking process. The model predictions are in good agreement with the experimental results. The developed model can be applied to investigate the influence of operating parameters and gas residence time in order to improve the process efficiency.
In the numerical investigation, dynamic studies for different bubble diameters were done. The simulations of volume flows as used in experiments show slug flow dynamics at highest volume flows. Generated bubble sizes vary with volume flow, but can be dispersed while dwelling in the reactor. Generally increasing the volume flow results in stronger bubble wobbling and appearance of strong swirls for higher volume flows. Another aspect, the splitting of methane bubbles into blebs via electromagnetic stirring was also explored. It has been shown that for bubbles with a sufficient high diameter bubbles can be split. Especially for greater bubbles around 9 mm, this split also result in greater residence times, which potentially lead to higher hydrogen yield. Instead of splitting the bubbles, it is generally interesting to avoid their coalescence. Several multi-inlets with ranges of bubble residence times as well as coalescence points were numerically discovered with regards to a possible scaling up of the system. The points of coalescence vary with bubble diameters and the relation between bubble diameter and reactor diameter.
The next steps in technological development will be scaling-up to prototype level and further investigations in view of increasing methane conversion efficiency. The possibility of efficient continuous carbon removal from the reactor may be the key factor for the practical implementation of the developed technology.