A CFD Design Approach for Industrial Size Tubular Reactors for SNG Production from Biogas (CO₂ Methanation)

Abstract

A tubular reactor based on the disk and doughnut concept was designed as an engineering solution for biogas upgrading via CO₂ methanation. CFD (Computational Fluid Dynamics) benchmarks agreed well with experimental and empirical (correlation) data, giving a maximum error of 8.5% and 20% for the chemical reaction and heat transfer models, respectively. Likewise, hot spot position was accurately predicted, with a 5% error. The methodology was used to investigate the effect of two commercially available coolants (thermal oil and molten salts) on overall reactor performance through a parametric study involving four coolant flow rates. Although molten salts did show higher heat transfer coefficients at lower coolant rates, 82% superior, it also increases, by five times, the pumping power. A critical coolant flow rate (3.5 m³/h) was found, which allows both a stable thermal operation and optimum pumping energy consumption. The adopted coolant flow range remains critical to guarantee thermal design validity in correlation-based studies. Due to the disk and doughnut configuration, coolant flow remains uniform, promoting turbulence (Re ≈ 14,000 at doughnut outlet) and maximizing heat transfer at hot spot. Likewise, baffle positioning was found critical to accommodate and reduce stagnant zones, improving the heat transfer. Finally, a reactor design is presented for SNG (Synthetic Natural Gas) production from a 150 Nm³ h⁻¹ biogas plant.

1. Introduction

The world is currently experiencing an energy revolution without precedents, motivated by social, environmental, and economic drivers. Still, an ambitious target of a 100% renewable energy supply will never be attainable due to the transient nature of wind and solar energy [1]. As an alternative, Power-to-Methane (PtM) stands out as a promising option to absorb and exploit surplus renewable energy in the form of “synthetic natural gas” (SNG). The possibility of using the existing natural gas infrastructure to store and transport the SNG to the end-user confers a critical advantage over other concepts. PtM systems consist of: (1) a H₂ source (water electrolyser), (2) CO₂ source, and (3) methanation reactor [2]. H₂ and CO₂, are converted in a methanation reactor into a gas mixture of CH₄ and H₂O through the Sabatier reaction (Equation (1)).

$$C{O}_2+4H_2⟷C{H}_4+2H_{2}O \mathsf{\Delta }{H}_r^0=-165.0\frac{kJ}{mol}$$

(1)

Potential sources for CO₂ are diverse: cement, iron and steel industries, flue gas from fossil-fuelled power plants, atmospheric air or biogas [3]. However, the abovementioned CO₂ sources, with the exception of biogas, require an additional carbon capture/upgrading step, which furthers reduces the energy efficiency and increases the costs significantly [4].

1.1. Current Trends in Reactor Design

Since CO₂ methanation is the core step in PtM, this topic has been the focus of extensive research with a particular interest in reactor design and optimization [5]. In methanation processes, multi-tubular (fixed bed) reactors appear as the standard concept for their simple design and straightforward manufacturing. Because of the high exothermicity of the Sabatier reaction, the shell and tube configuration appears as a practical and effective engineering solution for temperature control [6]. As depicted in Table 1, tubular reactors remain the focus of the most recent research works. The conventional design methodology for multi-tubular fixed bed reactors assumes that a single tube may be considered an appropriate representation of the whole system. Heat transfer is calculated assuming a constant heat transfer coefficient or a fixed temperature at wall boundaries [5,6,7,8]. Consequently, both coolant flow and temperature distributions along the shell-side are neglected. This approach leads to the assumption that all tubes are exposed to the same cooling conditions. From an operational point of view, this can only occur if the coolant flow is parallel and completely uniform to the reactor bundle. In addition, energy efficiency demands coolant flow minimization to reduce pumping energy consumption. Sun and Simakov [9] asserted that mineral oils are unsuitable for tubular reactor cooling systems due to their incapacity to handle high-temperature flows, hence molten salts are recommended. Indeed, most methanation studies consider this approach, incorporating molten salts as cooling systems [10]. Along with the problem mentioned above, most single tube models impose high external heat transfer coefficients [9,10,11,12]. Remarkably, no technical basis is given to that choice; neither is the nature of the cooling systems declared. Such high external heat transfer coefficients may only be found in boiling cooling systems, not conventional systems based on molten salts or thermal oil [13].

1.2. Research Scope and Contribution

To the best of our knowledge, no multi-tubular methanation reactor involving coolant side flow has been modelled in detail so far. Even in the field of reactor modelling, works of this kind are scarce. The work of Hukkanen et al. [14] is notable for being one of the first simulations which considers CFD modelling of the coolant flow. They correctly conclude that isothermal shell-side coolant behaviour was not observed, and distinctive coolant zones were identified through the reactor length. Nevertheless, individual tubes were not really included in the reactor. Instead, reaction heat was simulated through a simplified methodology [15]. Since the physical presence of tubes is neglected, it was not possible to model a realistic shell flow. Wu et al. [16], using a segmental baffled shell and tube reactor, proposed a correlation-based approach to evaluate the performance of the considered design. Correctly, this author highlights the importance of understanding the flow phenomenon inside the shell to assess the causes of weaknesses in the cooling system. However, there is no information regarding the performance of the reactor in terms of yield, conversion or selectivity of desirable products. Satisfactory heat transfer and flow features must be reflected in terms of reactor performance parameters. Jiang et al. [17] developed a CFD simulation of a multi-tubular reactor based on the disk and doughnut concept for the propylene to acrolein process. Remarkably, a whole bundle reactor comprising 790 tubes was modelled, and the actual effects of baffles and tubes on coolant flow were analysed. Although the superior performance of the disk and doughnut concept over the segmented baffle was probed, chemical reactions were not included. Instead, its effects were replaced through polynomial temperature distributions at reactive tube walls. Moon et al. [18] developed a CFD model of a multi-tubular reactor for oxidative dehydrogenation of butene to 1,3-butadiene, based on the standard segmented baffle design. Unlike Jiang [17], detailed Lan muir-Hinshelwood-Hougen-Watson (LHHW) kinetics was considered, but coolant flow, although included, was not studied. Thus, the necessity to account for the cooling jacket and the cooling medium, in order to get closer to real industrial-scale, remains a relevant topic in fixed bed reactor modelling, as highlighted by several authors [4,7].

One of the most critical problems in CO₂ methanation is the presence of a hot spot and the subsequent thermal runaway in reactor scale-up design and operation. A proper design of the cooling system is therefore of great significance. Furthermore, to understand the causes of weaknesses in the cooling system, the flow pattern inside the shell should also be considered. Therefore, a realistic approach is mandatory to study the key aspects that affect coolant performance and hot spot control. The traditional “single-tube” approach makes it possible to predict the general patterns of reactors’ behaviour and performance; however, a full CFD model can also visualise the coolant flow and temperature fields in detail under industrially relevant conditions. In this research, a novel modelling approach of multi-physics simulation of chemical reaction, fluid flow and heat transfer based on ANSYS Fluent is proposed for multi-tubular fixed bed methanation reactors. Relevant transport phenomena (turbulence, heat transfer and complex chemical kinetics) for both reactive and coolant flow are covered through a cost-effective computational method. From a methodological point of view, this study may be seen as a contribution to the work of Jiang et al. [17], adding a realistic interaction between coolant flow and an actual chemical reaction heat source. Consequently, coolant flow optimisation should be considered an essential step in reactor design, effectively improving the temperature uniformity and heat removal in reactive tubes. The focus of this study is the design of a tubular reactor, its cooling system and the subsequent performance evaluation. A design case study is developed, based on the “disk and doughnut” configuration, addressing the current needs of methanation technology: small units for decentralised plants based on biogas and an optimal coolant system for stable and energetic efficient operation.

Table 1. Summary of published fixed bed methanation models.

Modelling Approach	Reactor Type	Dim.	Code	Ref.
Pseudo-homogeneous models based on intrinsic kinetics	Tubular	2D	COMSOL	[19]
	Tubular	2D	Presto-Kinetics	[11]
	Tubular	2D	COMSOL	[20]
	Tubular	2D & 3D	Fluent	[8]
	Tubular	1D	Matlab & Athena	[21]
	Metallic-Honeycomb	2D	COMSOL	[22]
	Tubular	1D	FlexPDE	[23]
	Tubular	1D	N/A	[24]
	Tubular-Structured	2D	COMSOL	[25]
	Tubular	2D	OpenFoam	[26]
	Tubular	1D	Fortran 90	[27]
	Tubular	1D	AMPL & IPOPT	[28]
Pseudo-homogeneous models with effectiveness factor	Tubular	1D	Python	[29]
	Tubular	1D	Matlab	[9]
	Tubular	1D	CONOPT	[30]
	Tubular	2D	Matlab	[7]
	Structured-wall	1D & 3D	COMSOL 1D, Fluent 3D	[31]
	Tubular	1D	Fortran 90	[6]
	Tubular	1D	Matlab	[32]
	Tubular	1D	Matlab	[12]
	Tubular	1D	Matlab	[33]
	Tubular	1D	Matlab	[10]
	Tubular	1D	Matlab	[34]
	Tubular	3D	Fluent	[35]
	Tubular	1D	Matlab	[36]
	Tubular	0D, 2D, 3D	COMSOL	[37]
Heterogeneous models with intraparticle mass balance	Tubular, Structured	1D & 2D	Matlab (1D) COMSOL (2D)	[38]
	Catalytic wall reactor	1D & 2D	COMSOL	[39]
	Micro channel	3D	COMSOL	[40]
	Tubular-annular	1D & 2D	COMSOL	[41]
	Micro-structured	1D	gPROMS ModelBuilder	[42]
	Tubular	2D	Fortran 90	[43]
	Tubular	2D	Fortran 90	[44]
	Tubular, Fluidized bed	1D	Matlab	[45]
	Tubular	1D	Matlab	[46]
	Tubular, low dt/dp	3D (PRCFD-DEM)	COMSOL	[47]

Table 1. Summary of published fixed bed methanation models.

Modelling Approach	Reactor Type	Dim.	Code	Ref.
Pseudo-homogeneous models based on intrinsic kinetics	Tubular	2D	COMSOL	[19]
	Tubular	2D	Presto-Kinetics	[11]
	Tubular	2D	COMSOL	[20]
	Tubular	2D & 3D	Fluent	[8]
	Tubular	1D	Matlab & Athena	[21]
	Metallic-Honeycomb	2D	COMSOL	[22]
	Tubular	1D	FlexPDE	[23]
	Tubular	1D	N/A	[24]
	Tubular-Structured	2D	COMSOL	[25]
	Tubular	2D	OpenFoam	[26]
	Tubular	1D	Fortran 90	[27]
	Tubular	1D	AMPL & IPOPT	[28]
Pseudo-homogeneous models with effectiveness factor	Tubular	1D	Python	[29]
	Tubular	1D	Matlab	[9]
	Tubular	1D	CONOPT	[30]
	Tubular	2D	Matlab	[7]
	Structured-wall	1D & 3D	COMSOL 1D, Fluent 3D	[31]
	Tubular	1D	Fortran 90	[6]
	Tubular	1D	Matlab	[32]
	Tubular	1D	Matlab	[12]
	Tubular	1D	Matlab	[33]
	Tubular	1D	Matlab	[10]
	Tubular	1D	Matlab	[34]
	Tubular	3D	Fluent	[35]
	Tubular	1D	Matlab	[36]
	Tubular	0D, 2D, 3D	COMSOL	[37]
Heterogeneous models with intraparticle mass balance	Tubular, Structured	1D & 2D	Matlab (1D) COMSOL (2D)	[38]
	Catalytic wall reactor	1D & 2D	COMSOL	[39]
	Micro channel	3D	COMSOL	[40]
	Tubular-annular	1D & 2D	COMSOL	[41]
	Micro-structured	1D	gPROMS ModelBuilder	[42]
	Tubular	2D	Fortran 90	[43]
	Tubular	2D	Fortran 90	[44]
	Tubular, Fluidized bed	1D	Matlab	[45]
	Tubular	1D	Matlab	[46]
	Tubular, low dt/dp	3D (PRCFD-DEM)	COMSOL	[47]

2. Reactor Model

2.1. Reactor Description

The proposed design must suit the special technical constraints of a cooperative 150 Nm³ h⁻¹ biogas plant [48]. It was deducted that the reactor must be a compact modular unit of simple maintenance and operation. In addition, shell-side design demands an even heat transfer effect among tubes. Hence, a disk-doughnut configuration was adopted for such considerations. This concept makes extensive use of radial flow (crossflow) across tubes, maximizing heat transfer rates in the reactor [49]. The uniform flow rate distribution into the shell-side through the disk-doughnut configuration is the premise to ensure uniform cooling and better reactor performance. In addition, radial design minimizes both shell-side pressure drop and baffle-shell leakage [50]. Modular reactor dimensions are deducted from TEMA [51] standards and the Heat Exchanger Design Handbook [49]. A modular length of 1 m was selected to allow easy handling of the equipment and further plant scaling as needed. An inner tube diameter of 25 mm has been set to fulfil the recommended ratio of reactor internal diameter to catalyst particle >10 for a 2.5 mm spherical particle, thus reducing wall effects. Based on Jiang et al. [17] results, a co-current cooling system was selected as it shows better effectiveness in hot spot control than counter-current. According to these requirements and constraints, a disk-doughnut reactor concept is presented in Figure 1.

Since CO₂ methanation is an equilibrium reaction, high conversion rates (≈100%) may be extremely difficult to achieve. The final configuration (length and/or number of modules) will be determined in function of desired product composition. A restrictive target quality of CH₄ ≥ 95%, H₂ ≤ 5% and CO₂ ≤ 2.5% has been adopted, as suggested by Guilera et al. [52]. An operational temperature of 523.15 K was selected to reach a good compromise between kinetic restrictions and hot spots formation, with a maximum operational temperature of 650 K, which corresponds to catalyst admissible operating temperature [32]. Operational pressure is one of the most important parameters to define. It significantly influences investment and operational costs and reactor performance, favouring reaction kinetics and shifting the reaction equilibrium towards products. In addition, it has been reported that higher pressures prevent carbon formation. At 11 bar, carbon begins to appear just above 788 K for nickel-based catalysts [11]. Jürgensen et al. [53] reports that, for biogas-based carbon sources, methane may hinder carbon conversion greatly at pressures ≈ 1 bar. This effect can be easily compensated at pressures above 8 bar, while at 10 bar, the impact of methane content is nearly zero. Therefore, an operational pressure of 10 bar was chosen for its multiple benefits in terms of reactor performance and the availability of industrial equipment that may fulfil the technical requirements of hydrogen or biogas compression. CO formation was neglected as supposed by related works [9,24]. A non-stochiometric H₂/CO₂ ratio of 3.8 has been selected in order to easily meet product gas quality requirements following the experience of Jürgensen et al. [53]. Shell and tube characteristics and operating conditions are summarized in

Table 2. Design and operational parameters for simulated multi-tubular fixed-bed reactor. Biogas composition from [53].

Parameter	Unit	Value
Shell and tube dimensions
Shell inside diameter	mm	345
Tube outside diameter	mm	25
Tube length	mm	975
Number of tubes	-	20
Tube pitch	mm	40
Number of baffles	-	13
Baffle spacing	mm	70
Operational parameters
Reaction pressure	bar	10
Inlet temperature	K	523
Cooling temperature	K	523
GHSV	h−1	3200
Limit op. temperature	K	650
Gas flow per tube	Nm3/hr	1.5
Feed gas composition
CH4 inlet mole fraction	mol/mol	0.28
CO2 inlet mole fraction	mol/mol	0.146
H2 inlet mole fraction	mol/mol	0.55
H2O inlet mole fraction	mol/mol	0.015
O2 inlet mole fraction	mol/mol	0.004

.

2.2. Model Setup

In this work, CO₂ hydrogenation over a commercial nickel catalyst in a 3D multi-tubular shell and tube type reactor was investigated through CFD simulations. To accurately validate the adequacy of the CFD model (conservation equations and kinetic model [54]), experimental results from Gruber et al. [37] were used for comparison. In addition, due to the high relevance of the coolant flow in methanation reactors, simulation results were used to calculate shell-side heat transfer coefficients at tube walls. CFD results were further compared against the Gnielinski correlation [55] to check the validity of the shell-side heat transfer coefficients. Validation simulation was split into two benchmark problems and sub-domains, each of which concentrates on a particular phenomenon: (i) Chemical Kinetics (tube side) to verify validity of the model selected for the reaction kinetics, and (ii) Heat transfer-flow dynamics (shell side) to validate the calculated heat transfer coefficients. Once the heat transfer and kinetic model were satisfactorily validated, both simulations were merged into a whole model that accurately represented mass and energy transport. Thus, to save computation time while focusing on the relevant phenomena (heat transfer and chemical reactions), the simulation was divided into four steps which are described below.

• Step 1: Benchmark simulation for single tube: A single tube model was first developed to check and validate the LHHW kinetic model as adopted in Fluent via a User Defined Function (UDF). As only chemical reactions were of interest at this stage, a 2D simplified model was adopted. This Single Tube model was validated against experimental data reported by Gruber et al. [56]. Since heat transfer is the focus of this research, the experimental temperature profile was chosen as a reference variable for kinetic validation purposes at single tube level. Model geometry and simulation conditions were implemented to reflect both real experimental reactor geometry and conditions.
• Step 2: Benchmark simulation for shell-side: A non-reactive tube bundle was fully modelled to represent best the heat transfer interface between reactive tubes and coolant. This means that all structures that alter the flow were present (baffles, tube bundle, shell walls, coolant inlet and outlet). Since the focus was on heat transfer between individual tube walls and coolant, no chemical reactions were considered. Instead, as Jiang et al. [17] proposed, constant wall temperature distribution was incorporated as a boundary condition, emulating the presence of chemical reactions. Surface heat transfer coefficients at tube walls were then validated against Gnielinski [55] correlation for tube bundles.
• Step 3: Mesh independency: After validation of both phenomena (chemical kinetics and heat transfer), these were merged into a full model, and mesh independency was assessed for four mesh sizes. Boundary conditions considered for the mesh independency test reflect nominal operation of the proposed design, as declared in Table 2.
• Step 4: Study Case: A modular multi-tubular reactor was designed to fulfil the operational requirements of a small biogas plant. As already stated, it is in the best interest of an efficient reactor to reach the desired conversion under safe and stable thermal conditions while pumping work is minimized. The influence of coolant flow and type regarding pumping energy consumption and hot spot temperature control were analysed. Furthermore, the benefits of the disk-doughnut configuration were discussed from a thermal and hydrodynamic perspective. Finally, the number of required modules and final reactor configuration were determined to fulfil the needed SNG quality and biogas production.

2.3. Governing Equations

Governing equations (

Table 3. Governing and constitutive equations used in this CFD model.

Reactive flow (tubes):

Gas phase continuity: $\nabla ·\left(\phi {\rho }_g\overrightarrow{v}\right)=0 \left(2\right)$

Gas phase momentum: $\nabla ·\left(\phi {\rho }_g\overrightarrow{v} \overrightarrow{v}\right)=-\phi \nabla p+\nabla ·\left(\phi \stackrel{=}{\tau }\right)-\overrightarrow{S} \left(3\right)$

Gas-solid momentum exchange: $\overrightarrow{S}=\left(\frac{{\mu }_{\phi }}{\kappa }\overrightarrow{v}+C_2\frac{1}{2}\rho \left|\overrightarrow{v}\right|\overrightarrow{v}\right) \left(4\right)$

Gas phase Energy: $\nabla ·\left(\overrightarrow{v} \left({\rho }_g{E}_g+p\right)\right)=\nabla ·(\phi {\lambda }_{\phi }\nabla T_g-\sum _i{h}_{i }{j}_i)+S_{h,rxn} \left(5\right)$

${\lambda }_{\phi }=\phi {\lambda }_g+\left(1-\phi \right){\lambda }_s \left(6\right)$

Species: $\nabla ·\left({\rho }_g\overrightarrow{v} Y_i \right)=-\nabla ·\overrightarrow{J_l}+A_f·R_i \left(7\right)$,

Diffusive mass flux: $\overrightarrow{J_l}=\sum _{j=1}^{N-1}{\rho }_g{D}_{ij}\nabla Y_i \left(8\right)$,

Gas phase equation of state: ${\rho }_g=R{T}_g \left(9\right)$

Reynolds number based on the particle diameter:

$R{e}_{dp}=\frac{{\rho }_{g}U{d}_p}{{\mu }_g\left(1-\phi \right)} \left(10\right)$

Coolant flow (shell-side):

Fluid phase continuity: $\nabla ·\left({\rho }_f\overrightarrow{v}\right)=0 \left(11\right)$

Fluid phase momentum: $\nabla ·\left({\rho }_f\overrightarrow{v} \overrightarrow{v}\right)=-\nabla p+\nabla ·\left(\stackrel{=}{\tau }\right) \left(12\right)$

Fluid phase energy: $\nabla .\left(\overrightarrow{v} \left({\rho }_f{E}_f+p\right)\right)=\nabla ·({\lambda }_{eff}\nabla T) \left(13\right)$

${\lambda }_{eff}=\left(\lambda +{\lambda }_t\right) \left(14\right)$

Baffles & tube walls:

Solid phase energy: $\nabla ·({\lambda }_{st}\nabla T)=0 \left(15\right)$

) reflect the physical models that represent transport phenomena inside reactive tubes and shell-side flow. A porous media model (PMM) was adopted to represent the reactive tubes as a good compromise between accuracy and computational resources [57]. Interaction between phases is included through source terms that account for momentum, energy and mass transfer. Momentum balance between phases is modelled using a momentum exchange term, Equation (4), derived from the Ergun equation [58]. A mean porosity value was adopted in momentum (Equation (3)) and energy (Equation (5)) governing equations to account for the presence of voids in the continuum. Since all species concentrations are of the same order of magnitude, a multi-component diffusion model [58] was adopted to account for the diffusive mass flux, (Equation (8)). Radial dispersion effects were neglected as plug-flow conditions were assumed, that is, the ratio of catalytic bed length to catalyst particle (L_bed/d_p) > 150 and internal reactor to particle diameter, d_t/d_p > 10 [59]. An equilibrium model is adopted to represent the energy equation (Equation (5)) based on an effective thermal conductivity computed as the volume average of fluid and solid conductivity (Equation (6)). Since (d_t/d_p > 10), wall effects were discarded [59]. Radiation heat transfer remains negligible since operational temperatures do not surpass 400 °C [60]. Regarding mass transport limitation, a constant activity factor $A_f$ was incorporated to the species source term in Equation (7) following the experience of Matthischke et al. [12]. Chemical reactions were incorporated into Fluent as User Defined Functions (UDF) in C language [61]. Turbulence in the shell-side was modelled using the realizable k-ε turbulence model known to accurately represent flows involving curved geometries and swirling flows [62]. Gas flow in the porous medium exhibited a low Reynolds number (143, according to Equation (10)), which corresponds to a transition flow regime [63]. According to Jiang et al. [17] and Zhang et al. [26] turbulence effects in the transition regime are not very strong. Therefore, a laminar model was adopted in reactive tubes. Rigorous reaction kinetics based on the experimental work of Koschany et al. [54] with parameters of Matthischke et al. [12] were considered. Among current CO₂ hydrogenation kinetics, the model of Koschany et al. [54] remains remarkably useful in engineering applications since it was derived for a state-of-the-art catalyst and considering industrial operational conditions. Moreover, derived kinetics expressions cover a wide range of relevant operational conditions from the differential regime to the thermodynamic equilibrium [54]. For a detailed description of the kinetics adopted the reader is directed to Appendix A and Appendix B.

2.4. Numerical Methods

Pressure and velocity coupling was achieved through the SIMPLE algorithm. Second-order upwind discretization schemes were selected for momentum, energy and species equations, and a standard scheme for pressure. Variable gradients were discretised through the least-squares cell-based method. Under-relaxation factors for momentum and mass balances were set at 0.6. For the energy equation, the selected value was 0.4. Convergence was checked as the consecutive decrease in residuals by at least three orders of magnitude, 10⁻⁶ and 10⁻⁵ for the continuity, energy and mass balance, respectively. For simulations involving coolant flow, wall y+ was monitored and evaluated to guarantee standard wall functions requirements (y+ > 30) and heat transfer coefficient on tube walls. Additionally, CO₂ molar concentration at the tube’s outlet was checked for full model simulations. All steady-state governing equations, Table 3, were discretised and solved by the finite volume method using the commercially available CFD code ANSYS Fluent 19.

3. Results and Discussion

3.1. Benchmark Simulation for Single Tube

In order to verify the applicability of the adopted kinetics model [54], a 2D axisymmetric single tube model was implemented considering the geometry and operational conditions from Gruber et al. [56] (20 kW, 10 bar case): reactor inner diameter and length 30 mm and 1000 mm, respectively, operational pressure of 10 bar, reactants and coolant temperature of 523 K, gas hourly space velocity (GHSV_n) of 2250 h⁻¹ and H₂/CO₂ ratio of 4. A convection boundary condition was adopted at tube walls. Heat transfer coefficient was taken from [32], based on nucleate boiling water. A comparative plot between the experimental run and the CFD simulation for the individual tube is shown in Figure 2. There is a good match between results using the adopted kinetic model and the experimental temperature profile. The maximum discrepancy between simulation and experimental data is 20% at the tube inlet and 5% at the hot spot. This can partly be attributed to the lack of information regarding thermophysical catalyst properties (thermal conductivity and heat capacity) and to the adopted dilution profile. Simulation results also showed a similar trend with the experimental data (Figure 2). Gas temperature increases quickly until it reaches the hot spot at 80 mm from the inlet. Then, temperature steadily decreases, resulting in an outlet temperature of ≈523.15 K, which is coincident with the coolant temperature. The temperature at the hot spot was overestimated by 5% (860 K), but its position matches the experiment. It is concluded that an activity factor of 0.1 fits well the intrinsic kinetics to experimental data, and therefore it was adopted in subsequent simulations.

3.2. Benchmark Simulation for Shell-Side

Research and correlations that may suit disk-and-doughnut multi-tubular reactors are scarce and may not be adequate for systems with a small number of tubes. A proper alternative should represent the heat transfer from coolant bulk flow to reactor walls as accurately as possible. Due to the intrinsic symmetry of the flow in a disk-doughnut configuration, it is possible to assume that each tube section between two consecutive baffles behaves like a single row of tubes in crossflow. The Gnielinski [55] correlation for a single row of tubes stands out as the best model resembling this flow conditions, and therefore, was adopted for heat transfer CFD results verification. The Nusselt number for a row of tubes was calculated according to Equation (22), considering a laminar and a turbulent contribution. The associated Reynolds number was calculated from Equation (16), where $l$ corresponds to the “streaming length” (the length of the flow path traversed over a single tube), the void fraction,$\psi$, which in turn depends on the transverse pitch ratio in the row, was obtained from Equation (18). Calculation of $U_w$, is based on the mass conservation principle and the effective crossflow area, as specified by Slipcevic [49]. The heat transfer coefficient was evaluated through Equation (23), while numerical heat transfer performance was evaluated through Equation (24) by ANSYS CFD Post after the temperature field had been obtained.

$$R{e}_{\psi ,l}=\frac{U_{w}l}{\psi \nu }$$

(16)

$$10<R{e}_{\psi ,l}<{10}^6 Pr=\frac{\nu }{\alpha }$$

(17)

$$0.6<Pr<{10}^3\psi =1-\frac{\pi }{4(t_p/d_t)}$$

(18)

$$l=\left(\frac{\pi }{2}\right)d_t$$

(19)

$$N{u}_{l,lam}=0.664+\sqrt{R{e}_{\psi ,l}} \sqrt[3]{Pr }$$

(20)

$$N{u}_{l,tur}=\frac{0.037R{e}_{\psi ,l}{}^{0.8}Pr}{1+2.443R{e}_{\psi ,l}{}^{-0.1}\left(P{r}^{\frac{2}{3}}-1\right)}$$

(21)

$$N{u}_{O,row}=0.3+\sqrt{N{u}_{l,lam}{}^2+N{u}_{l,turb}{}^2}$$

(22)

$$h_{row}=\frac{N{u}_{O,row}{\lambda }_f}{l}$$

(23)

$$h_{num}=\frac{Q_w}{A\left(T_w-T_{bulk}\right)}$$

(24)

$$P_h=\frac{q\rho g{h}_w}{{3.6}^6}$$

(25)

$T_{bulk}$ is the flow weighted average temperature of the coolant fluid and $T_w$ the surface temperature near the external wall side $Q_w$ is the heat flux across tube walls and $A$ the tube wall surface area. Six test problems were performed for shell-side flow in the range of 1.7–4.6 kg/s. In Figure 3 the selected arrangement (mesh model and boundary conditions) is presented. Boundary conditions for shell-side test problems summarize as follows: Mass flow rate in kg/s at 473.15 K in coolant inlet (I), constant temperature at tube walls of 495 K (II) and adiabatic external walls (III). Thermally coupled interface boundary conditions were adopted when fluid (IV) and solid (V) cell zones (baffles) match. Operational pressure was set at 1 bar. The coolant chosen was thermal oil. Figure 4 gives the comparisons between correlation data and the shell-side simulation data. It reveals a good agreement between Gnielinski correlation and the CFD results. A maximum error of 10% between numerical and correlation values (at 4.5 kg/s) was obtained, which is acceptable for engineering applications. As a result, it can be concluded that the developed numerical approach was accurate enough to calculate tube to coolant heat transfer coefficients and thus effectively applicable to investigate the effect of the coolant flow on reactor performance in the considered range. For lower coolant flows, there was not a good agreement between Gnielinski correlation and CFD results. This becomes especially relevant for single tube simulations, in which the external heat transfer coefficient is calculated from correlations such as Fache [6]. In those situations, the range of validity of such correlations must be clearly stated.

3.3. Mesh Independency

Geometrical models and mesh systems were created with ANSYS Design Modeler and ANSYS Meshing. As shown in Figure 5, tetrahedral grids were adopted as the main mesh structure, and a non-conformal mesh approach has been considered to account for fluid (coolant, V), solid (baffles, VI) and porous media (reactor tubes VII) sub-domains. Different cell zones are coupled together through interface boundary conditions to allow heat transfer. Coolant cell zones (fluid domain) were modelled by means of tetrahedral grids, which were made finer in near wall regions through appropriate inflation layers. A structured (hexahedral) mesh was adopted for reactor tubes (porous media). These were made finer near tubes inlet to capture the chemical reaction phenomena and the heat coming from the catalyst bed to the fluid (coolant) domain. Refinement has been performed until the y+ value of the first layer grid meets standard wall functions requirements. Mesh independence tests are carried out considering four mesh densities (

Table 4. Mesh independence test results.

Mesh Count	Tube Average Outlet Temperature (K)	Tube Average CO2 Outlet (mol Fraction)
4.2 × 106	251.768	0.01430
4.7 × 106	251.764	0.01428
5.1 × 106	251.759	0.01425
6.8 × 106	251.748	0.01419

). Average CO₂ mole fraction at tube outlets and hot spot temperature inside tubes were calculated as representative variables. All mesh systems differ slightly in terms of maximum temperature and CO₂ mole fraction as well. As a result, the first mesh (4.2 × 10⁶) system was adopted to represent the physical model considering the balance between accuracy and workload. Boundary conditions for mesh independence summarize as follows: Thermal oil flow rate of 7 m³/h (1.7 kg/s) at 523.15 K in coolant inlet (I), thermally coupled interface boundary conditions at coolant-tube and baffle-coolant interfaces (II), adiabatic external walls (III) and pressure outlet (IV). Operational pressure has been set at 1 bar for coolant cell zones and 10 bar for tube cell zones.

3.4. Design Study

Parametric studies were conducted to investigate the effect of coolant flow rate and coolant type (molten salts and thermal oil). In addition, key reactor performance parameters were evaluated: (1) maximum tube temperature, (2) CO₂ conversion, (3) shell-side flow pressure drop and (4) average coolant temperature. CFD simulations were performed at coolant (molten salts and thermal oil) feed rates of: (1) 14, (2) 7, (3) 3.5 and (4) 1 m³/h. Scalable wall functions were implemented for near wall numerical calculation to avoid problems of successive mesh refinements due to different y+ values in each case. Furthermore, the impact of coolant type on reactor performance was considered in terms of consumed specific pumping power (kW of pumping power/kmol CH₄ produced) and hot spot temperature. Two commercially available coolants were selected: thermal oil [64] and molten salt [65]. Disk and doughnut potential in controlling the hot spot was also analysed by comparing two baffle configurations. Although numerical heat transfer results showed a good agreement with Gnielinski correlation in the range of 7–18 m³/h, lower flow rates (1–3.5 m³/h) were also considered for their relevance in hot spot generation and control. However, for the sake of design safety, results in the lower flow bracket should be taken only qualitatively.

3.4.1. Coolant Type Analysis

After the successful validation of the CFD model, the effect of coolant type and flow was investigated. First, by contrasting the coolant type and volumetric flow, the effects of these factors on reactor performance (heat transfer, CO₂ conversion and pumping power) were assessed. Performance parameters predicted by the CFD model are presented in

Table 5. Reactor Performance comparison between coolants (molten salt vs. thermal oil).

Coolant Flow and Type	Max Temperature (K)	Tube Average Temperature (K)	Coolant Average Temperature (K)	CO2 Conversion (%)	Shell-Side Pressure Drop (Pa)	Average Heat Transfer Coefficient (W/m2 K)	Specific Pumping Power (kW/kg mol CH4)
Molten salt Thermal oil (1 m3/h)	890.9 896.3	590.9 616.1	528.3 541.2	93.0 94.0	29 12	104 57	56 23
Molten salt Thermal oil (3.5 m3/h)	885.6 885.2	566.7 571.5	525.5 526.9	91.2 91.7	491 218	235 184	3310 660
Molten salt Thermal oil (7 m3/h)	880.8 880.7	561.5 563.4	524.3 525.1	90.2 90.6	1954 871	413 338	26,880 5360
Molten salt Thermal oil (14 m3/h)	880.1 879.5	558.3 558.9	523.7 524.1	89.4 89.4	7750 3480	760 645	215,500 43,250

. Regarding maximum (hot spot) temperature, there was a difference of 20 K between maximum and minimum flow rates, but no significant differences among coolant types for equal flows. Likewise, average tube temperature exhibits an expected difference between maximum and minimum flow rates (58 K). The lowest flow rate pair also shows a significant difference (max. 26 K) among coolant types. This behaviour is clearly explained from tube to coolant average heat transfer coefficient results, where molten salt exhibits a better performance in terms of heat transfer rate, especially for low coolant flowrates (5% superior to oil at 0.2 kg/s). This effect becomes less pronounced as the flow increases, since molten salt’s thermophysical properties represented in the Prandtl number, Equation (17), play an important role in increasing heat transfer at low Reynold numbers (Equations (20) and (21)). However, the most important difference between the selected coolants lies in the field of energy efficiency. The thermophysical properties which make molten salts more effective in low-flow regimes (kinematic viscosity and density) increase the pumping power required by five times (Equation (25)), with respect to thermal oil, for the same flow rate. Consequently, positive heat transfer properties of molten salts are greatly overcome by their unsatisfactory hydrodynamic performance. Therefore, thermal oil was chosen as the standard coolant for subsequent simulations.

3.4.2. Coolant Flow Rate Analysis

In this study, one of the most critical problems is the elucidation of the minimum thermal oil flow rate at which the temperature control problem is efficiently solved. CFD simulations were carried out at thermal oil feed rates of: (1) 0.2, (2) 1, (3) 1.7 and (4) 3.5 kg/s (1, 3.5, 7 and 14 m³/h respectively). As coolant mass flow varies, the minimum flow rate to control the temperature at hot spot position is thus determined. Figure 6a reports the simulated axial temperature profile along reactor length, while Figure 6b shows the radial temperature profile at hot spot. Based on simulation results, it is possible to position the hot spot at 75 mm from the tubes inlet, as seen in Figure 6a. As heat transfer coefficient increases, there is an obvious decrease on the reactor temperature profile. This improvement in reactor performance is less evident for faster coolant flow rates. The highest coolant flow rate (14 m³/h) shows only a marginal difference against the 7 m³/h temperature profile. The same trend applies for shell-side average temperature, where the fastest flow implies a nearly isothermal operation for the coolant, while at the lowest feed rate, an increase in nearly 50 K is appreciated. This suggests that subsequent increments in coolant flow rate beyond 3.5 m³/h does not promote hot spot control significantly and just slightly contribute to maintaining the average tube temperature at an enormous cost in additional pumping power. Only the slowest flow rate (1 m³/h) reaches a significantly higher outlet temperature (549 K), while all other cases tend to reach the operational coolant temperature ≈ 523.15 K. Figure 7, depicts the average heat transfer coefficient at tube surface and the hot spot temperature at Z = 400 mm from the inlet. While maximum temperature is well managed below 923 K in all cases, as coolant flow goes below 3.5 m³/h, a drastic increase in maximum temperature was observed, suggesting that an average heat transfer coefficient of less than ≈200 W/m²K promotes thermal run-away and therefore it can be considered as a limit to allow a safe and stable operation. This drastic decrease in temperature observed at the critical 3.5 m³/h flow (Reynolds 6700) may be attributed to the fact that this heat transfer rate closely overcomes the heat generation rate due to chemical reactions ($S_{h,rxn}$, Equation (7)). Although a flow of 3.5 m³/h (1 kg/s) appears as an optimal operational point, it was not considered for subsequent simulations (Flow field and Species distribution) since it lies out of Gnielinski correlation validity (Figure 4). Therefore, a coolant flow of 7 m³/h (1.7 kg/s) was adopted as the optimal design coolant flow.

3.4.3. Coolant Flow Field Analysis

Figure 8 describes the longitudinal coolant velocity field at reactor middle plane (yz) x = 0 for a feed rate of 7 m³/h (1.7 kg/s) and thermal oil as the fluid of choice. As coolant flows proceed across all baffles, it is possible to distinguish zones in crossflow (1) and vortex formation through doughnut openings (2). Crossflow zones, while maximizing heat transfer, also generate intermittent recirculation zones, which lead to the existence of stagnant areas of low heat transfer (3). Nevertheless, a regular flow pattern assures that all tubes are subjected to the same cooling conditions across the reactor length, guaranteeing performance uniformity among individual tubes. Indeed, coolant flow distribution suggests that an additional row of tubes may be added to make use of the parallel flow zones (4), contributing to process intensification in terms of CH₄ production. Figure 9 shows the characteristic transverse sections of coolant flow in the Disk and Doughnut reactor. Figure 9a shows the radial flow velocity vector at the hot spot position. Clearly, tubular reactors benefit from the Disk and Doughnut configuration since the baffle arrangement maximizes radial velocity at the first disk. This greatly improves the heat transfer capability if the first disk is correctly positioned to match the hot spot, as shown in Figure 10. Conversely, Figure 9b is positioned at the recirculation zone, showing poor flow features due to stagnation. Figure 9c shows a uniform crossflow zone among all tubes converging through the doughnut, while Figure 9d illustrates the expansion zone at the doughnut’s outlet. Both zones maintain a uniform crossflow pattern of an average of 0.2 m/s between tubes, which promotes turbulence (Re ≈ 14,000) and heat transfer as well.

Figure 10 illustrates the temperature distribution along the longitudinal section x = 0, for two alternative baffle positions at nominal operation conditions and thermal oil coolant flow of 7 m³/h. Figure 10b shows the reactor with baffle distancing as recommended by Thulukkanam [49], while Figure 10a shows an optimized baffle positioning to fit the hot spot location. The temperature distribution is visibly related to the baffle position, redirecting the main coolant flow radially as it expands from the inlet. Fluid velocity and heat removal are thus maximized, reducing the temperature gradient at the hot spot. Therefore, baffle positioning seems critical in the reactor design process to accommodate and reduce stagnant zones, improving the heat transfer effects.

3.4.4. Species Concentration Analysis

As it is characteristic of all exothermic processes, most of the reaction takes place in the first quarter of the tubular reactor. This can be appreciated from the drastic reduction in CO₂ concentration shown in Figure 11. The high reactants concentrations greatly promote reaction kinetics and reaction heat released. In this section, heat control is essential to guarantee a safe reactor operation and minimize carbon deposition and sintering.

At z = 0.575 m from the inlet, it is evident that as reactants concentration diminish, heat generation decreases accordingly due to slower reaction rates. In this reactor section, kinetics becomes less relevant as the reaction (Equation (1)) approaches equilibrium. Consequently, thermodynamic effects gain relevance as the reactor temperature reaches the coolant temperature. Therefore, it is critical to ensure that the chosen catalyst remains active at coolant temperature in this second section. While heat removal in the “kinetic” section is improved at lower coolant temperatures, in the “thermodynamic driven” section, it may hinder CO₂ conversion due to diminished catalyst activity.

Figure 12 and Figure 13 reports the resulting CH₄ and CO₂ mass fraction at the Y = 0 longitudinal cut at nominal conditions and a 7 m³/h thermal oil flow, showing a maximum methane mass fraction of 0.537 at reactor outlet, which means a methane mole fraction of 0.89 after water condensation as detailed in

Table 6. Outlet gas composition and SNG quality requirements (1 m reactor module).

Species (Inert Species Neglected)	SNG Target Composition % [52]	Outlet Comp. after H2O Cond. %	Mass Fraction at Outlet
CH4	≥95	89	0.54
CO2	≤2.5	2	0.036
H2	≤5	9	0.007
H2O	≈0	≈0	0.418

. Despite the reaction extent reached in a 1 m length reactor, it is not possible to achieve the quality required for the produced SNG to be injected into the gas grid. Therefore, an additional module is needed. Back to Table 5, it is important to notice that slower coolant flow rates promote CO₂ conversion up to 5% more than the nominal case (thermal oil 1 m³/h) due to a poorer heat removal capability, resulting in an increase of 50 K in the average reactor temperature. Although this may seem like an alternative to improve CO₂ conversion, thermal runaway risks make this option not recommendable. Therefore, a second 0.5 m reactor of the same configuration described in Table 2 is proposed to convert the remaining CO₂ and fulfil SNG quality requirements after an intermediate condensation step, as proposed by El-Sibai et al. [30] and Witte et al. [21]. For simplicity, all water is removed and product gases after the first reactor are assumed as the second reactor inlet composition. The rest of the boundary conditions remain the same. As most produced water is removed, reaction equilibrium is shifted to products, improving reaction performance in this second step as shown in Figure 14 and Figure 15. After this second methanation step, final SNG quality reaches acceptable levels. See

Table 7. Outlet gas composition and SNG quality requirements (0.5 m reactor module).

Species (Inert Species Neglected)	SNG Target Composition % [52]	Outlet Comp. after H2O Cond. %	Mass Fraction at Outlet
CH4	≥95	97.8	0.9475
CO2	≤2.5	0.42	0.0044
H2	≤5	1.68	0.0008
H2O	≈0	≈0	0.0472

below. The proposed design proved able to upgrade biogas to SNG and under thermally safe operational conditions.

4. Conclusions

In this study, the relevance of design features and coolant types on a fixed bed multi-tubular reactor were evaluated through a CFD-based methodology. A comprehensive disk and doughnut design was used for parametric studies involving thermal oil and molten salts for four different coolant flow rates. The adopted kinetics parameters can satisfactorily represent the hydrogenation of CO₂ under the considered operational conditions, with a maximum error of 20% at tube inlet. Moreover, hot spot position is accurately predicted with a 5% error between numerical and experimental data. CFD-calculated heat transfer coefficients coincide with Gnielinski correlation for a coolant mass flow range of 1.7–4.6 kg/s (7–18 m³/h for thermal oil) with a maximum error of 10%. Tube to coolant heat transfer coefficients calculated from empirical correlations should consider the range of validity of such correlations to represent the actual tube to coolant heat transfer. At the lowest coolant flow (1 m³/h), heat transfer coefficients obtained for molten salt cooling were 82% higher than thermal oil. Average temperature inside tubes for the molten salt-cooled reactor, were 26 K lower than the equivalent thermal oil-cooled reactor. Hot spot temperature difference between coolants was less pronounced (6 K). Regarding coolant maximum temperature, no significant increments were observed in all simulations (561 K max. for 1 m³/h thermal oil flow). Therefore, coolant decomposition is not expected. As the coolant flow increases, the heat transfer performance difference between coolants tends to diminish. A critical flow of 3.5 m³/h has been identified, after which no substantial improvements in heat transfer performance are observed. An increase in coolant flow from 3.5 m³/h to 7 m³/h implies an increment in eight times the pumping power in return for marginal improvements. An additional increment in cooling from 7 m³/h to 14 m³/h further increases the consumed pumping power, with almost no performance improvement. Among coolants, it is found that molten salt cooling increases the pumping power required per k_mol of methane produced by five times (in comparison with thermal oil) due to a greater viscosity and density. The disk and doughnut configuration proved highly efficient in terms of maintaining a uniform flow-field among tubes. All tubes are subjected to the same cooling conditions, no matter the position. Stagnation zones found behind baffles did not severely impact the cooling performance of the tubes; neither showed relevant increments in coolant temperature. Regarding hot spot temperature control, the first baffle position was critical to allowing the fastest coolant stream to match the hot spot’s intensifying heat transfer in the most demanded zone. CFD results demonstrated that faster coolant flow rates generated a sharper declining temperature profile in the reactor, diminishing CO₂ conversion (89.4% at max. cooling), while the slower coolant flow rate (1 m³/h) promotes the conversion up to 94%, due to a lesser heat removal capability resulting in an increase of 50 K in the average reactor temperature. Nevertheless, operation at such low flows is not encouraged due to the observance of thermal runaway behaviour. Although product gas at the reactor outlet does not fulfil SNG requirements, an intermediate condensation step and a second reactor of 0.5 m are incorporated, improving CO₂ conversion and enhancing CH₄ content. Final number of modules (1 m reactor + intermediate condensation + 0.5 m reactor) required to upgrade 150 Nm³/h of biogas are 30 and 1200 tubes in total.