A three-dimensional conjugate heat transfer model for methanol synthesis in a modular millireactor

(cid:1) A new modular millireactor (MMR) is designed and modeled by a novel approach. (cid:1) A conjugate heat transfer model is proposed to study methanol synthesis in MMR. (cid:1) Detailed analysis of heat transfer and methanol conversion in the MMR is studied. (cid:1) Methanol conversion is shown to vary between 9–23 % in the present MMR. (cid:1) The differences between 3D conjugate heat transfer and 1D/0D models are discussed.


Introduction
Methanol (CH 3 OH) is a versatile compound used both as a fuel and as a primary feedstock for the chemical industry (Ott et al., 2012).In combustion research, methanol is considered an alternative fuel for internal combustion engines (Verhelst et al., 2019) with the potential to be used in dual-fuel engines for marine shipping (Svanberg et al., 2018;Gadalla et al., 2021;Karimkashi et al., 2022).On the industrial scale, methanol is produced from natural gas or coal-based synthesis gas (syngas) (Bozzano and Manenti, 2016).Moreover, methanol can be produced in a sustainable manner using biomass-based syngas (Ohlström et al., 2074) or renewable hydrogen and captured carbon dioxide (Pontzen et al., 2011).However, currently, less than 1% of the global methanol demand, which reached 98 million tonnes in 2019, was supplied from renewables and biomass (International Renewable Energy Agency, 2021).The renewable methanol production process is not only a way to mitigate carbon dioxide emissions (Tran et al., 2016), but it is also an opportunity to store renewable electricity (König et al., 2015), utilize carbon dioxide (Soltanieh et al., 2012), and promote circular economy (Olah et al., 2009).
The largest industrial-scale methanol plants can reach over 9 kilotonnes in daily production levels (Bertau, 2014).Meanwhile, CO 2 -based plants are limited in size, the largest operating plant in Iceland produces 4 kilotonnes annually from renewable hydrogen and captured CO 2 (Harp et al., 2016).In 2023, a larger plant will start to operate using recycled CO 2 and H 2 from steel manufacturing with annual 100 kilotonnes methanol production (Carbon Recycling International, 2021).While syngas-based methanol production on an industrial scale is a technically mature and economical process, its synthesis from CO 2 faces drawbacks in various aspects (Zhong et al., 2020).Technical drawback derives from the excess water in the reactor due to hydrogenation of CO 2 (Maksimov et al., 2021).Moreover, economic disadvantage is a result of the combination of high price and low efficiency of water electrolysis for the production of hydrogen (Nyári et al., 2020).
Gas-phase methanol synthesis is conventionally carried out over copper-zinc (Cu/ZnO) based catalyst (Ali et al., 2015).The methanol synthesis is an exothermic process consisting of hydrogenation of CO and CO 2 and the reverse water-gas shift reaction.During methanol production, generated heat should be removed from the reactor for optimal catalyst activity, selectivity, and lifetime (Phan et al., 2016;Aartun et al., 2005;Yoshida et al., 2005).Thus, reactor design plays a critical role in high methanol yield (Rahimpour, 2008).Single-pass methanol yield is low, therefore, the process requires recycling of the non-reacted gases (Pontzen et al., 2011;Ott et al., 2012).At an industrial scale, this is a standard procedure, which is also reasonably cost-efficient.On the other hand, on small scale, any technical development in CO 2based synthesis could result in economic benefit.To achieve a higher single-pass yield, and therefore, increase efficiency, not only the catalyst needs to be tailored for CO 2 -based feed gas, but the reactor needs to be modified as well (Yang et al., 2017;Zhong et al., 2020).
Conventional industrial-scale methanol reactors can be divided into adiabatic and isothermal reactors (Bozzano and Manenti, 2016).A dual-type methanol reactor, a water-cooled reactor followed by a gas-cooled one, developed by the Lurgi Company is commonly used for large-scale production (Bertau, 2014).Most of the previous numerical works are based on optimization and modification of the syngas-based industrial-scale reactors, see Eksiri et al. (2020) and references therein.Note that most of the previous works assumed simplified one-dimensional models.Several works performed more detailed computational fluid dynamics (CFD) simulations.Mirvakili et al. (2018) performed full three-dimensional modelling of the industrial methanol synthesis reactor.Redondo et al. (2019) investigated various tubular reactor designs using 2D and 3D models.Son et al. (2019) investigated the performance of the proposed compact reactor in full 3D.Recently, Eksiri et al. (2020) investigated an axial-radial flow pattern in a modern plate reactor.On the other hand, only a few works investigated the CO 2 -based synthesis case.Arab et al. (2014) investigated the impact of catalyst structure on methanol production by comparing a packed-bed reactor to a monolithic one.They showed that the monolithic reactor shows better performance at high gas hourly space velocity while the packed-bed reactor is more appropriate at lower values.Samimi et al. (2018) compared CO 2 conversion to methanol in three different configurations, i.e., water-cooled, gascooled, and double-cooled tubular reactors.They observed that the gas-cooled reactor outperforms the other two configurations.
Microstructured reactors were proposed to overcome the restrictions of the conventional methanol synthesis reactors.These microchannel devices have a compact size with channel dimensions ranging from microns to millimetres.Moreover, these small-size reactors show efficient mass and heat transfer leading to improved energy efficiency and higher yield.Several experimental studies investigated methanol synthesis from syngas in different microstructured devices (Tonkovich et al., 2008;Bakhtiary-Davijany et al., 2011).However, there are only a few numerical studies on modelling the millireactors for methanol synthesis.Bakhtiary-Davijany et al. (2020) investigated the performance of an Integrated Micro Packed Bed Reactor-Heat Exchanger using a simplified 3D-pseudo homogeneous approach.Kyrimis et al. (2021) applied the CFD approach to study the role of direct CO hydrogenation in the fixed-bed reactor.Recently, Izbassarov et al. (2021) studied the effects of catalyst bed configurations under various conditions in a 2D fixed-bed reactor.The authors observed that the inline case showed the poorest performance among studied configurations.
Most of the previous studies considered the flow only in the packed bed, neglecting the coolant flow.Indeed, simulation is simplified when heat transfer due to cooling is modelled by imposing constant heat flux or temperature on the wall between the catalytic bed and the cooling channel.However, such simplification leads to inaccurate temperature distribution in the reactor and thus poor modelling of the overall system.
In the present work, a plate-shell type heat exchanger reactor with a simplified geometry, a modular millireactor (MMR) (Fig. 1), is proposed and compared to a conventional packed bed reactor (PBR).The present MMR is designed as a small laboratory-scale equipment and the physical dimensions are provided in Fig. 1.Due to the modular structure, the production scale can be increased when needed by adding more similar modules to the reactor.The compact size, modularity, and enhanced heat and mass transfer of the MMR are expected to be beneficial for CO 2 hydrogenation.We propose a full 3D computational fluid dynamics (CFD) model, a conjugate heat transfer (CHT) framework, for simulating the fully coupled heat transfer between cooling channels, stainless steel plates, and catalytic reactor.Finally, the performance of the MMR is investigated by utilizing the developed tool under various operating conditions.

Governing equations and discretization
In the present study, the conjugate heat transfer (CHT) approach is used to solve incompressible/compressible flow equations, species equations, and heat transfer equations in both, fluid and solid regions.The fluid in the cooling and reactor channels is modelled as incompressible and compressible, respectively.Note that constant thermophysical properties are assumed in the cooling channels and solid plates.In the reactor channel, the packed bed region is taken into account using the porous-media model due to the large tube-to-particle diameter ratio (N = 148).There are several assumptions considered in the porous-media model.Firstly, the exact particle geometry is not considered explicitly, but rather resulting flow resistance accounting for pressure drop is applied.Secondly, heat transfer is assumed to be under local thermal equilibrium between the fluid and the porous media.Thirdly, chemical reactions are assumed to occur volumetrically.Lastly, the resistance due to heat and mass transfer is ignored.Instead, such effects are compensated by incorporating the effective thermal conductivity k eff and diffusion coefficients D eff .The properties of the coolant and reactor bed are provided in Table 1.
In the fluid region, the continuity, momentum, species, and energy equations are given as: where u; p; Y k ; h s ; q; l; c p ; ; _ x k ; _ x h are the velocity, the pressure, the mass fraction of specie k, the sensible enthalpy, the fluid density, the fluid viscosity, the fluid heat capacity, the bed porosity, the reaction rate of species k, and the heat release rate, respectively.Note that is set to 1 everywhere except in the packed bed region.Here, I is the identity tensor and n s is the number of species.S is a source term added to account for the pressure drop due to porous medium.
Using the Darcy-Forcheimer law this term can be defined as follows: where d p is the diameter of the catalyst particles (Ergun, 1952).The effective properties of the catalytic bed consisting of gas and solid particles can be written as follows: where T is the fluid temperature, h t is the total enthalpy of the fluid defined as h t ¼ h s þ 0:5juj 2 ; k is the fluid thermal conductivity, k cat is the catalyst thermal conductivity, q cat is the catalyst density, c p;cat is the catalyst heat capacity, and s is the tortuosity factor of the bed.
Following Lanfrey et al. (2010), the tortuosity can be defined as where / is the particle shape factor.For all species, the diffusion coefficient is D ¼ k=ðqc p Þ as Lewis number is considered to be equal to unity.Note that the species formation/consumption rate ( _ x k ) and the heat source ( _ x h ) are defined only in the porous media and can be written as: where MW k is the molecular weight and h 0 f ;k is the enthalpy of formation of species k; a ik is the stoichiometric coefficient of component k and r i is the reaction rate of reaction i.Above, n r represents the number of reactions.Finally, the density of the mixture is calculated using the ideal gas law.The energy equation in the solid domain can be represented as: where T s ; q s ; c p;s ; k s are temperature, density, heat capacity, and thermal conductivity of the solid, respectively.Finally, the multiregion methodology requires appropriate boundary conditions at the solid-fluid interface:

Kinetic model of methanol synthesis
Following our previous work (Izbassarov et al., 2021), the formation of methanol was considered over a Cu/Zn/Al/Zr catalyst.The kinetic model for methanol synthesis derived by Graaf et al. (1988) is utilized here.The rate laws in their model describe the relationship between the rate of the chemical reaction and the fugacity of its reactants by a more complex Langmuir-Hinsel wood-Hougen-Watson (LHHW) type reaction rate.Graaf et al. (1988) model assumes that methanol production occurs by hydrogenation of CO 2 and CO linked together with the reverse water-gas shift reaction.The mechanism is based on three main reactions: Table 1 Main physical properties used in the current study.

Properties Value
Cooling oil SilOil M80.100 In our work, we apply the updated Graaf model (Graaf et al., 1988) as presented and validated by Kiss et al. (2016).The rate equations are: where f corresponds to fugacity [Pa] and the reaction rate r is in [kmol/kg cat s].All the constants are given in Table 2. Partial pressure is utilized instead of fugacity, due to the ideal gas assumption.The same LHHW kinetic model is implemented in Aspen Plus Ò and OpenFOAM as well.The details of the implementation can be found in Refs.(Izbassarov et al., 2021;Nyári et al., 2020).Note that the formation of by-products, e.g.dimethyl ether are not considered due to their insignificant amounts (Ushikoshi et al., 2000).

Overview of the computational setup
The governing equations in the fluid and solid domains are discretized using a second-order accurate finite volume method (FVM) using the open-source C++ library OpenFOAM (Weller et al., 1998).The default CHT solver chtMultiRegionFoam is modified to take into account the porous-media model.The thermal conductivity and diffusion coefficients are modified in energy and species equations, respectively.Moreover, transient porous media effects are accounted for in the continuity, energy, and species equations.The pressure-velocity coupling is done via the PIM-PLE algorithm, which is the combination of the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) and the pressure implicit splitting of operators (PISO) algorithms.Similar to our previous works (Izbassarov et al., 2021;Tekgül et al., 2020;Kahila et al., 2018), the flux-limited Gamma scheme of Jasak et al. (1999) is used to discretize the convective terms, while the diffusion terms are approximated using central differences.Finally, a second-order implicit scheme is utilized for temporal discretization.
Various modules of the present numerical tool have been utilized in different applications by the authors.Peltonen et al. (2019) studied non-reacting heat transfer without CHT in plate and pin fin heat exchangers confined in a pipe flow.Later, the non-reacting CHT solver was used by the present group in Laitinen et al. (2020) for cooling of high power electronics.Recently, Izbassarov et al. (2021) developed a novel computational tool for CO 2 and H 2 -based methanol synthesis.The method was used for particle-resolved numerical simulation of methanol production in a 2D fixed bed reactor.Here, we present a numerical tool that takes into account methanol synthesis in porous media within the CHT framework.
The non-reacting CHT solver was validated in our previous work (Laitinen et al., 2020).More recently, the catalytic solver was developed and validated for particle-resolved simulations of methanol synthesis (Izbassarov et al., 2021).In the present work, the catalytic solver is modified and combined with a porousmedia model to be used within the CHT solver.Firstly, we perform a grid convergence study.Next, we validate the CHT solver against the standard benchmark case.Finally, we study methanol production in the MMR under various conditions.
Altogether 52 simulations are performed to explore the effects of pressure and temperature on methanol production.Firstly, 43 simulations for the 1D PBR and thermodynamic equilibrium conditions are carried out.Following our previous work (Izbassarov et al., 2021), PBR is modeled using Aspen Plus Ò .The PBR is modelled as a plug flow reactor.The flow is assumed to be onedimensional, steady, and frictionless inside the reactor.The PBR has the same dimensions as the porous media region.The radius of the PBR is 37 mm and a length of 5 mm.The main PBR simulations are related to assuming that the reactor is always isothermal i.e. the produced heat and the lost heat balance one another.The sensitivity of this assumption is assessed in the validation section where also the adiabatic wall condition is considered.Additionally, 0D reference cases are used throughout the paper assuming thermodynamic equilibrium under isothermal conditions.
Secondly, 9 CFD simulations on the 3D MMR are performed using OpenFOAM.The main CFD simulations are made by assuming a full CHT model taking into account the reactants in the porous media, the solid part conduction, and the cooling channels.However, as part of model validation, we also investigate the possibility of applying a simpler constant heat transfer rate boundary condition.The CHT assumption is obviously the most accurate one since the heat transfer is directly computed by CFD and no extra assumptions are needed on the heat transfer surfaces.
The CHT method is employed to investigate the effects of pressure and temperature in reacting section on methanol production in the full MMR (Fig. 2).Here, the flow dynamics in the lower and upper cooling channels and the reactor catalytic channel are taken into account.Moreover, the heat transfer between channels and solid plates is considered.For both cooling and reacting sections, Dirichlet and Neumann boundary conditions are applied for the inlet and outlet for all variables.At all the walls, no-slip conditions are applied for velocity, while zero-gradient conditions are assumed for the remaining variables, except for temperature at heat transfer surfaces.Note that at the solid/fluid interfaces, the heat flux (Eqs.12) and the temperature (Eqs.13) are equal at solid and fluid phases.

Grid convergence study of the 3D computational setup
Firstly, grid convergence study is performed for the present CHT solver.The full 3D CHT simulations are performed for three different meshes as shown in Table 3 for the baseline case.Following the previous works (Kiss et al., 2016;Maksimov et al., 2020), the baseline case for the present catalytic reactor is defined via the following conditions: SV = 1.93 m 3 /kg cat h, P = 5 MPa, R = H 2 :CO 2 = 3, Re ¼ 1000; T in ¼ 503 K.The conditions for the cooling section are fixed as Re ¼ 1800 and T in ¼ 423 K. Here, the subscript in stands for the inlet.Mass fraction for the methanol at the outlet of the MMR are presented in Table 3 for the three meshes.As can be seen in the table, the results are very close for all three meshes considered in this work.Moreover, the average CH 3 OH profile in the porous bed and temperature profile positioned 0.02 m away from the center of the inflow region in the positive x-direction are plotted in Fig. 3.A small deviation in the average methanol mass fraction is observed for M2 at a radial distance less than 10 mm.Otherwise, all three meshes show very close results.It is decided to choose the finest mesh M4 for all results presented in the current study.

Validation: 3D MMR versus 1D PBR
Next, the numerical validation of the developed 3D reactive CHT model is done against a conventional packed-bed reactor (PBR).Two extreme cases for methanol production are considered in the PBR simulations: adiabatic and isothermal.Simulations are conducted at a fixed inlet composition R = H 2 :CO 2 = 3, T in = 503 K, P ¼ 5 MPa, and space velocity SV = 1.93 m 3 /kg cat h.For a consistent comparison between MMR and PBR, the validation takes into account the CHT from the reactants to the coolant while a constant cooling heat transfer rate is imposed from the cold plates.The CHT solver is validated against adiabatic (q = 0) and isothermal (q = À37.5 W) PBR predictions.Here, q is obtained from the Aspen Plus Ò simulation for the corresponding isothermal case and it is then used as input value to the CHT to provide a meaningful comparison between the two approaches.Comparison of mass fraction Y at the outlet of PBR and MMR reactors is shown in Fig. 4. As can be seen, the excellent agreement between PBR and MMR for both isothermal and adiabatic cases justifies the present 3D CHT CFD model.
Next, a more detailed analysis of the MMR is performed.Average quantities are calculated on the cross-sectional plane of   various radial positions.Mean temperature and methanol mass fraction profiles are shown in Fig. 5 for the adiabatic (q = 0) and q = À37.5 W cases.The profile can be divided into two separate zones: zone inlet (R 69 mm) and zone channel (R > 9 mm).Inside zone inlet , both T and Y CH 3 OH profiles decrease, reaching a local minimum at R % 9 mm.For q = 0, both mass fraction profiles monotonically increase in the zone channel reaching a plateau.Note that this increase is similar to the one typically observed in adiabatic PBR.
On the other hand, for q = À37.5 W, T and Y CH 3 OH profiles monotonically decrease and increase, respectively.Note that this continuous reduction in T is due to external heat flux resulting in a radial enhancement in methanol production.Indeed, high conversion of CO 2 occurs at low reaction rates due to exothermic and reversible methanol synthesis reactions.On the other hand, the reaction rate increases at high temperatures but the conversion is limited by thermal equilibrium.Therefore, the continuously decreasing temperature profile along the reactor is aimed for high conversion, due to the high reaction rate at the inlet of the reactor and enhanced thermodynamic equilibrium conversion at the outlet.
For further insight, the volumetric distribution of mean mass fraction Y and temperature T inside the reactive section of the MMR is studied.Volumetrically averaged radial profiles of T and Y CH 3 OH are presented in Fig. 6 and Fig. 7, respectively, for both adiabatic (q = 0) and constant heat transfer rate (q = À37.5 W) cases.In the former case, more than 90 % of the volume consists of hightemperature gases close to T ¼ 530 K which has a relatively low mass fraction of methanol (Y CH 3 OH = 0.085).In contrast, for the lat-ter case (q = À37.5 W), the temperature distribution is relatively flat in the range 490 6 T 6 530 K. Also, the mass fraction Y CH 3 OH shows a negatively skewed distribution varying in the range 0 6 Y CH 3 OH 6 0:15.Furthermore, CO 2 conversion (X CO 2 ) and CH 3 OH yield (Y CH 3 OH ) are studied inside the reactor.The normalized radial profiles of X CO 2 and Y CH 3 OH are defined as where subscript in denotes inlet.Next, X CO 2 and Y CH 3 OH profiles are plotted against the T profile for q = À37.5 W (Fig. 8).For comparison reasons, results obtained from Aspen Plus Ò for isothermal (1D) and thermodynamic equilibrium (0D) conditions are also provided in the figure.Conversion and yield values in the PBR increase with temperature (T < 510 K) due to an increased reaction rate.At T > 510 K, PBR results are equilibrium-limited.On the other hand, the radial profiles obtained from the MMR (3D) follow the thermodynamic equilibrium curve in the range of 490 < T < 520 K with an approximately constant offset.It is worth noting that in the CFD simulations the temperature decreases in the radial direction because of cooling.Thus the lowest temperature values are obtained at the outlet.As can be seen, the MMR outperforms the PBR locally at the exit.Note that at T > 520 K, there is a nonmonotonic trend which is related to the highly varying heat transfer close to the flow impingement point at the inlet.).Adiabatic (left) and isothermal (right) PBR cases are compared against MMR configuration with cooling at adiabatic (q = 0) and constant heat transfer rate (q = À37.5 W), respectively.Here, a constant heat transfer rate q = À37.5 W is used on the coolant side to allow consistent comparison against PBR.(SV = 1.93 m 3 /kg cat h, P = 5 MPa, R = H 2 :CO 2 = 3, Re ¼ 1000; Tin ¼ 503 K).

3D MMR:
Effects of pressure and temperature

Global flow characteristics
In this subsection, the results are presented for the MMR model.Unlike the previous subsection, besides considering a reacting channel with solid plates, cooling channels are also taken into account.The overall temperature distribution inside the MMR for the base case is shown in Fig. 3. Non-uniform distribution is observed due to heat removal through the cooling channels.The highest and lowest temperatures are observed in the middle and outer reactor channels, respectively.Heat conduction from the reacting channel to the cooling compartments takes place through solid plates.Fig. 9 (left) presents the dynamics of cooling oil in the upper channel.The oil enters from the left pipe, impinges and swirls near the surrounding walls.Tortuous streamlines result in the formation of a vortex system.The vortex structures disappear towards the outlet mainly due to the expanding region and low Reynolds number (Re = 1800) of the flow system.Fig. 9 (right) shows the mean temperature contours at the bottom of the upper cooling channel.The cooling fluid temperature increases as it moves along the channel.The impingement of the cooling fluid results in effective cooling of the surface.The hot spot occurs in the wake of the cylinder.In the lower cooling channel, the oil advances from the right to the left.Note that similar dynamics are observed for the lower cooling channel.
Fig. 10 shows flow dynamics inside the reacting channel.The left side of Fig. 10 displays the top view of the reactor coloured with average non-dimensional velocity contours.Reacting gas flows from the top towards the fixed bed and leaves from the side channels.The flow expands and contracts in the inlet and outlets of the fixed bed, respectively.The magnitude of the velocity of the gas decreases as it moves radially inside the fixed bed, resulting in lowvelocity regions.The velocity increases as it contracts at the entrance to narrow outlet channels.A high-temperature zone in the middle of the reactor can be observed inside the bed due to exothermic methanol synthesis reactions.As can be seen, the local temperature is inversely proportional to the local flow velocity.This behaviour can be attributed to a larger flow residence time at lower speeds.Larger residence time results in enhanced cooling and thus lower local temperature.To investigate the problem in more detail, the contours of mean temperature in the XY-plane are presented in Fig. 10 (middle).The low-temperature regions are observed near the edge of the reactor.It could be attributed to the effective heat transfer in low-velocity regions.As a result, the reduced radial distribution of the average temperature is observed.
Next, the contours of the mean methanol mass fraction Y CH 3 OH are plotted in Fig. 10 (right).The sample XY-plane is chosen in the middle of the reacting section.It is visible that methanol is produced only in the fixed-bed region.Indeed, as the gas mixture passes through the porous media, the Y CH 3 OH increases radially, reaching the highest value in the low-velocity zones.As the average temperature decreases along the radial direction inside the fixed bed, the conversion increases, and the reaction rate Fig. 6.Volumetric distribution of temperature in the porous media of MMR with cooling at adiabatic (q = 0) (left) and constant heat transfer rate (q = À37.5 W) (right).(SV = 1.93 m 3 /kg cat h, P = 5 MPa, R = H 2 :CO 2 = 3, Re ¼ 1000; T in ¼ 503 K).Fig. 7. Volume fraction of methanol mass fraction in the porous media of MMR with cooling at adiabatic (q = 0) (left) and constant heat transfer rate (q = À37.5 W) (right).(SV = 1.93 m 3 /kg cat h, P = 5 MPa, R = H 2 :CO 2 = 3, Re ¼ 1000; Tin ¼ 503 K).
decreases.On the other hand, local residence time increases as the velocity decreases along the streamline, approaching thermodynamic equilibrium similar to Fig. 8.

Temperature effects in the reacting section of the MMR
Next, the effects of the inlet feed temperature T in on the performance of the MMR are investigated.We perform simulations by varying the inlet temperature in the range of 483 6 T in 6 533 K while the other parameters are the same in the base case.Fig. 11 shows the average methanol mass fraction and temperature profiles for feed temperatures T in = 483, 503, and 533 K.As mentioned above, the profile can be divided into zone inlet (R 69 mm) and zone channel (R > 9 mm).Overall, the average temperature is nearly constant within the zone inlet , while it monotonically decreases inside the zone channel .The average T increases with the feed temperature in the inlet region, while it has a steeper decrease at higher feed values reaching approximately the same temperature towards the end of the zone channel .The ability to control the variation in temperature difference along the reactor has several advantages.Indeed, reduction of local temperature helps to operate at safer operating conditions, increase the lifetime of the catalyst, and possible enhancement of selectivity (Aartun et al., 2005;Yoshida et al., 2005;Phan et al., 2016).Moreover, it is observed that the feed temperature weakly affects local methanol production.The average Y CH 3 OH shows non-monotonic spatial variation, i.e. a decrease (zone inlet ) followed by an increase (zone channel ).Note that Y CH 3 OH changes non-monotonically with the inlet temperature at the entrance of the zone inlet .On the other hand, the maximum   value of Y CH 3 OH , at the end of the zone channel , increases with the feed temperature.
Fig. 12 shows X CO 2 and Y CH 3 OH profiles against T at various feed temperatures.Current results are compared with isothermal PBR and thermodynamic equilibrium predictions obtained with Aspen Plus Ò .Note that each point in the PBR case represents values in the outlet of a separate simulation at the corresponding temperature, whereas in the MMR, each point corresponds to the average quantity at the corresponding radial position of the same simulation at a particular feed temperature.Thermodynamic predictions for X CO 2 and Y CH 3 OH decrease with T, while results for the isothermal PBR case show a non-monotonic trend as explained above in Section 3.2.Both X CO 2 and Y CH 3 OH profiles for MMR monotonically decrease with T in the channel region (zone channel ), i.e.T < 500, 515 and 529 K at the feed temperatures 483, 503 and 533 K, respectively.On the other hand, the remaining non-monotonic region represents the inlet region (zone inlet ) as explained above.Within zone channel , both X CO 2 and Y CH 3 OH values are above the isothermal PBR case but limited by the equilibrium predictions.It is worth noting that both X CO 2 and Y CH 3 OH profiles shift upwards when the feed temperature increases.
Fig. 13 shows the X CO 2 and Y CH 3 OH outlet values for the present MMR (3D), the isothermal PBR (1D), and the thermodynamic equilibrium (0D) cases.Fig. 13 indicates that the conversion results for the MMR are similar to the corresponding PBR cases, except for T in = 483 K, where the conversion of the MMR model is larger than the corresponding PBR.Next, the effects of feed temperature on the methanol yield in the MMR and the isothermal PBR cases are discussed.The MMR surprisingly outperforms the PBR case for the range of feed temperatures studied in the current work.Note that the methanol synthesis reactions are equilibrium-limited.Indeed, the methanol yield for both the MMR and the PBR is bounded by the equilibrium data at inlet temperatures of 483 and 503 K.At the highest feed temperature (533 K), the methanol yield for the PBR is also limited by the equilibrium predictions.Surprisingly the results obtained with the MMR predict higher values than the thermodynamic equilibrium data at T in = 533 K and P = 5 MPa.It could be attributed to the fact that the thermodynamic equilibrium data is obtained assuming constant temperature, while the actual temperature distribution inside the MMR case is non-constant.The yield values for the MMR case slightly change with the feed temperature, unlike the non-monotonic variation of the PBR.At T in = 533 K and P = 5 MPa, Y CH 3 OH obtained with the proposed MMR predicts 1.71 more yield than the corresponding PBR case.Current results suggest that working at a high feed temperature (T in = 533 K) would be beneficial for the MMR.High-temperature results in an improved reaction rate at the inlet, the temperature gradually decreases along the reactor leading to enhanced methanol yield (see Fig. 12).

Pressure effects in the reacting section of the MMR
Finally, the effects of pressure inside the reacting section are investigated.Simulations are performed at P = 2, 5, and 6 MPa, while other parameters are the same as in the base case.Overall temperature distribution (not shown here) is similar to the one observed above (see Fig. 10), i.e., the temperature decreases radially with the lowest values in the low-velocity zone.Moreover, high pressure leads to an elevated temperature due to exothermicity.Thus, the hot temperature zone expands radially with pressure.Detailed inspection of the zone reveals that the maximum temper-  ature increases with P (Fig. 14).Overall, the methanol mass fraction increases radially.It is worth reminding that high pressure favours the reaction and the equilibrium results in higher methanol production.
Radial profiles for T and Y CH 3 OH at various P are plotted in Fig. 14.The temperature stays nearly constant in the inlet region, while it monotonically decreases in the channel zone due to cooling.On the other hand, Y CH 3 OH decreases in the inlet region, while it increases in the channel region due to higher conversion at a lower temperature.Note that lower local T results in a lower reaction rate, on the other hand, due to low local velocity the reaction time is enhanced resulting in higher Y CH 3 OH .It can be seen that both T and Y CH 3 OH profiles shift upward with pressure.The T distribution is qualitatively similar for all P = 2, 5, and 6 MPa.While the distribution of Y CH 3 OH is more affected by the pressure.It is observed that at P = 2 MPa, the Y CH 3 OH is almost constant throughout the reactor due to low pressure.On the other hand, distinctive non-monotonous distribution of the Y CH 3 OH is observed at higher pressures.
Finally, Fig. 15 shows the comparison of CO 2 conversion and CH 3 OH yield in the outlet of the reactor for P = 2, 5, and 6 MPa.The MMR results (3D) are compared with the corresponding isothermal PBR (1D) and thermodynamic equilibrium (0D) predictions.Both CO 2 conversion and CH 3 OH yield increase with pressure, as explained above.The predictions for both MMR and PBR cases are generally below the corresponding equilibrium values.CO 2 conversion values obtained with the MMR model are below the PBR case at P ¼ 2 MPa, while they approach the PBR case at P = 5 and 6 MPa.On the other hand, methanol yield (MMR) values are generally larger than the PBR ones but less than the equilibrium predictions.It could be attributed to the varying temperature distribution resulting in local different equilibrium values.Indeed, as can be seen in Fig. 14, the T decreases radially enhancing thermodynamic equilibrium conversion.

Conclusions
In the present work, a new numerical setup has been developed to study methanol production in a proposed modular millireactor.The model combines a porous media model for the reacting section and a conjugate heat transfer model for the solid parts.Compressible and incompressible flow solvers are utilized in the reactant and coolant parts.The developed solver was validated by comparing the results against the conventional PBR (1D) obtained with the Aspen Plus Ò .Overall results at the outlet of the reactor were found in good agreement.Next, the performance of the MMR was investigated by a detailed analysis of heat and mass transfer in the reactor.Simulations were performed by taking into account the fully coupled heat transfer between cooling channels, stainless steel plates, and catalytic reactor.The effects of pressure and feed temperature on the performance of the MMR were studied.The MMR (3D) results were compared with PBR (1D) and thermodynamic equilibrium data (0D).The main findings can be summarised as: 1.It was observed that the average temperature along the reactor part of the MMR was consistently affected by the coolant temperature.In particular, the thermodynamic equilibrium (0D) conversion was shown to act as a theoretical upper bound for the reactions.On the other hand, the equilibrium conversion was spatially fixed for the PBR due to the applied isothermal condition.In general, methanol obtained at the outlet of the MMR was higher than in the corresponding isothermal PBR but was mainly bounded by the equilibrium data.2. The effects of the feed temperature on the methanol production in the outlet of the isothermal PBR and the MMR were noted to differ from one another.The methanol yield was nearly constant at feed temperatures 483, 503, and 533 K for the MMR, while it showed a non-monotonic trend for the PBR.For cases at 483 and 503 K, the yield for the MMR was between the PBR and the equilibrium data.Surprisingly at 533 K, the yield for the MMR was observed higher than the thermodynamic equilibrium prediction.The observed difference was attributed to the non-uniform temperature distribution inside the MMR case.Indeed, it was shown that the local methanol yield inside the MMR was limited by the equilibrium prediction.3. The effects of pressure on the methanol production at the outlet of the isothermal PBR and MMR showed a similar overall trend.
The methanol yield increased with pressure as high pressure is favoured by the reactions.On the other hand, increased pressure resulted in increased temperature due to the exothermic character of the reaction.Due to cooling, the average temperature decreased along the reactor.Note that, the high temperature increases the rate of reaction, while the lower temperature increases the equilibrium predictions.Thus, the yield obtained with the MMR was higher than the isothermal PBR values.
4. The developed methodology could be effectively utilised for a comprehensive analysis of the performance of a catalytic reactor and its development.The benefits of 3D CFD modelling of the reactor design were presented in the current work.Note that the PBR is a one-dimensional case, while the MMR is three-dimensional.Thus, unlike in the PBR case, in the MMR case, there is a strong interaction of flow with chemical reactions.Detailed analysis revealed that the local mass and heat transfer were significantly affected by the flow dynamics.It was shown that at 533 K and 5 MPa, the MMR produced methanol yield by a factor of 1.71 more than the case of the isothermal PBR.
Finally, current work clearly showed the advantages of 3D CHT modeling of geometrically complex methanol reactors.However, further developments in the modelling of the porous media and the reactions would be required to better understand the modeling sensitivities.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.A schematic of the modular millireactor is presented in the present work.

Fig. 8 .
Fig. 8. Average temperature vs CO 2 conversion (left) and CH 3 OH yield (right) profiles (red solid lines) for the MMR with cooling for the CHT model.Here, a constant heat transfer rate q = À37.5 W is used on the coolant side to allow consistent comparison against PBR.The dotted and dashed lines represent results obtained for the isothermal PBR (1D) and thermodynamic equilibrium conditions (0D), respectively.(SV = 1.93 m 3 /kg cat h, P = 5 MPa, R = H 2 :CO 2 = 3, Re ¼ 1000; T in ¼ 503 K).

Fig. 9 .
Fig. 9. Upper cooling channel of the MMR.(left) Velocity streamlines in the upper cooling channel.The flow goes from left to right.(right) Average temperature contours on the bottom of the channel.(Re ¼ 1800 and Tin ¼ 423 K).

Table 3
Mesh convergence.Comparison of mass fraction obtained with 3D (OpenFOAM) in the outlet of the reactor.Below, NC stands for the total number of cells.