Computational Fluid Dynamic (CFD) and Reaction Modelling Study of Bio-oil Catalytic Hydrodeoxygenation in Microreactors

A Computational Fluid Dynamic (CFD) model was derived and validated, in order to, investigate the hydrodeoxygenation reaction of 4-propylguaiacol, which is a lignin-derived compound present in bio-oil. A 2-D packed bed microreactor was simulated using pre-sulphided NiMo/Al 2 O 3 solid catalyst in isothermal operation. A pseudo-homogeneous model was first created to validate the experimental results from literature. Various operational parameters were investigated and validated with the experimental data, such as temperature, pressure and liquid flow rate; and it was found that the CFD findings were in very good agreement with the results from literature. The model was then upgraded to that of a detailed multiphase configuration; and phenomena such as internal and external mass transfer limitations were investigated, as well as, reactant concentrations on the rate of 4-propylguaiacol. Both models agreed with the experimental data, and therefore confirm their ability for applications related to the prediction of the behaviour of bio-oil compounds hydrodeoxygenation.


Introduction
Biomass presents numerous advantages as a renewable feedstock for bio-fuels. It contains low sulphur and nitrogen; and presents its self as a lucrative alternative due to its lack of net carbon dioxide (CO 2 ) emissions to the environment 1 . Resources from biomass consist of a vast range of materials, such as organic waste products, forest and agricultural residues and energy crops 2 . Biomass feedstock's which contain cellulose, hemicellulose and lignin possess a high-energy content, and are often converted to oil using fast pyrolysis 2, 3 .
Fast pyrolysis is the process of rapidly heating biomass under moderately high temperature conditions of around 500 o C and short reaction times of 2 seconds, in the absence of oxygen. In return, biomass degrades to produce mainly vapours, aerosols and some solid char. After processing, a dark brown liquid is obtained which has a heating value of approximately half of that of conventional fuel oil averaging at about 30 MJ kg -1 4 . Biomass derived bio-oil has several disadvantages such as a low heating value, high viscosity and a high oxygen content, which all restrict its application as a liquid fuel. Therefore, further upgrading of bio-oil by hydrodeoxygenation (HDO) is required 5 .
The HDO process converts the oxygen containing compounds such as acids, aldehydes, alcohols and phenol to oxygen-free hydrocarbon fuels 6 . Bio-oil obtained from the fast pyrolysis of lignin contains approximately 39% of guaiacol and its derivatives. Amongst these constituents, guaiacol is often regarded as a representative model for bio-oil derived from lignin because it has two types of C-O bonds (Csp 2 OH and Csp 2 OCH 3 ) within its molecular structure 7 . Based on this, majority of studies have utilised the compound guaiacol as a model compound and have investigated the HDO of guaiacol using catalysts such as NiMo/Al 2 O 3 and CoMo/Al 2 O 3 , precious metal catalysts, such as platinum, ruthenium and rhodium, and nickel (Ni) catalyst 8 .
Phenolic compounds such as guaiacol, anisol and phenol have been extensively modelled in past studies due to their significant present in bio-oil 6 . Despite this, not much attention has been dedicated to lignin-derived compounds. Therefore, there exists a limited understanding of the reaction pathways and kinetics of the HDO reaction. The study of 4-propylguaiacol HDO has been selected as a basis to produce a model representing lignin-derived compounds 3 . 4-Propylguaiacol represents some of the key lignin-derived components present in bio-oil such as benzene, phenol, guaiacol, anisole, propyl anisole, propylphenol, and propylbenzene. The existence of phenolic compounds in the bio-oil is the origin of polymerisation and coke formation during HDO at elevated temperatures greater than 300 o C 3 .
to a variety of hydrocarbons, such as cyclohexane. They found that cyclohexane was the predominant product and the yield increased with increasing number of acid sites. In order to obtain bio-oil with the maximum yield of cyclohexane and alkylated cyclohexanes, the SI/AL molar ratio should be adjusted to balance the Pt particle-induced hydrogenation with acid site-induced methyl group transfer.
Patil et al. 9 studied the HDO of the model compound guaiacol in bio-oil to produce fuel grade oil, using a bimetal catalyst No-Mo/ZrO 2 -Al 2 O 3 . The results showed that the conversion of guaiacol and the product yield of phenol and cyclohexane was found to increase with increasing Mo content (10%, 15%, and 20%) at continuous Ni (4%) loading and Ni content (2%, 4%, and 6%) at continuous Mo (20%) loading. However, there was a lower difference in guaiacol conversion and phenol and cyclohexane yield at variable Ni loading. Guaiacol conversion was 100% at 330°C and 30 bars for the improved catalyst.
Taghvaei and Rahimpour 10 investigated the catalytic HDO of guaiacol via a combination of dielectric barrier discharge (DBD) and catalyst. It was found that the highest conversion of guaiacol (92%) and deoxygenation degree of 65% are achieved in the presence of Pt-Cl/Al 2 O 3 and Pt-Re/Al 2 O 3 catalysts respectively. The predominant products obtained were BTX, phenol, methylphenols and dimethylphenols. It was concluded that the difficulties of using hydrogen for the HDO reaction can be overcome by using the catalytic DBD reactor.
Liu et al. 5 investigated the HDO of bio-oil model compounds over amorphous NiB/SiO 2 -Al 2 O 3 catalyst. The performance of the catalysts was evaluated in an oil-water biphasic system using anisole and guaiacol as the selected model compounds of bio-oil. The results showed that increasing the reaction temperature or reaction time would enhance the conversion of guaiacol and anisole. The HDO pathways of guaiacol and anisole were studied which provided a reference for the HDO mechanism of bio-oil.
Microreactors present numerous advantages to investigate the reaction of HDO. Their enhanced surface-area-to-volume-ratio leads to a much-improved mass and heat transfer, in addition to, shorter residence time. Therefore, reactions which contain unstable intermediates are better suited to these type of reactors because of their stability and high degree of control 11 . These collective advantages of microreactors mean that they are greener and environmentally sustainable 12 .
In order to understand the effects of HDO on the further processing of bio-oil several studies have identified mathematical modelling exercises to fully simulate its interaction with compounds under various conditions [13][14][15][16] . The majority of HDO studies have focused on sulphided cobalt and nickel-based molybdenum (CoMo and NiMo) based catalysts for the separation of sulphur, nitrogen and oxygen from petrochemical feedstocks 17 .
The catalytic HDO reaction of 4-propylguaiacol to 4propylphenol using pre-sulphided NiMo/Al 2 O 3 catalyst is investigated and presented in this study using a packed-bed plug flow microreactor. A pseudo homogeneous model was produced, and good model validation was obtained with the experimental data. In this study, a 2-D Computational Fluid Dynamic (CFD) model was created to improve the understanding of the mass transfer and catalytic reactions taking place within the microreactor and provide an insight into any potential improvements that could be made by investigating these parameters. A validation of the model with the experimental data from literature 3 is shown, and further investigations such as, internal and external mass transfer are conducted herein.

Reaction kinetics
The reaction considered for the CFD mathematical modelling is the HDO of 4-propylguaiacol to 4-propylphenol (eq. 1). The kinetic studies from the experiment suggest that the kinetics for the reaction are determined by the surface reaction step which represent the competitive reaction with non-dissociative adsorption of hydrogen 3 . The kinetics were determined using the Langmuir-Hinshelwood-Hougen-Watson (LHHW) method, thus giving the rate equation shown below: (1) 10  where k is the kinetic rate constant of each reaction pathway (L mol -1 ) and C is the molar concertation at each stage of reaction mechanism. Readers are referred to Joshi and Lawal, 2013 for the reaction scheme. From the rate expression, it can be deduced that the reaction is pseudo-first-order with respect to 4-propylguaiacol and hydrogen. Table 1 displays the kinetic constants utilised for the study, and Table 2 conveys the preexponential factors, activation energies and heats of adsorption upon which the CFD modelling is based.

Pseudo-homogeneous model
The microreactor models were simulated using CFD to demonstrate the particle-fluid transport phenomena. Experimental work is typically costly and time-consuming, while multiphase CFD studies can deliver comprehensive information on the spatiotemporal variations in species flows, concentrations and temperatures within the reactor and with minimal effort. Therefore, CFD is a favourable approach/methodology predicting the parameters thus, enabling a detailed study of the physico-chemical processes involved 16 . CFD is an integrated methodology in the package that was used.

Please do not adjust margins
Please do not adjust margins  A 2-D microreactor model was created based on the assumption that the concentration and temperature gradients occur only in the axial direction. The only transport mechanism operating in this direction is the overall flow itself, and this is of plug flow type. Further assumptions which the model was founded upon are (a) Application of steady-state and isothermal conditions; (b) there is a pressure drop of 0.07 MPa along the length of the microreactor; (c) Henrys law applies is valid for the gas-liquid interface; (d) the ideal gas law applies for the fluids in the gas phase; and that (d) there is a constant axial fluid velocity in the gas phase with uniform physical properties and transport coefficients. The microreactor has a height of 0.762 mm, and a length of 180 mm. The catalyst used is a pre-sulphided NiMo/Al2O3 in the form of solid spherical particles. Figure 1 shows a schematic diagram of the model used for the simulation.

Mass balances in microreactor
Modelling of the liquid chamber was established by connecting the gas phase reaction zone with the bulk liquid phase. The mass balance equation for the 4-propylguaiacol in the liquid phase is expressed as: where c 4PG is the concentration of 4-propylguaiacol in the bulk liquid phase in mol/m 3 , D i is the molecular diffusion coefficient (m/s) in the bulk liquid, u x is the liquid velocity in the axial direction.
The mass balance in the gas phase reaction zone is expressed as: (4) where χ is the association factor of 4-propylguaiacol, 1 for nonassociated solvents, V H 2 is the molar volume of hydrogen at normal boiling point in m 3 /kmol 21

Activation energies and heats of adsorption
The mass balance equations coupled with the appropriate boundary conditions were solved using COMSOL Multiphysics ® software version 5.3. The finalized geometry comprised of a mesh consisting of 4400 domain elements and 400 boundary elements, and 9,090 degrees of freedom was used, and the results were found to be mesh independent with a computational time of 5 seconds.

Detailed multiphase model
The previous pseudo-homogeneous model shown in section 2.2 was further enhanced to incorporate the catalyst particles on which the gas phase reaction occurs. The assumptions for the pseudo-homogeneous model are applicable to this detailed model, with the exception that the reaction zone in figure 1 is packed with solid spherical catalyst particles of the same size and shape. The mass balance equation for the species in the catalyst bed is expressed as: where, D i,A and D i,T are the axial and transverse dispersion coefficients respectively, J i is the molar flux of i into the catalyst particles, S b is the specific surface area of the particles exposed to the reacting fluids in the packed bed (assuming randomly packed spherical particles) and is given by 22 : where, is the fractional voidage of the packed bed and S is the specific surface area, in m, of the particles. For spherical particles this is given by: where, r p is the catalyst particle radius.
An assumption of the film condition is realised at the pellet-fluid interface. The mass flux across this pellet-fluid interface into the pellet is potentially rate determined by the resistance to mass transfer on the bulk fluid side. This resistance can be expressed in terms of the external mass transfer coefficient: where, c i,ps is the concentration of i at the catalyst particle surface and h i is the external mass transfer coefficient. Sc is the Schmidt number, μ is viscosity of 4-propylguaiacol and ρ is the density of 4-propylguaiacol. Re is the particle Reynolds number 23 . Sh is the Sherwood number, which is based upon the Frössling correlation 24 .
The chemical reaction which occurs inside (within) the pellets is incorporated into the mass balances with the Reactive Pellet Bed feature in COMSOL ® . This feature has a predefined 1-D extra dimension on the normalised radius of the catalyst pellet particle . The mass balance inside the catalyst ( = / ) pellet is obtained by performing a shell balance across a spherical shell: where r is the catalyst particle radius (dimensionless), D i,eff is the effective diffusion coefficient of chemical species i in the catalyst pores, c i,p is the concentration of chemical species i in the catalyst particle in mol/m 3 . R i,p is the reaction source term (rate of reaction per unit volume of catalyst particle). The effective diffusivities of 4-propylguaiacol and hydrogen into the pores of the catalyst pellet are calculated by relating the diffusion coefficient to either the bulk or Knudsen diffusivity.
is the bulk diffusivity of chemical species , is , the pellet porosity, is the constriction factor and is the tortuosity. Typical values of the constriction factor, the tortuosity, and the pellet porosity are, respectively, = 0.8, = 3.0 and = 0.4 25 . Boundary conditions used were as per the following: at y = l 1 ;

Effects of temperature
The results obtained from the CFD pseudo homogeneous and detailed models were compared to the experimental data. Both models were compared to experimental results for a microreactor operating at a pressure of 300 psig, and temperatures of 250-450 o C. The comparison between the simulated results and the experimental results give an indication towards the validity and robustness of the model. Figure 2 shows the comparison between the modelling and experimental results. The graph depicts the conversion of 4propylguaiacol obtained for a temperature range of 250-450 o C with a good agreement between the experimental and obtained CFD (i.e. pseudo-homogeneous and detailed model) results. The results show that as the temperature increases, the conversion of 4-propylguaiacol also increases. For temperatures below 400 o C, there is a very good validation between the results; however, there is a slight deviation in the results at temperatures, which exceed 400 o C. This could be due to the fact that secondary side reactions are occurring within the microreactor 26 , these reactions are not considered within the CFD modelling due to the lack of reaction kinetic data available. The predominant products from the HDO of 4-propylgauaicol

Effects of hydrogen partial pressure
The effect of hydrogen partial pressure on the conversion of 4propylguiacol was investigated using the models developed in this work and were compared with the results from literature. The range of pressures used for the study varied from 240-480 psig at a constant reactor temperature of 400 o C with a constant residence time. The results in Figure 3 show that as the hydrogen partial pressure increases, the conversion of 4propylguaiacol also increases. However, at pressures of approximately 400 psig, the conversion remains relatively constant. This is because at higher temperatures there is a maximum adsorption of hydrogen on the surface of the catalyst resulting in a stable conversion. The comparison in results between the CFD modelling and experimental show a very good agreement.

Effects of liquid flow rate of 4-propylguaiacol
The effect of liquid flow rate on the conversion of 4propylguaiacol was also studied using the proposed models. The liquid flow rate of 4-propylguaiacol varied between 0.03-0.15 mL min-1, and the remaining operating parameters were kept constant. The studies were performed at a temperature and pressure of 400 o C and 300 psig respectively. Figure 4 shows that as the liquid flow rate increases, the conversion of 4propylguaiacol decreases. The reacting fluids spend a shorter time within the microreactor as the liquid flow rate decreases, hence the conversion declines. Both models predict a decreasing conversion profile with the liquid flow rate which is consistent with the dependence of conversion of a reactant upon space velocity. 4-propylguaiacol is only a reactant in the reaction network presented by the experimentalists 3 , which cannot justify the existence of a maximum in the profile. This is also consistent with the findings and conclusions of the experimentalists 3 , who also describe their experimental conversion profile as decreasing. We believe that the maximum in the experimental profile is an artefactual one and caused by small experimental errors. The deviation of the experimental conversion values is within in the experimental error range.
As the results obtained from the models are in agreement with the experimental findings 3 , they demonstrate a good validation of our models. Figure 5 shows the concentration profiles of 4-propylguaiacol in the bulk fluid phase across the transverse direction at different axial positions. It can be observed that the concentration of 4propylguaiacol is highest at the initial axial positions of the microreactor, this is where the 4-propylguaiacol first encounters the hydrogen at the interface (at x=0; y=0). As the axial position of the microreactor progresses, the concentration decreases until eventually becoming stable. This is because, as the concentration of hydrogen increases along the axial direction of the microreactor, it will encounter greater concentrations of the 4-propylguaiacol, which results in a stable reactant conversion.

Rate analysis
Please do not adjust margins Please do not adjust margins The effect of hydrogen concentration on the rate of disappearance of 4PG was investigated using the detailed model. The hydrogen concentration varied between 0-0.15 mol/L, and the concentration of 4-propylguaiacol remained constant at 1.1 mol/L. Three different temperatures of 350 o C, 300 o C and 250 o C were investigated to obtain the CFD modelling results. All other parameters were kept constant, and the conversion of the reactants was limited to a maximum of 10% in order to solely determine the reaction rates founded on the initial concentrations of the reacting fluids. Figure 6 shows the results obtained from the study. It can be observed that as the concentration of hydrogen increases, the rate of 4propylguaiacol also increases.
The effect of 4-propylguaiacol concentration on the rate of disappearance of 4-propylguaiacol was then investigated. The 4-propylguaiacol concentration varied from 0.1-2.1 mol/L, at a constant hydrogen pressure of 208 psig. The three temperatures used for the study were 350 o C, 300 o C and 250 o C. Figure 7 shows the results obtained. It can be observed that at lower concentrations of 4-propylguaiacol the rate of 4propylguaiacol appears to increase. However, at higher concentrations, the rate remains relatively constant.

Internal and External mass transfer limitations
The heterogeneous detailed model incorporates the solid catalyst, and so it can ascertain the internal and external mass transfer resistances within the microreactor. The model can demonstrate which parameters lead to the reaction being diffusion limited or surface-reaction-limited. Using the CFD models to study the internal mass transfer limitations can allow the determination of which factors enhance or diminish the mass transfer rates (affecting the apparent rate of reaction); therefore, an understanding of the optimisation of the HDO reaction can be achieved.

Internal mass transfer limitations
The concentration profiles of 4-propylguaiacol are presented in figure 8. The profiles were obtained at x = 90 mm, and varying reactor heights of y = 0.7; 0.5 and 0.2 mm, for a catalyst particle size of 75-150 µm. The internal mass transfer resistance is responsible for the concentration gradient within the catalyst particle. It can be deduced that the difference in concentration from the surface of the particle, to within the particle, is less than 5%. In addition, the model was simulated with constant reactor properties, and catalyst particle sizes which were halved and quartered. It was observed that there was no significant difference in the conversion of 4-propylguaiacol (less than 2%). The internal mass transfer limitation was further validated by calculating the Thiele modulus ( ) for a particle size ranging Φ from 70-160 µm by assuming the reaction to be pseudo-firstorder with respect to 4-propylguaiacol and hydrogen 25 .
where ρ p is the density of the catalyst pellet. The Thiele modulus was found to be approximately 0.5. For this value of the Thiele modulus, the effectiveness factor is found to be 1; which suggests that the reaction is surface-reaction-limited (overall rate of reaction is equal to the rate of reaction obtained from within the catalyst particle maintaining the same conditions as the particle surface). Therefore, from this study, it can be concluded that there is negligible internal mass transfer resistance. This conclusion coincides with that obtained from the experimental results. The size of the catalyst particle was doubled and quadrupled to introduce pore diffusion limitations. Figure 9 shows the concentration profiles obtained at x = 90 mm, and varying reactor heights of y = 0.7; 0.5 and 0.2 mm. The introduction of the intraparticle transport resistances is responsible for the steeper concentration profile. This makes microreactors desirable due to their shorter lengths, which allows the use of smaller catalyst particles, which may not usually be feasible in conventional macroscopic reactors. The comprehensive heterogeneous model was used to obtain the internal effectiveness factor. This provides an indication of the relative importance of diffusion and reaction limitations.
The effectiveness factor is regarded as the ratio between the actual overall rate of reaction, and the rate of reaction that would occur if the interior surface of the pellet were exposed to the external surface conditions 25 . The effectiveness factor for a spherical catalyst particle following a first-order reaction can be obtained as: Figure 10 shows a plot of the effectiveness factor against the Thiele Modulus using the detailed model. For the catalyst pellet sizes between 75-150 µm, the Thiele Modulus is approximately less than 1, this corresponds to an effectiveness factor of unity which suggests negligible internal mass transfer resistances. However, increasing the particle sizes results in a decrease of the effectiveness factor, suggesting that the reaction is becoming diffusion-limited within the pellet.  Figure 11 shows a plot of the conversion of 4-propylguaiacol and the liquid flow rate for the pseudo-homogeneous model, detailed model and the detailed model with a catalyst size of 220 µm. It can be observed that the reactant conversion decreases with the larger catalyst sizes. From the modelling results, the overall rate of reaction could be enhanced by decreasing the catalyst particle size; increasing the internal surface area; increasing the temperature; and increasing the concentration. A study of comparison between the pseudo-homogeneous and detailed model was conducted. Figure 12 shows the conversion of 4-propylguaiacol against the liquid flow rate for the detailed model, comprising of catalyst particle sizes ranging between 220-500 µm, and the pseudo-homogeneous model. The results show that as the pellet sizes increase, the reactant conversion decreases because of the diffusion limitations within the particle. The reaction rates obtained from the pseudohomogeneous model were multiplied with the effectiveness factor to obtain the results shown in figure 12. There is a good agreement between the results which demonstrates the validity and robustness between the two models. As a result, either model can be utilised to demonstrate the catalytic hydrodeoxygenation of bio-oil to produce desirable results.
The external mass transfer resistance was then investigated. The hydrodeoxygenation reaction involves the diffusion and mass transfer of hydrogen gas into the 4-propylguaiacol liquid phase, and subsequently diffusion through the liquid phase to the proximate area of the catalyst particle. In order to investigate the resistance to the diffusion across the boundary layer, the concentration surrounding the catalyst particle should be compared to the concentration on the surface of the particle. Figure 13 shows the bulk concentrations of 4propylguaiacol compared to the concentrations obtained from the surface of the catalyst particle. The concentrations of 4propylguaiacol in the bulk are found to be less than 2% when compared to the concentrations on the surface of the catalyst particle. It can be deduced that there is a negligible resistance to mass transfer, and this agrees with the experimental results obtained from Joshi and Lawal (2013).

Conclusions
The modelling results obtained from the CFD study for the catalytic HDO of 4-propylguaiacol to 4-propylphenol have demonstrated a sound validation with the experimental data obtained from literature 3 . The slight deviations in results at temperatures greater than 350 o C could owe to the fact that the results closely follow the % yield of other major reaction products. The secondary side reactions are not considered in this work (due to the lack of kinetic data available). Therefore, the models developed could have slight limitations at higher temperatures because of the previously mentioned reasons; however, this would be overcome by the development of more complex kinetic models which account for the side reactions interacting with the HDO of 4-propylguaiacol to 4-propylphenol. Further model validations were demonstrated by the effects of pressure and liquid flow rate illustrate a good validation with the experimental data. It can be observed that for this specific reaction, the conversion of 4-propylguaiacol increased with increasing temperature and pressure; however, increasing the liquid flow rate appears to decrease the reactant conversion.
The detailed model further allowed the investigation of the concentrations of hydrogen and 4-propylguaiacol on the rate of 4-propylguaiacol. The model results demonstrated that as the concentration of these reactants increased, the rate of disappearance of 4-propylguaiacol also increased. In addition, the detailed model gave rise to the characterisation and evaluation of the reaction-coupled transport phenomena taking place within the catalyst bed in the microreactor. The internal and external mass transfer limitations were investigated by obtaining concentration profiles within the individual catalyst particle. It was concluded that there were negligible internal and external mass transfer resistances, which agreed with the experimental results. In addition, the study of pore diffusion limitations suggested a good agreement between the pseudohomogeneous and detailed heterogeneous models. It can be observed that both models performed similarly when being compared to the experimental data. This indicates the validity and robustness between the models; hence, either model has the ability to predict the catalytic HDO of bio-oil in microreactors. The heterogeneous model has allowed the investigation of pore diffusion limitations by predicting a range of values of the Thiele Modulus at which this occurs. It has shown how this affects the reaction; and gives rise to the study of particle fluid transport phenomena, which aids the understanding of internal and external mass transfer resistances. Performing numerical simulation studies is valuable as it provides an understanding of parameter optimisation and predicting the HDO of other various compounds present in biooil, which would be time consuming and costly to do on an experimental basis.
As the HDO reaction of lignin-derived compounds has not been studied extensively, this model provides a basis to predict and enhance the comprehensive understanding of the HDO reaction in various other lignin-derived compounds in bio-oil. Furthermore, microreactors have demonstrated various benefits 11,27,28 compared to conventional reactors, and so future research should be directed towards investigating the scalability of microreactors to be used on an industrial scale.

Conflicts of interest
There are no conflicts to declare.