Jason Williams ,Hussameldin Ibrahim,*,Nima Karimi ,Kelvin Tsun Wai Ng

1 Clean Energy Technologies Research Institute,Process Systems Engineering,Faculty of Engineering and Applied Science,University of Regina,Regina SK S4S 0A2,Canada

2 Environmental Systems Engineering,Faculty of Engineering and Applied Science,University of Regina,Regina SK S4S 0A2,Canada

Keywords:Hydrogen,crude glycerol Autothermal reforming Numerical analysis Fixed-bed reactor

ABSTRACT A mathematical model for the catalytic autothermal reforming(ATR)reaction of synthetic crude glycerol to hydrogen in a fixed bed tubular reactor(FBTR)and over an in-house developed metal oxide catalyst is presented in this work.The heterogeneous model equations account for a two-phase system of solid catalyst and bulk feed gas.Also,the ATR of crude glycerol reaction scheme and intrinsic kinetic rate model over an active,selective,and stable nickel-based catalyst were integrated in the developed model.Also,the model was validated using experimental data generated in our labs for the ATR of synthetic crude glycerol.The modelling results adequately described the detailed gas product composition and distribution,temperature profiles,and conversion propagation in the axial direction of the fixed bed reactor over a wide range of reaction temperature (773-923 K) and mass-time (12.71-158.23 g cat·min·(mol C)-1).The crude glycerol conversion predicted with the model showing a close resemblance to those obtained experimentally with an average absolute deviation (AAD) of less than 8%.The maximum crude glycerol conversion and hydrogen yield were found to be 92% and 3 mol hydrogen/mol crude glycerol,respectively.Also,the gas product concentration profile in the reactor was adequately described(90%)accuracy with a hydrogen concentration of 39% (volume).

1.Introduction

The provision of sustainable and affordable energy solutions for the world’s growing population poses a major challenge for governments as it is estimated that the world’s population will reach approximately 9 billion by the year 2050[1].Along with the growing population,energy demands will exponentially increase with the need for clean,sustainable sources.Additionally,the outlook of hydrogen production indicates a positive global impact on the world [2].Being the most abundant element on Earth,hydrogen has been identified as a potential energy carrier,most notably for its combustibility releasing a considerably high amount of energy per unit mass.Apart from its industrial application,there are potential pathways for hydrogen utilization in residential homes and as a transportation fuel.Douet al.has discussed the use of various composite materials for storage hydrogen and studied the sorption kinetics to investigate its application during syngas production [3].To date,the most commercialized industrial process for hydrogen production is steam reforming,and even with its drawback of high energy consumption,it boasts high H2yields.However,different technologies have been researched over the years along with different feed sources.Recent work by Douet al,discussed the issues and challenges of thermal conversion using a biomass feedstock [4].Different pathways were explored that could increase the efficiency of the conversion process making it commercially viable,however,limitations exist based on catalyst selection and varying biomass sources.Crude glycerol,a byproduct of biodiesel production,provides an alternative,especially when there is a glut in the market,from methane (CH4),which is the most common feedstock used in H2production [5,6].

Autothermal reforming(ATR)is a combination of steam reforming and partial oxidation,and although steam reforming is the most common method for industrial hydrogen production,an autothermal reaction where steam reforming and partial oxidation occur simultaneously should require no external energy [7].Studies associated with the ATR of methane done by Dauenhaueret al.[8]showed that high steam-to-carbon(SC)ratios resulted in favorable hydrogen selectivity and yield.Other studies showed that with reactions occurring at temperatures in the range of 773 to 1073 K,and under atmospheric pressure,hydrogen yields were higher than those observed with steam reforming.Furthermore,these studies were all predicated on the selection of the appropriate catalyst.An unsteady state model for the ATR or methane was developed by Chanet al.[9],incorporating effectiveness factors to account for diffusion into the catalyst pores.Tanget al.[10] modelled both the SMR and ATR of methane,developing a twodimensional steady-state heterogeneous model,comparing the two methods with the model focusing on varying parameters to achieve full conversion.

This work presents a model developed to simulate the conversion of CG to H2viaan ATR process developed by Abdul Ghani[11].Abdul Ghani’s experimental work included evaluation of the catalytic performance of different compositions of a Ni catalyst over a CeZrM (M=Ca,Mg,Gd) support.With Ca,the highest CG conversion and H2selectivity were observed.The kinetic experiments were performed under conditions to eliminate resistance to heat and mass transfer.After initial experiments,it proves beneficial to transition to a simulation environment,where results can be reproduced at a lower cost and faster rate.A heterogeneous model was developed for this study,which can prove to be more comprehensive,although more computationally expensive than a pseudo-homogeneous model [12].

2.Materials and Methods

2.1.Background

The physical experiments for the ATR of synthetic crude glycerol were done under the portfolio of the Advanced Green Energy Systems Research Group at the University of Regina [11].The experimental setup is shown in Fig.1,with steam and oxygen added to the reactor in different ratios to the amount of synthetic crude glycerol.Synthetic crude glycerol was prepared by mixing components that would be found in crude glycerol generated from biodiesel production,resulting in a mixture with an average molecular weight of 69.08 g·mol-1and an average molecular formula of C2.5H7O2[11].The reactor was enclosed in an electric furnace used to heat the catalyst bed 45 mm in length and 12.7 mm in diameter.After the reactants were passed through the cylindrical reactor packed with spherical catalyst pellets,the outlet gas was passed through a condenser,removing the water,where the remaining gas was analyzed with a gas chromatograph (GC).Experiments were performed under plug flow criterion conditions to eliminate heat and mass transfer limitations [11].

2.2.ATR reaction and kinetic model

Kinetic rate model development for the ATR of synthetic crude glycerol was detailed in[11].Given the different components used to make the mixture,an empirical power-law rate model was suggested for an overall synthetic CG ATR reaction and given as:

This power-law model developed from the experimentation was chosen over generating a mechanistic model due to several complex reactions associated with crude hydrocarbon reforming.From Abdul Ghani’s work,a 95% predictability was achieved [11].

2.3.Model assumptions

The assumptions made in the derivation of the numerical models are as follows:

Fig.1.Schematic diagram of the experimental setup [11].

(1) The continuity equations being used for this model are implemented for cylindrical coordinates,neglecting angular variations,making it two-dimensional.As shown in Fig.2,the reactor tube is represented as a rectangular domain.

(2)The feed mixture entering the reactor is at an elevated temperature (greater than or equal to 773 K) and atmospheric pressure.Under these conditions of temperature and pressure,intermolecular interactions can be neglected,and the feed mixture is assumed to behave as an ideal gas [14,15].Density and heat capacity are assumed to be constant throughout the reactor.

(3)The axial velocity is assumed to be constant and equal to the inlet velocity.

(4)Under the plug flow criterion,dispersion in the radial direction is dominant over convective flux,making radial convection terms for both heat and mass negligible.

2.4.Numerical modelling

The fixed bed tubular reactor(FBTR)model was defined by a set of mass and energy balance partial differential equations (PDEs)with the domain for the model is depicted as shown in Fig.2.

The modelling equations are described as follows [16]:

Bulk Fluid Phase

Mass Balance:

Energy Balance:

For the solid,spherical catalyst phase,the set of equations can be given as follows:

Balance:

Energy Balance:

The initial and boundary conditions are applied for Eqs.(2)-(6)given the modelling domain in Fig.2 are given as:

Fig.2.Schematic of the cross-section of the packed bed reactor [13].

2.5.Numerical solutions

The equations were solved with COMSOL Multiphysics?version 5.2A.This new version of the software has a Reactive Pellet submodule,which makes developing the heterogeneous model easier.The user can enter the physical properties for the pellet and the software platform generates an extended dimension giving a one dimensional (1-D) radial profile of the pellet,as shown in Fig.2.

The Heat Transfer in Fluids module,which is commonly used to generate profiles does not incorporate this feature and assumes that the temperature of the bulk fluid is in equilibrium with the pellet.Allain [17] modelled a 2-D heterogeneous model utilizing the General Extrusion coupling feature,generating temperature profiles for a solid pellet.Given the continuity equations are of the same form,each with a convective and conductive term,a simpler way of generating temperature profiles proposed in this work was to utilize the Transport of Diluted Species module for this heat transfer problem.

3.Results and Discussion

3.1.Model input data

The input values for the model are given in Table 1.

3.2.Model validation

The experimental and predicted conversions of synthetic crude glycerol (CG) are shown in the parity plot of Fig.3 to validate the model’s predictive efficiency.The conversion in the reactor was determined by the following formula:

The average absolute deviation (AAD) of the heterogeneous model was very good at 7.56%.The higher absolute deviations(ADS) were observed in the runs at a lower temperature(T=773 K),which can be attributed to side reactions of components in the synthetic crude glycerol not accounted for (see Table 2).

The outlet concentrations from the model were compared with those provided by the GC analysis.The mole fraction values are given in Table 3 and an acceptable AAD of 10.84% is obtained.A considerably lower methane volume concentration (5.82%) was obtained in the model which may also be due to side reactions for components that are not accounted for in the simplified stoichiometric equation used for the model.

Table 1Operating conditions and parameters used for the FBTR numerical models

Table 2Experimental and prediction conversion results

Table 3Outlet concentration profile of the FBTR at a feed temperature of 923 K and W/FA0 of 127.42 g cat·min·(mol C)-1

3.3.Comparison of the pseudo-homogenous and heterogeneous model

Fig.4 shows the fluid temperature profile for the heterogeneous model,as well as a previously built pseudo-homogeneous model developed by Afabor [13] for the ATR of synthetic crude glycerol.As shown,the inclusion of governing equations for the solid catalyst shows a reduction in the temperature of the reaction zone.The heterogeneous model has a reaction zone temperature ofapproximately 5 K less than that of the pseudo-homogeneous model.

Fig.3.Comparison of measured and predicted crude glycerol conversion for temperature and mass-time ranges of 773-923 K and W/FA0 of 12.71-158.23 g cat·min·(mol C)-1,respectively.

Fig.4.Bulk fluid temperature profile along length of reactor at 873 K and W/FA0 of 12.71 g cat·min·(mol C)-1.

For the given concentration profile in Fig.5,the differences in concentration profile between the models are negligible.Similar results were reported in simulations done by Iordanidi [21] and Routet al.[22],indicating that the reaction may be more sensitive to temperature than to the concentration of the reactants.It can also be observed that the inclusion of the governing equations for the catalyst particle into the model will highlight the existence of inter-phase resistance.

Fig.5.Bulk fluid CG concentration profile along length of reactor at 873 K and W/FA0 of 12.71 g cat·min·(mol C)-1.

When developing the kinetic model,the intent is to remove interphase limitations of heat and mass and this was done by using a smaller pellet size of 0.8 mm.As shown in Fig.6,simulations done based on a smaller pellet size gave a negligible difference between the temperatures of the two phases.For a larger pellet size of 3.8 mm,a larger difference is observed between the temperature of the two phases,where a 2 mm lag in the position of the peak temperature for the fluid phase is observed.This shows that there is a greater interphase resistance with the increase in particle size.The larger pellet size can be favorable when considering the cost of an industrial operation because there is a reduction in pressure drop.However,this may reduce the amount of conversion of the reactants and product yield.

3.4.Effect of temperature and W/FA0 on CG conversion and H2 yield

Fig.6.Temperature profile along the length of the reactor with spherical pellet diameters of(a)0.8 mm and(b)3.8 mm;at a feed temperature of 873 K and W/FA0 of 76.47 g cat·min·(mol C)-1.

Fig.7.Effect of temperature of H2 yield at different W/FA0.

When operating an FBTR,it is important to determine the optimum conditions for maximizing the yield of desired products,minimizing coking within the reactor,which shortens the ‘life’ of the catalyst[23].The mass-time,W/FA0is commonly adjusted to affect the conversion of reactants and hydrogen yield in the reactor.Mass-time in the experiments performed elsewhere [11] was varied by changing the mass of the catalyst while keeping the feed flowrate constant.As shown in,Fig.7,the predicted H2yield at the reactor output increases with increasing temperature,similar observation was made in a thermodynamic analysis of the ATR reaction [24].The simulation showed that the change in H2yield provided diminishing returns as the mass-time or temperature increases.With a feed temperature of 773 K,the rate of conversion shown in Fig.8 decreases significantly after 25 mm into the catalyst bed,and with a higher temperature,this rate would decrease at a shorter distance which may lead us to conclude that there are limits to increasing the catalyst bed length depending on the optimum conditions for achieving a significant H2yield.

Fig.8.Crude glycerol conversion along the length of the reactor with a feed temperature of 773 K at a W/FA0 range from 12.71-101.94 g cat·min·(mol C)-1.

4.Conclusions

A 2-D heterogeneous model was developed for the autothermal reforming of synthetic crude glycerol to hydrogen in this study using FEM in COMSOL Multiphysics.The model was validated against measured experimental data with an acceptable AAD of synthetic crude glycerol conversion of less than 8%.The simulations show that the increase in mass-time and temperature can increase H2yield,greater than 40% in some cases depending on the temperature.Overall,the model proved useful in predicting temperature and concentration variations in the reactor for the ATR process.

Nomenclature

Accross-sectional area of reactor tube,m2

avparticle surface area to volume ratio,m-1

Ci/Ci0concentration of species i,mol·m-3

Cpfspecific heat capacity of bulk fluid,kJ·J-1·K-1

DABbinary diffusivity between two species,m2·s-1

Dereffective radial diffusivity in the gas phase,m2·s-1

Deseffective diffusivity in the solid phase,m2·s-1

Dezeffective axial diffusivity in the gas phase,m2·s-1

Dimdiffusivity of a species in a mixture,m2·s-1

Dtreactor internal diameter,m

dpdiameter of catalyst particle,mm

EAactivation energy,J·mol-1

FCG0inlet molar flowrate of crude glycerol,mol·s-1

Gmass flux,kg.m-2·hr-1

hfseffective heat transfer coefficient between solid and gas,W·m-2·K-1

Hrxheat of reaction,kJ·mol-1

JdColburn-Chilton j factor

kfseffective mass transfer coefficient between solid and gas,m·s-1

Llength of reactor,mm

Maveverage molecular weight,kg·mol-1

Nmolar flux,mol·s-1

Qvolumetric flowrate,m3·s-1

Rgideal gas constant,J·mol-1·K-1

Rradial distance,mm

rC2.5H7O2reaction rate,mol·m-3·s-1

rpradius of catalyst particle,mm

T/To/Twtemperature,K

Utwoverall heat transfer coefficient,W·m-2·K-1

uzaxial superficial velocity,m·s-1

Vstoichiometric coefficient

W/FA0weight-time,g cat·min·(mol C)-1

zaxial distance,mm

εbbulk fluid porosity

εssolid pellet porosity

λer/feffective radial conductivity,W·m-1·K-1

λeseffective solid phase conductivity,W·m-1·K-1

λez/feffective axial conductivity,W·m-1·K-1

μ viscosity,Pa·s

ρbcatalyst bulk density,kg·m-3

ρggas density,kg·m-3

ρssolid density,kg·m-3

τ tortuosity

Subscripts and superscripts

f fluid phase

ps particle surface

s solid phase

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.

Acknowledgements

The authors are grateful for the financial support provided by the Natural Science and Engineering Research Council of Canada(NSERC) and the Canada Foundation for Innovation (CFI).The authors would like to also acknowledge the Clean Energy Technologies Research Institute(CETRI)for providing access to the Simulation Laboratory to conduct this work.The views expressed herein are those of the writers and not necessarily those of our research and funding partners.

Appendix A

The density of the mixture can be calculated as [15]:

whereMaveis average molecular weight of the feed mixture in kg·mol-1;Rgis the ideal gas constant in J·mol-1·K-1.

For effective radial and axial diffusivity,the following correlations can be used [16]:

Qis the volumetric flow rate into the reactorin m3·s-1.NREis the Reynolds number given by:

Diffusion into the catalyst is governed by different phenomena than the diffusivity of the bulk fluid in a heterogeneous system.Equimolar counter-diffusion takes place in the catalyst dependent upon the pore geometry.The effective diffusivity for the solid catalystDesiis described by the following equation [23]:

whereDABis the ordinary bulk diffusivity in a binary gaseous mixture in m2·s-1,τ is tortuosity path of the particle assumed to be 2.The Fuller-Schettler-Giddings correlation given in Eq.(A6) [22]:

With the binary diffusivity between the two species in a particle,the diffusivity of a species within a mixture is given as [25]:

Mis the Molecular weight of the componentin g·mol-1;Vis the diffusion volume;Nis the molar flux mol·s-1.

The particle surface area to volume ratio,av(m-1) [25]:

Effective mass transfer coefficient between solid and gas,kfs(m·s-1) [26]:

Effective heat transfer coefficient between solid and gas,hfs(W·m-2·K-1) [26]: