Catalyst Distribution Optimization Scheme for Effective Green Hydrogen Production from Biogas Reforming

Abstract

Green hydrogen technology has recently gained in popularity due to the current economic and ecological trends that aim to remove the fossil fuels share in the energy mix. Among various alternatives, biogas reforming is an attractive choice for hydrogen production. To meet the authorities’ requirements, reforming biogas-enriched natural gas and sole biogas is tempting. Highly effective process conditions of biogas reforming are yet to be designed. The current state of the art lacks proper optimization of the process conditions. The optimization should aim to allow for maximization of the process effectiveness and limitation of the phenomena having an adverse influence on the process itself. One of the issues that should be addressed in optimization is the uniformity of temperature inside a reactor. Here we show an optimization design study that aims to unify temperature distribution by novel arrangements of catalysts segments in the model biogas reforming reactor. The acquired numerical results confirm the possibility of the enhancement of reaction effectiveness, coming from improving the thermal conditions. The used amount of catalytic material is remarkably reduced as a side effect of the presented optimization. To ensure an unhindered perception of the reaction improvement, the authors proposed a ratio of the hydrogen output and the amount of used catalyst as a measure.

1. Introduction

The hydrogen industry is currently at the beginning of an undeniable peak. The energy industry has started to pursue technology related to the production and consumption of hydrogen, to meet the ecological regulations introduced by the governmental organizations [1]. The trend gains its strength not only due to a pure interest in the fuel cells technology itself but also because of a rapidly growing necessity of energy storage. The demand for energy storage comes from its periodic oversupply from renewable sources [2]. The leading technologies used for production of hydrogen are hydrocarbons reforming and water electrolysis [3,4]. Water electrolysis from renewables would be the first choice, considering the ecological aspect. However, the technology is still under development and requires further research before its use would be economically encouraged [5]. The present state of the art makes the use of the water electrolysis justified only in applications implementing renewable energy sources, in which no other possibility of the surplus energy consumption is available [6]. Therefore, the majority of the hydrogen produced worldwide comes from the reforming of hydrocarbons [7]. The reforming sovereignty was established due to a relatively low cost of hydrogen production and the reaction placing second, considering the lowest environmental impact [8]. Production of hydrogen on the way of the reforming reaction is assured in the nearest future due to a lack of better alternatives. Furthermore, the size of the natural gas establishment induces no possibility of its immediate shutdown [9,10]. Even though the governmental regulations call for a termination of fossil fuel exploitation, it is impossible to achieve instantly. Thus, natural gas production is inevitable in the nearest future. Although, emissions coming from the energy production using natural gas may be heavily limited, if the natural gas would be converted into hydrogen, on the way of the steam reforming reaction. Considering the current economic sentiments, the reform of hydrocarbons remains the most feasible way of hydrogen production. To decrease the conversion process’ environmental impact, a carbon-capture system can be introduced [11,12]. Furthermore, reform of biogas-enriched natural gas may be a reasonable solution for even better addressing of the ecological issue of the process [13], as biogas is qualified as a renewable energy source [14]. However, biogas production may often require additional energy input, limiting its feasibility. An interesting take on the clean biomass transformation is presented by Gonclaves et al. [15]. Processing of sole biogas, and its integration with a carbon capture system would reduce the reforming’s environmental impact to a minimum [16]. Based on the presented literature review, scientific efforts on the development of the reforming process can not yet be forfeited. Despite its almost century-old history, the reforming technology remains flawed. Carbon deposition and degradation of the catalytic material, due to the thermal stress induced by temperature gradients, are the most remarkable issues [17,18]. Proper maintenance of the thermal conditions during the reaction has been already proven to alleviate the catalyst degradation [19]. Moreover, proper moderation of the temperature field inside the reactor has been reported to increase the overall process’ efficiency [20]. Our team has developed an original strategy to address the issue of temperature gradients occurring in the reactor [21]. The presented article proposes an application of the macro-patterning concept. The concept predicts the division of the reactor’s volume into separate segments, and their further filling with alternating catalytic and non-catalytic materials [21]. The main goal of macro-pattering is to allow for the moderation of the thermal conditions inside the reforming reactor. The strategy is combined with an evolutionary optimization to find the most optimal design of the catalytic insert for the reforming of biogas. The application of macro-patterning is predicted to improve the thermal conditions during the process. The conditions improvement would result in enhancement of the reaction thermodynamics, effecting in reduction of the carbon deposition and premature catalyst degradation [19,22]. The model biogas being the subject of this analysis is assumed to consist only of CH${}_4$ and CO${}_2$. The influence of other biogas components is out of the scope of the presented research [23]. The analysis prepared for the needs of the presented research included:

• Preparation of the numerical simulation considering reforming of model biogas
• Application of the macro-patterning concept for the reactor’s geometry
• Combining the prepared numerical model with a genetic algorithm to find the most optimal catalyst insert design

2. Mathematical Model

2.1. Chemical Reactions

The presented research included formulation of an adequate mathematical model. The analysis of the biogas reforming process can be considered as a quasi three-dimensional, due to the reactor’s geometry and its axial symmetry. Thanks to the axial symmetry, the computational domains may be simplified to two-dimensions, when proper conditions are constituted within the model (Section 3.1) [24]. The geometry implemented into the model of biogas reforming process is presented in Figure 1.

Due to the process character, two main chemical reactions are implemented in the mathematical model [25,26]. According to the previously conducted research, inclusion of the methane-steam reforming (MSR) (Equation (1)), dry reforming (DRY) (Equation (2) and water–gas-shift reaction (WGS) (Equation (3)) is sufficient for acquiring relevant results [27,28]. The reactions are described using the following equations:

$${\mathrm{CH}}_4+{\mathrm{H}}_2\mathrm{O}\to 3{\mathrm{H}}_2+\mathrm{CO},\Delta H_{\mathrm{MSR}}=206.1\frac{\mathrm{kJ}}{\mathrm{mol}},$$

(1)

$${\mathrm{CH}}_4+{\mathrm{CO}}_2\rightleftharpoons 2{\mathrm{H}}_2+2\mathrm{CO},\Delta H_{\mathrm{DRY}}=247\frac{\mathrm{kJ}}{\mathrm{mol}},$$

(2)

$$\mathrm{CO}+{\mathrm{H}}_2\mathrm{O}\rightleftharpoons {\mathrm{H}}_2+{\mathrm{CO}}_2,\Delta H_{\mathrm{WGS}}=-41.15\frac{\mathrm{kJ}}{\mathrm{mol}}.$$

(3)

The enthalpy change $\Delta H$ are taken from literature [21,29]. To allow the inclusion of the reactions into the model, knowledge of their rates is essential. According to the research conducted by Brus et al. [30], the effective rate of MSR and DRY reactions can be expressed with a common equation:

$$R_{\mathrm{eff}}={\dot{w}}_{\mathrm{cat}}{A}_{\mathrm{MSR}}exp\left(-\frac{E_{\mathrm{a}}}{\overline{R}T}\right)p_{{\mathrm{CH}}_4}^{\alpha }{\left(p_{{\mathrm{H}}_2\mathrm{O}}+p_{{\mathrm{CO}}_2}\right)}^{\beta }.$$

(4)

The individual reaction rates for the MSR and DRY reactions can be distinguished as follows:

$$R_{\mathrm{MSR}}=R_{\mathrm{eff}}\frac{p_{{\mathrm{H}}_2\mathrm{O}}}{p_{{\mathrm{CO}}_2}+p_{{\mathrm{H}}_2\mathrm{O}}},$$

(5)

$$R_{\mathrm{DRY}}=R_{\mathrm{eff}}\frac{p_{{\mathrm{CO}}_2}}{p_{{\mathrm{CO}}_2}+p_{{\mathrm{H}}_2\mathrm{O}}}.$$

(6)

The WGS reaction has a more unpredictable nature. Thus, preparation of a formula, returning proper values regardless of the process conditions, is not plausible. However, according to Ahmed and Föger, the WGS reaction can be assumed to maintain equilibrium under specific conditions [31]. The correctness of the described approach is confirmed by other research [32,33,34]. The boundary conditions applied in the presented research, described in Section 3.1, are designed to satisfy the equilibrium assumption. Therefore, the WGS rate is calculated using the following procedure:

$$K_{\mathrm{WGS}}={\displaystyle \frac{k_{\mathrm{WGS}}^{+}}{k_{\mathrm{WGS}}^{-}}}={\displaystyle \frac{p_{{\mathrm{CO}}_2}{p}_{{\mathrm{H}}_2}}{p_{\mathrm{CO}}{p}_{{\mathrm{H}}_2\mathrm{O}}}}=exp\left(-{\displaystyle \frac{\Delta G_{\mathrm{WGS}}^0}{\overline{R}T}}\right).$$

(7)

Combining of the Equations (4) and (7) allows for formulation of the WGS reaction rate equation, as follows:

$$R_{\mathrm{WGS}}=k_{\mathrm{WGS}}^{+}{p}_{\mathrm{CO}}{p}_{{\mathrm{H}}_2\mathrm{O}}+k_{\mathrm{WGS}}^{-}{p}_{{\mathrm{H}}_2}{p}_{{\mathrm{CO}}_2}.$$

(8)

The values of partial pressures, included in the Equations (4)–(8), are derived from an analysis of the reactions’ stoichometry [30].

2.2. Heat and Mass Transfer

The prepared mathematical model has the fundamental transport equations incorporated. Considering the defined computational domain, the equations are formulated for a two-dimensional case. Therefore, the model’s equations are calculated for the longitudinal and radial directions only. The fluids taking part in the reaction are considered Newtonian and incompressible. Thus, the continuity equation (Equation (9)) takes the following form [35]:

$$\frac{\partial \left({\rho }_0{U}_x\right)}{\partial x}+\frac{1}{r}\frac{\partial \left(r{\rho }_0{U}_r\right)}{\partial r}=0.$$

(9)

The materials serving as the catalytic insert are porous. Therefore, the terms of the momentum conservation equation have to contain parameters depending on the structure morphology [36]. A separate equation for each dimension is formulated (Equations (10) and (11)):

$$\begin{array}{c}\hfill \begin{array}{c}\frac{{\rho }_0}{{\epsilon }_0^2}\left(U_x\frac{\partial U_x}{\partial x}+U_r\frac{\partial U_x}{\partial r}\right)=\hfill \\ \hfill -\frac{\partial P}{\partial x}+\frac{\mu }{{\epsilon }_0}\left[\frac{{\partial }^2{U}_x}{\partial x^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial U_x}{\partial r}\right)\right]\\ \hfill -\frac{\mu }{K_{\mathrm{p}}}{U}_x-\frac{{\rho }_0{c}_{\mathrm{ine}}}{\sqrt{K_{\mathrm{p}}}}{U}_x\sqrt{U_x^2+U_r^2},\end{array}\end{array}$$

(10)

$$\begin{array}{c}\hfill \begin{array}{c}\frac{{\rho }_0}{{\epsilon }_0^2}\left(U_x\frac{\partial U_r}{\partial x}+U_r\frac{\partial U_r}{\partial r}\right)=\hfill \\ \hfill -\frac{\partial P}{\partial r}+\frac{\mu }{{\epsilon }_0}\left[\frac{{\partial }^2{U}_r}{\partial x^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial U_r}{\partial r}\right)-\frac{U_r}{r^2}\right]\\ \hfill -\frac{\mu }{K_{\mathrm{p}}}{U}_r-\frac{{\rho }_0{c}_{\mathrm{ine}}}{\sqrt{K_{\mathrm{p}}}}{U}_r\sqrt{U_x^2+U_r^2},\end{array}\end{array}$$

(11)

The use of metallic foams implies a necessity of application of a relevant model, allowing us to calculate the effective thermal conductivity coefficient of the material ${\lambda }_{\mathrm{eff}}$ [37,38]. Knowledge of a proper ${\lambda }_{\mathrm{eff}}$ value is vital for the calculation of the energy conservation equation (Equation (12)) [36]. An adequate model is prepared, described in detail within our previous article [39].

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\rho }_0{C}_p\left(U_x\frac{\partial T_{\mathrm{loc}}}{\partial x}+U_r\frac{\partial T_{\mathrm{loc}}}{\partial r}\right)=\frac{\partial }{\partial x}\left({\lambda }_{\mathrm{eff}}\frac{\partial T_{\mathrm{loc}}}{\partial x}\right)\\ \hfill +\frac{1}{r}\frac{\partial }{\partial r}\left(r{\lambda }_{\mathrm{eff}}\frac{\partial T_{\mathrm{loc}}}{\partial r}\right)+Q_{\mathrm{s}},\end{array}\end{array}$$

(12)

where the heat sources/sinks $Q_{\mathrm{s}}$ depend on enthalpy changes $\Delta H$ and rates of the MSR, DRY and WGS [25,30], and are defined using following equations:

$$Q_{\mathrm{MSR}}=-\Delta H_{\mathrm{MSR}}{R}_{\mathrm{MSR}},$$

(13)

$$Q_{\mathrm{DRY}}=-\Delta H_{\mathrm{DRY}}{R}_{\mathrm{DRY}},$$

(14)

$$Q_{\mathrm{WGS}}=-\Delta H_{\mathrm{WGS}}{R}_{\mathrm{WGS}}.$$

(15)

The species conservation is calculated using molar fractions of species taking part in the reaction (Equation (16)). The formulated equation was derived from the Fick’s law of diffusion [40]:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\rho }_0\left(U_x\frac{\partial Y_j}{\partial x}+U_r\frac{\partial Y_j}{\partial r}\right)=\frac{\partial }{\partial x}\left({\rho }_0{D}_{j,\mathrm{eff}}\frac{\partial Y_j}{\partial x}\right)\\ \hfill +\frac{1}{r}\frac{\partial }{\partial r}\left(r{\rho }_0{D}_{j,\mathrm{eff}}\frac{\partial Y_j}{\partial r}\right)+S_j.\end{array}\end{array}$$

(16)

The mass sources and sinks $S_j$ depend on the MSR, DRY and WGS rates and molar masses of the species taking part in the reaction [41,42]. The exact equations defining the values of $S_j$ are described in the

Table 1. Mass generation.

Species	Mass Generation MSR	Mass Generation WGS	Mass Generation DRY	Summarized Generation
${\mathrm{H}}_2$	$3R_{\mathrm{MSR}}{M}_{{\mathrm{H}}_2}$	$R_{\mathrm{WGS}}{M}_{{\mathrm{H}}_2}$	$2R_{\mathrm{DRY}}{M}_{{\mathrm{H}}_2}$	$3R_{\mathrm{MSR}}{M}_{{\mathrm{H}}_2}+R_{\mathrm{WGS}}{M}_{{\mathrm{H}}_2}$ + $2R_{\mathrm{DRY}}{M}_{{\mathrm{H}}_2}$
CO	$R_{\mathrm{MSR}}{M}_{\mathrm{CO}}$	$-R_{\mathrm{WGS}}{M}_{\mathrm{CO}}$	$2R_{\mathrm{DRY}}{M}_{\mathrm{CO}}$	$R_{\mathrm{MSR}}{M}_{\mathrm{CO}}-R_{\mathrm{WGS}}{M}_{\mathrm{CO}}$ + $2R_{\mathrm{DRY}}{M}_{\mathrm{CO}}$
${\mathrm{CO}}_2$	0	$R_{\mathrm{WGS}}{M}_{{\mathrm{CO}}_2}$	$-R_{\mathrm{DRY}}{M}_{{\mathrm{CO}}_2}$	$R_{\mathrm{WGS}}{M}_{{\mathrm{CO}}_2}$$-2R_{\mathrm{DRY}}{M}_{{\mathrm{H}}_2}$
${\mathrm{CH}}_4$	$-R_{\mathrm{MSR}}{M}_{{\mathrm{CH}}_4}$	0	$-R_{\mathrm{DRY}}{M}_{{\mathrm{H}}_2}$	$-R_{\mathrm{MSR}}{M}_{{\mathrm{CH}}_4}$$-R_{\mathrm{DRY}}{M}_{{\mathrm{H}}_2}$
${\mathrm{H}}_2\mathrm{O}$	$-R_{\mathrm{MSR}}{M}_{{\mathrm{H}}_2\mathrm{O}}$	$-R_{\mathrm{WGS}}{M}_{{\mathrm{H}}_2\mathrm{O}}$	0	$-R_{\mathrm{MSR}}{M}_{{\mathrm{H}}_2\mathrm{O}}$$-R_{\mathrm{WGS}}{M}_{{\mathrm{H}}_2\mathrm{O}}$

.

3. Numerical Model

Proceeding with computations of the prepared mathematical model, requires its implementation within a relevant numerical procedure. To allow doing so, a discretization of the computational domain is necessary. The governing equations are applied into the numerical model using the Finite Volume Method [24,43]. The generalized transport equation, representing the partial differential Equations (9)–(12) and (16), takes the following form [24,44]:

$${\Psi }_x{\displaystyle \frac{\partial \varphi }{\partial x}}+{\Psi }_r{\displaystyle \frac{\partial \varphi }{\partial r}}={\displaystyle \frac{\partial }{\partial x}}\left(\Gamma {\displaystyle \frac{\partial \varphi }{\partial x}}\right)+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r\Gamma {\displaystyle \frac{\partial \varphi }{\partial r}}\right)+\overline{S}.$$

(17)

The source terms $\overline{S}$ in Equation (17) are gathered in the

Table 2. Source terms in Equation (17).

Equation	$\overline{\mathit{S}}$
(10)	$-{\displaystyle \frac{\partial P}{\partial x}}-{\displaystyle \frac{\mu }{K_{\mathrm{p}}}}{U}_x-{\displaystyle \frac{{\rho }_0{c}_{\mathrm{ine}}}{\sqrt{K_{\mathrm{p}}}}}{U}_x\sqrt{U_x^2+U_r^2}$
(11)	$-{\displaystyle \frac{\partial P}{\partial r}}-{\displaystyle \frac{\mu }{K_{\mathrm{p}}}}{U}_r-{\displaystyle \frac{{\rho }_0{c}_{\mathrm{ine}}}{\sqrt{K_{\mathrm{p}}}}}{U}_r\sqrt{U_x^2+U_r^2}-{\displaystyle \frac{\mu U_r}{\epsilon r^2}}$
(12)	$Q_{\mathrm{s}}$
(16)	$S_j$

. The values presented in the Table 2 are relevant for segments filled with the catalytic material. When a non-catalytic segment is considered, the source terms values are set to 0. The chemical reactions are assumed to be suppressed, just after the gases mixture enters the non-catalytic region of the reformer [21]. To acquire the discrete transport equation (Equation (18)), integration of the Equation (17) is conducted. After the application of simple, mathematical transformations, the discrete transport equation is expressed as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\left[{\left({\Psi }_x\varphi \right)}_{\mathrm{e}}-{\left({\Psi }_x\varphi \right)}_{\mathrm{w}}\right]r_{\mathrm{m}}\Delta r+\left[{\left(r{\Psi }_r\varphi \right)}_{\mathrm{n}}-{\left(r{\Psi }_r\varphi \right)}_{\mathrm{s}}\right]\Delta x=\hfill \\ \hfill \left({\Gamma }_{\mathrm{e}}\frac{{\varphi }_{\mathrm{E}}-{\varphi }_{\mathrm{P}}}{\delta x}-{\Gamma }_{\mathrm{w}}\frac{{\varphi }_{\mathrm{P}}-{\varphi }_{\mathrm{W}}}{\delta x}\right)r_{\mathrm{m}}\Delta r\hfill \\ \hfill +\left(r{\Gamma }_{\mathrm{n}}\frac{{\varphi }_{\mathrm{N}}-{\varphi }_{\mathrm{P}}}{\delta r}-r{\Gamma }_{\mathrm{s}}\frac{{\varphi }_{\mathrm{P}}-{\varphi }_{\mathrm{S}}}{\delta r}\right)\Delta x+\overline{S}{r}_{\mathrm{m}}\Delta r\Delta x,\end{array}\end{array}$$

(18)

$$a_{\mathrm{P}}{\varphi }_{\mathrm{P}}=a_{\mathrm{E}}{\varphi }_{\mathrm{E}}+a_{\mathrm{W}}{\varphi }_{\mathrm{W}}+a_{\mathrm{N}}{\varphi }_{\mathrm{N}}+a_{\mathrm{S}}{\varphi }_{\mathrm{S}}+b,$$

(19)

$$a_{\mathrm{P}}=a_{\mathrm{E}}+a_{\mathrm{W}}+a_{\mathrm{N}}+a_{\mathrm{S}},$$

(20)

$$b=\overline{S}{r}_{\mathrm{m}}\Delta r\Delta x.$$

(21)

The coefficients $a_j$ represent the fluxes crossing the faces of the control volumes. The fluxes are calculated using the Power Law scheme [24]. The subscripts N, E, S, W correspond to the compass directions, and represent the faces above, right, below and left to the currently analyzed node.

The SIMPLE algorithm is implemented into the numerical model, for calculation of the pressure corrections. The use of the SIMPLE algorithm is necessary to calculate fluid velocity in the specific regions of the reactor [24,45]. The systems of equations created during the discretization processes are further applied to the numerical procedure, and solved using the Gauss-Seidl iterative method [46].

3.1. Boundary Conditions

The presented numerical procedure is considered a quasi three-dimensional case. Thanks to the axial symmetry of the reactor, the computational domain is reduced to two dimensions [47]. The symmetry boundary conditions are applied at the axis of the reactor, simultaneously allowing for an acceleration of the computations. Afterward, thirty separate segments are designated. The segments divide the reactor into separate zones in the longitudinal direction. Each of the segments has th exact same dimensions and is predicted to be filled with a catalytic nickel/yttria-stabilized-zirconia (Ni/YSZ) composite or with a stainless steel metallic foam serving as a non-catalytic segment. The presence of the metallic foam segments serves to reheat the gases mixture before entering the subsequent catalytic segment. Therefore, the gases are carrying more energy, allowing for an intensification of the reforming process inside the proceeding catalytic region of the reactor. The prepared numerical procedure requires the definition of x adequate values of the thermal conductivity coefficients $\lambda$ for the two materials. Considering the process temperature taking values between 800 and 900 K, the $\lambda$ value for the Ni/YSZ is set at 22 W m${}^{-1}$ K${}^{-1}$ [48] and for the metallic foam at 30 W m${}^{-1}$ K${}^{-1}$ [49]. The reactors length L is equal to 0.3 m and the inlet velocity $u_{\mathrm{in}}$ is set at 0.15 m s${}^{-1}$. The reactor’s wall is assumed to be made of stainless steel and to have a thickness of 2 × 10${}^{-3}$ m. The non-slip and Neumann boundary condition, with the heat flux equal to 0, are set at the reactor’s wall. The feedstock entering the reactor is considered to acquire the reactor’s temperature immediately. The properties of specimens being the subject of the analysis are taken from the literature [50]. The numerical grid’s dimensions are set at 150 elements in the longitudinal direction, and 25 elements in the radial direction. A sensitivity analysis of the grid resolution and model’s falsification procedures are described in detail in our previous work [21]. All of the described boundary conditions are summarized in Figure 2.

4. Optimization Procedure

The optimization conducted in the presented analysis predicts the finding of the most optimal alignment of the catalytic and non-catalytic materials. The pursued optimal alignment has to lead to a unified distribution of temperature field without any significant losses in the produced hydrogen yield. The defined optimization problem is an arduous task due to the amount of possible segments configurations. Each segment can be distinguished as a catalytic or non-catalytic. Moreover, every single segment may vary in the values of porosity $\epsilon$ and average pore size $d_{\mathrm{p}}$. The large variety of possible solutions induces the presence of numerous local extrema. Due to the problem characteristics, the genetic algorithm (GA) is chosen as the optimization technique. The genetic algorithm is a mature optimization concept, with a stochastic nature [51]. The GA is reported to return satisfying results for optimizations regarding adversely conditioned search spaces [52,53,54]. The genetic algorithm is successfully applied in researches considering the optimization of the steam reforming reaction [39,55,56]. Due to the algorithm having its roots in the process of natural selection, vocabulary originating from genetics is introduced to the nomenclature. A single solution is referred to as a specimen. Each specimen consists of segments, going by the name of chromosomes. A single segment has an individual set of parameters, called genes. In the case of the presented research, genes represent if a single segment has catalytic properties and its values of $\epsilon$ and $d_{\mathrm{p}}$. A single iteration of the genetic algorithm includes the calculation of the reforming simulation for a set of specimens. The set is referred to as a population [57]. A numerical procedure is prepared for the needs of the presented research. The procedure starts with the preparation of the initial population. The chromosomes for the first set of solutions are initialized randomly. The algorithm is coupled with an in-house solver of the reforming process [39]. After randomization of the initial parameters, their values are loaded into instances of the reforming simulation. Afterward, the computation is commenced and proceeds until the convergence criteria for each specimen are met. The numerical procedure is assumed to converge after differences in values of the heat sources, acquired for two subsequent iterations, are lower than ${10}^{-5}$. The sufficiency of the given threshold is proven in our previous research [21,39]. After converging each specimen, evaluation of the results is performed, based on user-defined fitness functions [51]. The higher the overall fitness f of a specific reactor, the higher are its chances to pass its chromosomes to the proceeding population. Each subsequent population is prepared using the crossover procedure. First, a set of parent specimens is chosen, basing on the roulette rule [39,51]. Afterward, the algorithm randomly selects a pair of reactors from the parents’ set and translates their chromosomes’ values into a single binary string. Then, the strings are split at the randomly selected crossover point, and their corresponding parts are interchanged. The operation results in the creation of a pair of specimens, inheriting the chromosomes from two different parents. The crossover procedure is repeated until a population of the same size is created. After each crossover, a mutation procedure is performed [58]. The procedure predicts flipping of a single bit in the binary string, occurring with a probability specified by the user. The mutation is introduced to prevent the GA from finding a local extremum instead of the global one. Another benefit of the mutation is the possibility of introducing new genes, which acquiring would be impossible in the way of the sole crossover. After having the new population defined, the GA calls the reforming simulation over the created specimens and its procedure is repeated. The presented analysis predicted thirty subsequent iterations of the genetic algorithm. Each population is constrained to include thirty specimens. The number of reactors in a single generation was reduced to be only thirty, due to an enormous computation cost of the numerical simulation of the reforming process [39]. The goal of the presented optimization procedure is to define catalyst distribution, allowing for unification of the temperature field with simultaneous maintenance of the methane conversion, at levels similar to the conventional reforming reactors. To define the robustness of the acquired optimization results, two different fitness functions are defined. The first function analyzes the amount of methane converted during the process. The amount of methane processed during the reaction is calculated basing on the methane (CH${}_4$) fraction at the inlet and the outlet of the reformer, according to the following formula:

$$f_{{\mathrm{CH}}_4}=\frac{fra{c}_{{\mathrm{CH}}_4\mathrm{in}}-fra{c}_{{\mathrm{CH}}_4\mathrm{out}}}{fra{c}_{{\mathrm{CH}}_4\mathrm{in}}},$$

(22)

where $fra{c}_{{\mathrm{CH}}_4\mathrm{in}}$ and $fra{c}_{{\mathrm{CH}}_4\mathrm{out}}$ stand for the methane fraction at the inlet and the outlet of the reactor, respectively. The second function considers minimization of the temperature gradients occurring inside the reformer. The algorithm analyzes the local temperature differences for each of the created control volumes (Section 3). The procedure iterates over each of the nodes and computes values of the temperature T differences between the analyzed node P and its neighbors: N, E, S, W. When all the differences are calculated, the highest difference is set as $\Delta T_{{\mathrm{loc}}_{i,j}}$ of an analyzed control volume, where i and j represent the column and row number of the numerical grid (Equation (23)).

$$\Delta T_{{\mathrm{loc}}_{i,j}}=\max \left\{|T_{{\mathrm{P}}_{i,j}}-T_{{\mathrm{N}}_{i,j}}|,|T_{{\mathrm{P}}_{i,j}}-T_{{\mathrm{E}}_{i,j}}|,|T_{{\mathrm{P}}_{i,j}}-T_{{\mathrm{S}}_{i,j}}|,|T_{{\mathrm{P}}_{i,j}}-T_{{\mathrm{W}}_{i,j}}|\right\}.$$

(23)

After having calculated $\Delta T_{{\mathrm{loc}}_{i,j}}$ values for each of the control volumes, the highest found is set as the temperature difference representing the whole reactor $\Delta T$, according to the following formula:

$$\Delta T=\max \left\{T_{{\mathrm{loc}}_{1,1}},T_{{\mathrm{loc}}_{1,2}},\dots ,T_{{\mathrm{loc}}_{i,j}}\right\}.$$

(24)

The overall temperature fitness $f_T$ is calculated using the following equation:

$$f_T=1-\frac{\Delta T}{\Delta T_{\max }},$$

(25)

where $\Delta T_{\max }$ represents the maximal temperature difference, occurring in the reactor with continuous and homogeneous catalyst. The $\Delta T_{\max }$ for the needs of the presented research is set at 25 K. After having computed the $f_{{\mathrm{CH}}_4}$ and $f_T$ values, the overall fitness f is calculated using Equation (26)

$$f={\omega }_1·f_{{\mathrm{CH}}_4}+{\omega }_2·f_T,$$

(26)

where ${\omega }_1$ and ${\omega }_2$ represent the weights. The weights values are 0.6 for ${\omega }_1$ and 0.4 for ${\omega }_2$. The summarized procedure of the genetic algorithm is presented as a block schema in Figure 3.

5. Numerical Results

The presented analysis included computation of the genetic algorithm procedure, for three different compositions of the inlet gas. The compositions are described in detail in

Table 3. Variants of the inlet gas compositions.

No.	CH${}_4$	CO${}_2$	H${}_2$O	SC	CC
(1)	23%	30%	46%	2.0	1.3
(2)	20%	40%	40%	2.0	2.0
(3)	18%	50%	32%	2.0	2.9

. The different ratios of chemical specimens are defined in the numerical code by setting corresponding values of steam-to-carbon (SC) and carbon-to-carbon (CC) ratios [23]. The reference cases (Figure 4a–c) are reactors with homogenous and continuous catalytic material. The porosity of the catalytic material is set at 50% for each of the reference cases. The temperature field for the cases has a noticeable decrease in the temperature values at the inlet of the reactor (Figure 4a–c). The decrease occurs due to the reaction’s activation. The temperature distribution becomes more uniform, closing to the end of the reactor. The unification is induced by the reaction progress. As there is less methane for conversion left, the reaction decreases in its intensity, consuming lower amounts of energy locally. The main goal of the carried optimization, is to unify the temperature field in the whole reactor, simultaneously preventing significant reduction of the methane conversion rate. The comparison of the temperature fields for subsequent generations in presented in Figure 5, Figure 6 and Figure 7. To provide insight into the whole GA progress, the temperature fields of example reactors for different biogas compositions, from the initial generation are presented in Figure 5a1, Figure 6b1 and Figure 7c1. A single reactor is chosen for each composition. The optimization process is carried out for thirty subsequent generations. Temperature fields for the most optimal solutions acquired for the 30th, 30th and 30th generations are presented in Figure 5a2–a4 (composition (1)), Figure 6b2–b4 (composition (2)) and Figure 7c2–c4 (composition (3)). The information on fitness values is presented in

Table 4. Fitness values and conversion acquired for the most optimal specimens in the subsequent generations and for the reference cases.

Gen.	Composition (1)			Composition (2)			Composition (3)
	$f_T$	$f_{{\mathrm{CH}}_4}$	$f$	$f_T$	$f_{{\mathrm{CH}}_4}$	$f$	$f_T$	$f_{{\mathrm{CH}}_4}$	$f$
REF	0.06	0.84	0.53	0.07	0.86	0.54	0.10	0.88	0.56
INIT	0.01	0.66	0.40	0.01	0.72	0.43	0.02	0.62	0.38
10th	0.17	0.57	0.41	0.04	0.61	0.38	0.29	0.52	0.42
20th	0.33	0.55	0.46	0.43	0.61	0.54	0.44	0.49	0.47
30th	0.61	0.43	0.51	0.67	0.50	0.57	0.64	0.52	0.58

. According to the data presented in Figure 5, Figure 6 and Figure 7 and in Table 4, the algorithm manages to enhance the thermodynamic conditions inside the reformer. A visible improvement is reported for each of the inlet gas compositions. Detailed analysis of Table 4 gives information that methane conversion is decreased for each biogas compositions, when comparing with the reference cases. The methane conversion’s decrease is directly caused by the reduction of the amount of the catalyst used for the most optimal solutions found (

Table 5. The amount of catalyst used in comparison with the reference cases.

Gen.	Composition (1)		Composition (2)		Composition (3)
	$\iota$	$\zeta$	$\iota$	$\zeta$	$\iota$	$\zeta$
REF	100%	0.40	100%	0.38	100%	0.29
30th	17%	0.94	10%	1.8	6%	2.17

$\zeta$—green hydrogen productivity (-); $\iota$—amount of the catalyst used.

). The exact values of chemical specimens conversion are presented in

Table 6. Conversion of methane and carbon dioxide.

Gen.	Composition (1)		Composition (2)		Composition (3)
	${\mathrm{CH}}_4$	${\mathrm{CO}}_2$	${\mathrm{CH}}_4$	${\mathrm{CO}}_2$	${\mathrm{CH}}_4$	${\mathrm{CO}}_2$
REF	84%	19%	86%	21%	88%	23%
30th	43%	2%	50%	4%	52%	8%

$\zeta$—green hydrogen productivity (-); $\iota$—amount of the catalyst used.

. Considering the scale of the catalyst use limitation, forces to analyze the results from a different perspective. The use of the catalyst is reduced by circa 90% for each case, while the methane conversion fell only by circa 40%. A green hydrogen productivity $\zeta$ is introduced for a better insight on the increase of the effectiveness of the biogas reforming. The $\zeta$ parameter is an exact ratio of the hydrogen output and amount of the used catalyst $\iota$ for a specific reactor (Equation (27)).

$$\zeta =\frac{{\mathrm{H}}_{2_{\mathrm{output}}}}{\iota }$$

(27)

The observed phenomenon clearly indicates improvement of the overall process effectiveness (green hydrogen productivity). Simultaneously confirming observations reported by Ricca et al. and Settar et al. [59,60]. Therefore, the presented results prove the macro-patterning concept, to be a valid strategy for an optimization of the biogas steam reforming. The acquired results indicate the necessity of defining of the minimal catalyst amount in the algorithm to prevent the limitation of the catalyst use of the reported magnitude. The methane conversion could be maintained by changing of the reactor dimensions [21]. However, an additional optimization procedure for defining the optimal radius and length of the reactor would be necessary.

6. Conslusions

The presented paper focuses on the optimization of the temperature distribution inside a biogas reforming reactor. The research pursues optimal distribution of the catalytic material to ensure enhancement of the reaction effectiveness. The optimization procedure succeeds in modification the reactors’ design for the creation of a more unified thermal field. The reported results indicate a decrease in the hydrogen yield. However, the amount of the catalytic material used for the most optimal cases was reduced by: 83% for composition (1), 90% for composition (2) and 94% for composition (3). Considering the ratio of the hydrogen yield and the reduction of the used catalyst amount, the research confirms the improvement of the thermodynamic conditions inside the reactor, to be a valid strategy for increasing the biogas reforming effectiveness. The capabilities of the presented optimization strategy are most visible in the case of the composition (3). The optimization carried out for the composition (3) not only leads to an improvement in the thermal conditions but also manages to maintain conversion rate at a similar level as for the initial population. Following the presented results, the optimization procedure requires improvement before future calculations. For acquiring results of a higher quality, a sensitivity analysis of the evaluation procedure has to be conducted. The reaction effectiveness is ranked by evaluation of the biogas conversion rate. Hydrogen output at the end of the process should be considered instead of the methane fraction. The influence of magnitudes of the ${\omega }_1$ and ${\omega }_2$ values has to be investigated. Currently, the calculation of a thirty-generation algorithm run takes about 500 h. Therefore, the robustness of the numerical procedure has to be improved before the conduction of the sensitivity analysis. During the optimization, the overall amount of the catalyst used in the reactors is reduced significantly. A constrain defining a minimal amount of the catalyst used should be added to the algorithm’s procedure to allow the unification of the temperature distribution, preventing losses in the amount of hydrogen produced.