Numerical solution of nonlinear and non-isothermal general rate model of reactive liquid chromatography

Abstract

A nonlinear general rate model (GRM) of liquid chromatography is formulated to analyze the influence of temperature variations on the dynamics of multi-component mixtures in a thermally insulated liquid chromatographic reactor. The mathematical model is formed by a system of nonlinear convection–diffusion reaction partial differential equations (PDEs) coupled with nonlinear algebraic equations for reactions and isotherms. The model equations are solved numerically by applying a semi-discrete high-resolution finite volume scheme (HR-FVS). Several numerical case studies are conducted for two different types of reactions to demonstrate the influence of heat transfer on the retention time, separation, and reaction. It was found that the enthalpies of adsorption and reaction significantly influence the reactor performance. The ratio of density time heat capacity of solid and liquid phases significantly influences the magnitude and velocity of concentration and thermal waves. The results obtained could be very helpful for further developments in non-isothermal reactive chromatography and provide a deeper insight into the sensitivity of chromatographic reactor operating under non-isothermal conditions.

GRAPHICAL ABSTRACT

Keywords:

• Chromatographic reactor
• LC
• mass and heat transfer
• non-isothermal operation
• TLC

Introduction

High-performance liquid chromatography is a selective separation technique which is based on differences in the distribution of chemical species between the phases of the separation column. The technique is mainly used for the separation of components from complex mixtures, such as fine chemicals, pharmaceutical, food additives, and biological products. In recent years, reactive chromatography has gained ample interest in chemical engineering research. Chromatographic reactor is a multi-functional unit in which chemical reaction and chromatographic separation lead to an integrated process that improves the conversion of reactants and purity of products. In chemical engineering, the reactive chromatography is available as an invaluable laboratory tool and as an industrial unit for multiple reasons.^[1] First, continuous removal of at least one of the reaction products for an equilibrium limited reactions shift the equilibrium in a direction that reduces by-product formation and increases conversion. Second, when a catalyst is used as packing material in the column, it stays inside the column and its removal from the product stream is not required. Third, complex mixtures can be separated with great precision using reactive chromatography, such as separation of proteins mixture or production of biodiesel. In reactive chromatography, it is often easier to meet the purity requirements of the products as compared to other separation methods. Thus, reactive chromatography can be utilized to obtain delicate products and harsh process conditions can be easily avoided.

A reaction inside the chromatographic reactor can be catalyzed homogeneously or heterogeneously. In the homogeneously catalyzed reaction, a catalyst for the reaction is included in the mobile phase. Thus, removal of catalyst from the product is needed at the end of the process. On the other hand, in the heterogeneously catalyzed reaction, the stationary phase is used simultaneously as a catalyst for the reactant and as an adsorbent for the separation of mixture components. Reactive chromatography can be applied to several classes of reactions. Some typical applications include esterification by ion-exchange resins^[2–7] or by immobilized enzymes,^[8] transesterifications,^[9,10] hydrogenation,^[11] and sugar involving reactions.^[12,13]

Elution order of components and the type of reaction has a strong influence on the reactive chromatography. Mostly equilibrium limited reaction of the form $A\hspace{0.17em}\rightleftarrows \hspace{0.17em}B,A\hspace{0.17em}\rightleftarrows \hspace{0.17em}B+C$ and $A+B\hspace{0.17em}\rightleftarrows \hspace{0.17em}B+C$ are extensively investigated in the literature for isothermal case.^[14–19] The elution order constraint is important to gain high purity products. In consecutive or parallel reactions, the formation of by-products leads to a reduction in product purities.^[20] To understand the basic concept of chromatographic batch reactors, consider a single chromatographic column in which an equilibrium limited reversible reaction of the form $A\hspace{0.17em}\rightleftarrows \hspace{0.17em}B+C$ takes place. In such a reaction process, the reactant A is injected as a rectangular pulse into the chromatographic column constituting an adsorbent of lower affinity toward the reactant A and the product C and high affinity toward the product B. The reactant A reacts with the surface of catalyst to form the product B and C. Both product species will move at different velocities through the reactor, with B stays behind the reaction front as being retained strongly than A and C. Pure fractions of products could be collected separately on the basis of their elution time. Such separation facilitates the reversible reaction to go beyond the thermodynamic equilibrium by continuously separating the products from the reaction zone. Therefore, high purity product could be withdrawn from these systems due to their separation from the reactants.^[21–23]

Thermal effects such as heat generated by chemical reaction, heat of adsorption, and mixing are mostly overlooked in the analysis of liquid chromatographic reactors. The available experimental, theoretical, and numerical investigation in the literature is limited to isothermal conditions. A few contributions regarding thermal effects in liquid chromatography are also available in the literature.^[24–29] However, thermal effects have been discussed extensively for gas chromatography.^[30–34] Some authors have discussed the effects of temperature on the constant equilibrium under constant inlet pressure.

Mathematical modeling is widely used as an important tool for studying the chromatographic process. In this regard, a number of mathematical models with different complexity have been developed in the literature, which provides detailed information about the mass transfer and partitioning process inside the chromatographic column. These models include the equilibrium dispersive model (EDM), the non-equilibrium lumped kinetic model (LKM), and the general rate model (GRM).^[35–38] Analytical solutions are possible under linear adsorption conditions only. In the case of non-isothermal reactive chromatography, we have to linearize the adsorption isotherm and reaction term for solving the model equations analytically.^[19] Thus, the same nonlinear and non-isothermal reaction behavior presented in this article cannot be accommodated in the analytical framework. The analytical results presented for linearized isotherms and reactions can only be used to simulate chromatography processes considering small volumes of the injected or diluted samples. Such solutions are only valid and meaningful for small values of the enthalpy of adsorption and they give over-predicted solutions for larger values of the enthalpy of adsorption.

The aim of this work is to demonstrate the significant influence of thermal effects on separation and conversion and to quantify the magnitude of thermal effects on the performance of non-isothermal chromatographic reactor. A three-component non-isothermal GRM of reactive liquid chromatography with nonlinear heterogeneous reaction term is considered along with nonlinear adsorption conditions to simulate the process of chromatographic reactor. The reactive non-isothermal GRM can be described by a coupled system of nonlinear convection-diffusion reaction-type partial differential equations (PDEs) for mass and energy balances coupled with algebraic equations describing the adsorption isotherms and reaction. The high-resolution finite volume scheme (HR-FVS) of Koren is recommended to solve the model equations. Several test problems of practical applicability are considered to explain the coupling between the concentration and temperature profiles in the chromatographic reactor. The significant effect of thermal gradient on separation and conversion is demonstrated and the influence of key parameters is analyzed. It is observed that analogous to gas chromatography, density time heat capacity ratio of solid to liquid phase has a significant influence on the reactor performance. The current nonlinear model is more general and flexible for the simulation of reactive liquid chromatography process considering a variety of samples.

This article is further arranged in the following manner. In Section 2, a three-component GRM is formulated for non-isothermal liquid chromatographic reactor. In Section 3, a semi-discrete HR-FVS is derived and implemented for solving the model equations. Section 4 presents basic formulation and procedure to check the consistency of the results. Section 5 presents some numerical case studies. Finally, conclusions are drawn in Section 6.

Non-isothermal GRM of reactive chromatography

The GRM of column chromatography incorporates several sorption and transport processes. Molecules in the mobile phase are transported between the chromatographic beads by convection and are dispersed due to eddies and other inhomogeneities of flow. The film mass transfer equilibrates concentration gradient between the particle bulk phase and the stagnant film around the particle surface. The molecules are then diffused in the pores of the particles (mobile phase) and, hence, adsorbed selectively at their inner surface. Many properties of sorption dynamics and equilibria play a major role to separate the individual components from a complex mixture.

In the derivation of the model equations, it is assumed that column is insulated thermally and is packed homogeneously, the mobile phase is incompressible, the volumetric flow rate is constant and the solid phase is acting as a catalyst for heterogeneous reaction, the effect of temperature variation is negligible on the physical properties (e.g. viscosity, density, and heat capacity) and transport coefficients (e.g. axial dispersion and heat-conductivity), and the axial dispersion and the heat-conductivity coefficients are independent of the flow rate.

The governing mass balance equations for the bulk-fluid phase are expressed as (1) $$\frac{\partial c_{b,i}}{\partial t}+u\frac{\partial c_{b,i}}{\partial z}=D_{b,i}\frac{{\partial }^2{c}_{b,i}}{\partial z^2}-\frac{3}{R_p}{F}_b{k}_{\text{eff},i}\left(c_{b,i}-c_{p,i}{|}_{R=R_p}\right),\hspace{1em}i=1,2,\dots ,N_c$$(1) where N_c represents the number of components in the mixture introduced at the column inlet, $c_{b,i}(t,z)$ represents the solute concentration of the ith component in the bulk-fluid phase, u is the interstitial velocity, $D_{b,i}$ denotes the axial dispersion coefficient, the phase ratio $F_b=(1-{ϵ}_b)/{ϵ}_b$ represents the interstitial bulk volume over the total bead volume, ${ϵ}_b$ is the bulk porosity, $k_{\text{eff},i}$ is the effective mass transfer coefficient, R_p is the radius of spherical particles, the factor $\frac{3}{R_p}$ denotes the surface to volume ratio of spherical particles, $c_{p,i}(t,z,r)$ is the i-th component solute concentration in the particle pores, and $(c_{b,i}-c_{p,i}{|}_{R=R_p})$ represents the concentration differences between the extra particular mobile phase in the external film and the intra particular mobile phase at the surface of the particle.

The corresponding differential mass balance for the sorbent particle pores, assuming diffusion and sorption of components within the sorbent pores and reaction in the stationary phase, is expressed as (2) $${ϵ}_p\frac{\partial c_{p,i}}{\partial t}+(1-{ϵ}_p)\frac{\partial q_{p,i}}{\partial t}=D_{\text{eff}}\left[\frac{1}{R^2}\frac{\partial }{\partial R}\left(R^2\frac{\partial c_p}{\partial R}\right)\right]+(1-{ϵ}_p){\nu }_i{r}^{\text{het}}$$(2) where $q_{p,i}=q_{p,i}(t,z,r)$ is the local equilibrium concentration of solute in the stationary phase for i-th component, $D_{\text{eff}}={ϵ}_{\mathrm{p}}{D}_{\mathrm{p}}$ is the effective internal pore diffusivity, $r^{\text{het}}$ denotes the heterogeneous reaction rate in the solid phase and ν_i are the corresponding stoichiometric coefficients of i-th component. In general, the stoichiometric coefficient ν_i is negative for reactants and positive for products.

The corresponding energy balance for an adiabatic chromatographic reactor is given as (3) $$\frac{\partial T_b}{\partial t}+u\frac{\partial T_b}{\partial z}=\frac{{\lambda }_{\text{eff},\mathrm{z}}}{{\rho }^L{c}_p^L}\frac{{\partial }^2{T}_b}{\partial z^2}-\frac{3F_b}{R_p{\rho }^L{c}_p^L}{h}_{\text{eff}}\left(T_b-T_p(R=R_p)\right)$$(3)

In the above equations, T_b is temperature of the bulk fluid, T_p is fluid temperature inside the particle pores, ${\lambda }_{\text{eff},\mathrm{z}}$ represents effective axial heat conductivity coefficient, ρ^L is density per unit volume in the liquid phase, $c_P^L$ is heat capacity for liquid phase and $h_{\text{eff}}$ is effective heat transfer coefficient quantifying the rate of heat exchange between mobile and stationary phases.

An energy balance law for the radial temperature profile inside particles pores is expressed as (4) $$\overline{\rho c_p}\frac{\partial T_p}{\partial t}-(1-{ϵ}_p)\sum _{j=1}^N(-\Delta H_{A,j})\frac{\partial q_{p,i}}{\partial t}={\lambda }_{\text{eff},\mathrm{e}}\left[\frac{1}{R^2}\frac{\partial }{\partial R}\left(R^2\frac{\partial T_p}{\partial R}\right)\right]+(1-{ϵ}_p)(-\Delta H_R)r^{\text{het}}$$(4) here $\Delta H_{A,j}$ represents the enthalpy of adsorption of j-th component, $\Delta H_R$ is the enthalpy of reaction, $r^{\text{het}}$ is the heterogeneous reaction rate, ${\lambda }_{\text{eff},\mathrm{e}}$ denotes the internal heat diffusivity coefficient. Moreover, (5) $$\overline{\rho c_p}={ϵ}_p{\rho }^L{c}_p^L+(1-{ϵ}_p){\rho }^S{c}_p^S$$(5) where ρ^S is the density per unit volume of solid phase and $c_p^S$ is the corresponding heat capacity. The considered density and heat capacity ρ^L, ρ^S, $c_p^L,$ and $c_p^S$ are not depending on temperature.

Consider a simple case of non-reactive, non-dispersive, and equilibrium chromatography, i.e. $c_{b,i}=c_{p,i},$ $q_{b,i}=q_{p,i},$ and all dispersion and reaction terms are zero. Then, the speeds of concentration profiles $u_c^k$ and of the thermal wave u_T inside the column can be approximated from Equations (1)–(4) as^[29] (6) $$u_c^k\cong \frac{u}{1+F\frac{\partial q_{b,k}}{\partial c_{b,k}}} ,\hspace{1em}{u}_T\cong \frac{u}{1+F\frac{{\rho }^{\mathrm{S}}{C}_{\mathrm{p}}^{\mathrm{S}}}{{\rho }^{\mathrm{L}}{C}_{\mathrm{p}}^{\mathrm{L}}}} ,\hspace{1em}k=1,2,\dots ,N$$(6)

It is clear from Equation (6)(6) $$u_c^k\cong \frac{u}{1+F\frac{\partial q_{b,k}}{\partial c_{b,k}}} ,\hspace{1em}{u}_T\cong \frac{u}{1+F\frac{{\rho }^{\mathrm{S}}{C}_{\mathrm{p}}^{\mathrm{S}}}{{\rho }^{\mathrm{L}}{C}_{\mathrm{p}}^{\mathrm{L}}}} ,\hspace{1em}k=1,2,\dots ,N$$(6) that the speed $u_c^k$ is depending on the respective local gradients of adsorption isotherms and u_T is influenced by the ratio of density times heat capacity. For less adsorbed components or for smaller gradients of isotherms $\frac{\partial q_k}{\partial c_k}$ and for a larger ratio of $\frac{{\rho }^{\mathrm{S}}{C}_{\mathrm{p}}^{\mathrm{S}}}{{\rho }^{\mathrm{L}}{C}_{\mathrm{p}}^{\mathrm{L}}},$ the speeds $u_c^k$ of concentration fronts are higher than the speed u_T of thermal wave.

The equilibrium adsorption isotherms are assumed to be nonlinear in terms of the concentrations and temperature and are described by van’t Hoff type relation using the enthalpy of adsorption (7) $$q_{p,i}=\frac{a_i^{\text{ref}}{c}_{p,i}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-\Delta H_{A,i}}{R_g}\left(\frac{1}{T_p}-\frac{1}{T^{\text{ref}}}\right)\right)}{1+\sum _{j=1}^{N_c}{b}_j^{\text{ref}}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-\Delta H_{A,j}}{R_g}\left(\frac{1}{T_p}-\frac{1}{T^{\text{ref}}}\right)\right)c_{p,j}}$$(7) where $a_i^{\text{ref}}$ denotes the Henry’s constant and $b_i^{\text{ref}}$ represents the nonlinearity coefficient of the ith component at the reference temperature, R_g is the gas constant, and $T^{\text{ref}}$ represents the reference temperature. The chemical reaction in the chromatographic reactor could occur either in the liquid phase (i.e., homogenous reaction) or in the particle phase (i.e., heterogeneous reaction) or in both phases. Heterogeneously catalyzed reaction is usually considered for esterification where the same ion exchange resin acts as a catalyst for the reaction and as an absorbent for the separation. In this study, we consider only the heterogeneous (solid phase) reaction. Based on the laws of conservation, the reaction rate of three components model reaction ($A\hspace{0.17em}\rightleftarrows \hspace{0.17em}B+C$) is given as (8) $$r^{\text{het}}=k^{\text{het}}(T_p)\left(q_{p,A}-\frac{q_{p,B}{q}_{p,C}}{K_{eq}^{\text{het}}}\right)$$(8) here $k^{\text{het}}$ represent the heterogeneous forward reaction rate constant and $K_{eq}^{\text{het}}$ represent reaction equilibrium constant, respectively. The Arrhenius equation uses activation energy $E_A^{\text{het}}$ to describe the temperature effect on the chemical reaction rate $k^{\text{het}}(T_s)$ as an exponential function of absolute temperature: (9) $$k^{\text{het}}(T_p)=k^{\text{het}}(T^{\text{ref}})\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-E_A^{\text{het}}}{R_g}\left(\frac{1}{T_p}-\frac{1}{T^{\text{ref}}}\right)\right)$$(9)

To simplify the notations and reduce the number of parameters appearing in the model equations, the following dimensionless quantities are introduced: (10) $$\begin{array}{c}x=\frac{z}{L},\hspace{1em}\tau =\frac{ut}{L},\hspace{1em}r=\frac{R}{R_p},\hspace{1em}P{e}_{b,i}=\frac{Lu}{D_{b,i}},\hspace{1em}P{e}_T=\frac{Lu{\rho }^L{c}_P^L}{{\lambda }_{\text{eff},\mathrm{z}}},\hspace{1em}{B}_{c,i}=\frac{k_{\text{eff}}{R}_p}{D_{\text{eff}}},\hspace{1em}{B}_T=\frac{h_{\text{eff}}{R}_p}{{\lambda }_{\text{eff},\mathrm{e}}},\hspace{1em}{\eta }_{c,i}=\frac{D_{\mathrm{p},\mathrm{i}}L}{R_p^{2}u},\hspace{1em}{\eta }_T=\frac{{\lambda }_{\text{eff},\mathrm{e}}L}{R_p^{2}u{\rho }^L{c}_P^L},\\ \hspace{1em}{\xi }_{c,i}=3B_{c,i}{\eta }_{c,i}{F}_b,\hspace{1em}{\xi }_T=3B_T{\eta }_T{F}_b \end{array}$$(10) here L denotes the column length, while $P{e}_{b,i}$ and Pe_T is the peclet numbers, $B_{c,i}$ and B_T is the Biot numbers for mass and energy, while ${\eta }_{c,i},$ η_T, ${\xi }_{c,i},$ and ξ_T are the dimensionless constants. After using the above dimensionless parameters in the mass and energy balances (c.f. Equations (1)–(4)), we obtain (11) $$\frac{\partial c_{b,i}}{\partial \tau }+\frac{\partial c_{b,i}}{\partial x}=\frac{1}{P{e}_{b,i}}\frac{{\partial }^2{c}_{b,i}}{\partial x^2}-{\xi }_{c,i}\left[c_{b,i}-c_{p,i}(r=1)\right]$$(11) (12) $$\frac{\partial c_{p,i}}{\partial \tau }+F_p\frac{\partial q_{p,i}}{\partial \tau }={\eta }_{c,i}\left[\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial c_{p,i}}{\partial r}\right)\right]+F_p\frac{L}{u}{\nu }_i{r}^{\text{het}}$$(12) (13) $$\frac{\partial T_b}{\partial \tau }+\frac{\partial T_b}{\partial x}=\frac{1}{P{e}_T}\frac{{\partial }^2{T}_b}{\partial x^2}-{\xi }_T\left[T_b-T_p(r=1)\right]$$(13) (14) $$\left(1+F_p\frac{{\rho }^S{c}_p^S}{{\rho }^L{c}_p^L}\right)\frac{\partial T_p}{\partial \tau }-F_p\sum _{j=1}^{N_c}\frac{(-\Delta H_{A,j})}{{\rho }^L{c}_p^L}\frac{\partial q_{p,i}}{\partial \tau }={\eta }_T\left[\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial T_p}{\partial r}\right)\right]+F_p\frac{(-\Delta H_R)}{{\rho }^L{c}_p^L}\frac{L}{u}{r}^{\text{het}}$$(14) where $F_p=\frac{1\hspace{0.17em}-\hspace{0.17em}{ϵ}_p}{{ϵ}_p}.$ Equations (c.f. Equations (11)–(14)) are solved using the appropriate initial and boundary conditions (BCs).

The initial conditions are given as (15) $$c_{b,i}(0,x)=0, T_b(0,x)=T_b^{\text{init}}, 0\le x\le 1$$(15) (16) $$c_{p,i}(0,x,\rho )=0, T_p(0,x,\rho )=T_p^{\text{init}}, 0\le x\le 1, 0\le \rho \le 1, 0\le {\rho }_p\le 1$$(16) here $T_b^{\text{init}}$ and $T_p^{\text{init}}$ represent the initial temperature in the bulk and particle phase of the column, respectively.

BCs are required at the column inlet and the column outlet, at the center of the particles and for the stagnant film. In this study, the following Robin type BCs, also known as Danckwert BCs,^[39] are considered at the column inlet. These BCs are based on the flux conversation which are expressed as (17a) $${-\frac{1}{P{e}_{b,i}}\frac{\partial c_{b,i}}{\partial x}+c_{b,i}|}_{x=0}=\{\begin{array}{c}{c}_{b,i}^{\text{inj}},\hspace{1em}\text{if} \hspace{1em}0\le \tau \le {\tau }^{\text{inj}} \hfill \\ \hspace{1em}0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\tau >{\tau }^{\text{inj}} \hfill \end{array}$$(17a) (17b) $${-\frac{1}{P{e}_T}\frac{\partial T_b}{\partial x}+T_b|}_{x=0}=\{\begin{array}{c}{T}_b^{\text{inj}},\hspace{1em}\text{if} \hspace{1em}0\le \tau \le {\tau }^{\text{inj}} \hfill \\ \hspace{1em}0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\tau >{\tau }^{\text{inj}} \hfill \end{array}$$(17b) here $c_{b,i}^{\text{inj}}$ is concentration and $T_b^{\text{inj}}$ is temperature of the injected sample. In this work, we have taken $T_b^{\text{inj}},T_b^{\text{init}},$ and $T^{\text{ref}}$ the same. At the column outlet (x = 1), the zero Neumann BCs are utilized: (17c) $${\frac{\partial c_{b,i}}{\partial x}|}_{x=1}=0,\hspace{1em}{\frac{\partial T_b}{\partial x}|}_{x=1}=0$$(17c)

For Equations (12)(12) $$\frac{\partial c_{p,i}}{\partial \tau }+F_p\frac{\partial q_{p,i}}{\partial \tau }={\eta }_{c,i}\left[\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial c_{p,i}}{\partial r}\right)\right]+F_p\frac{L}{u}{\nu }_i{r}^{\text{het}}$$(12) and (14)(14) $$\left(1+F_p\frac{{\rho }^S{c}_p^S}{{\rho }^L{c}_p^L}\right)\frac{\partial T_p}{\partial \tau }-F_p\sum _{j=1}^{N_c}\frac{(-\Delta H_{A,j})}{{\rho }^L{c}_p^L}\frac{\partial q_{p,i}}{\partial \tau }={\eta }_T\left[\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial T_p}{\partial r}\right)\right]+F_p\frac{(-\Delta H_R)}{{\rho }^L{c}_p^L}\frac{L}{u}{r}^{\text{het}}$$(14) , radial BCs at r = 0 and r = 1 are expressed as (17d) $${\frac{\partial c_{p,i}}{\partial r}|}_{r=0}=0,\hspace{1em}{\frac{\partial c_{p,i}}{\partial r}|}_{r=1}=B_{c,i}(c_{b,i}-{c_{p,i}|}_{r=1})$$(17d) (17e) $${\frac{\partial T_p}{\partial r}|}_{r=0}=0,\hspace{1em}{\frac{\partial T_p}{\partial r}|}_{r=1}=B_T(T_b-{T_p|}_{r=1})$$(17e)

This completes the derivation of non-isothermal GRM for chromatographic reactors. In the next section, the suggested semi-discrete HR-FVS is applied to solve the model equations.

Numerical scheme

Various numerical procedures have been introduced in the literature for approximating the model equations for chromatographic models.^{[37,38,40–42]} In this article, a semi-discrete high resolution flux-limiting finite volume method of Koren is applied to solve current non-isothermal GRM.^[41,42] A second-order TVD-RK method is used for solving the resulting ODEs system.^[41,42] The order of accuracy of this scheme has been verified analytically and numerically in one of our previous articles.^[42] It was found that this scheme is second-to-third order accurate. Moreover, the scheme has capability to resolve sharp fronts and peaks in the solutions accurately. In this section, a complete derivation of the proposed numerical scheme is presented for a three-component reaction.

In compact form, the above system of equations (c.f. Equations (11)–(14)) can be rewritten as (18) $$\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial \tau }+\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial x}=\mathbf{P}\frac{{\partial }^2{\mathbf{c}}_{\mathbf{b}}}{\partial x^2}-\mathbf{\xi }({\mathbf{c}}_{\mathbf{b}}-{\mathbf{c}}_p{|}_{r=1})$$(18) (19) $$\mathbf{J}\frac{\partial {\mathbf{c}}_p}{\partial \tau }=\frac{\mathbf{\eta }}{{ϵ}_p}\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial {\mathbf{c}}_p}{\partial r}\right)+F_p\frac{L}{u}\mathbf{R}{r}^{\text{het}}$$(19) where (20) $$\begin{array}{c}\begin{array}{cc}\hfill & \begin{array}{c} \mathbf{c}=\left[\begin{array}{c}{c}_{b,1}\\ c_{b,2}\\ c_{b,3}\\ T_b\end{array}\right],\hspace{1em}{\mathbf{c}}_p=\left[\begin{array}{c}{c}_{p,1}\\ c_{p,2}\\ c_{p,3}\\ T_p\end{array}\right],\hspace{1em}\mathbf{P}=\left[\begin{array}{cccc}\frac{1}{P{e}_{b,1}}& 0& 0& 0\\ 0& \frac{1}{P{e}_{b,2}}& 0& 0\\ 0& 0& \frac{1}{P{e}_{b,3}}& 0\\ 0& 0& 0& \frac{1}{P{e}_T}\end{array}\right],\\ \mathbf{\xi }=\left[\begin{array}{cccc}{\xi }_{c,1}& 0& 0& 0\\ 0& {\xi }_{c,2}& 0& 0\\ 0& 0& {\xi }_{c,3}& 0\\ 0& 0& 0& {\xi }_T\end{array}\right], \mathbf{R}=\left[\begin{array}{cccc}{\nu }_{c,1}& 0& 0& 0\\ 0& {\nu }_{c,2}& 0& 0\\ 0& 0& {\nu }_{c,3}& 0\\ 0& 0& 0& \frac{-\Delta H_R}{{\rho }^L{c}_p^L}\end{array}\right], \eta =\left[\begin{array}{cccc}{\eta }_1& 0& 0& 0\\ 0& {\eta }_2& 0& 0\\ 0& 0& {\eta }_3& 0\\ 0& 0& 0& {\eta }_T\end{array}\right],\end{array}\hfill \end{array}\\ \mathbf{J}=\left[\begin{array}{cccc}1+F_p\frac{\partial q_{p,1}}{\partial c_{p,1}}& F_p\frac{\partial q_{p,1}}{\partial c_{p,2}}& F_p\frac{\partial q_{p,1}}{\partial c_{p,3}}& F_p\frac{\partial q_{p,1}}{\partial T_p}\\ F_p\frac{\partial q_{p,2}}{\partial c_{p,1}}& 1+F_p\frac{\partial q_{p,2}}{\partial c_{p,2}}& F_p\frac{\partial q_{p,2}}{\partial c_{p,3}}& F_p\frac{\partial q_{p,2}}{\partial T_p}\\ F_p\frac{\partial q_{p,3}}{\partial c_{p,1}}& F_p\frac{\partial q_{p,3}}{\partial c_{p,2}}& 1+F_p\frac{\partial q_{p,3}}{\partial c_{p,3}}& F_p\frac{\partial q_{p,3}}{\partial T_p}\\ F_p\frac{\Delta H_A}{{\rho }^L{c}_p^L}{\sum }_{j=1}^3\frac{\partial q_{p,j}}{\partial c_{p,1}}& F_p\frac{\Delta H_A}{{\rho }^L{c}_p^L}{\sum }_{j=1}^3\frac{\partial q_{p,j}}{\partial c_{p,2}}& F_p\frac{\Delta H_A}{{\rho }^L{c}_p^L}{\sum }_{j=1}^3\frac{\partial q_{p,j}}{\partial c_{p,2}}& 1+F_p\frac{{\rho }^S{c}_p^S}{{\rho }^L{c}_p^L}+F_p\frac{\Delta H_A}{{\rho }^L{c}_p^L}{\sum }_{j=1}^3\frac{\partial q_{p,j}}{\partial T_p}\end{array}\right].\end{array}$$(20)

The model Equations (18)(18) $$\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial \tau }+\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial x}=\mathbf{P}\frac{{\partial }^2{\mathbf{c}}_{\mathbf{b}}}{\partial x^2}-\mathbf{\xi }({\mathbf{c}}_{\mathbf{b}}-{\mathbf{c}}_p{|}_{r=1})$$(18) and (19)(19) $$\mathbf{J}\frac{\partial {\mathbf{c}}_p}{\partial \tau }=\frac{\mathbf{\eta }}{{ϵ}_p}\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial {\mathbf{c}}_p}{\partial r}\right)+F_p\frac{L}{u}\mathbf{R}{r}^{\text{het}}$$(19) ) constitute a system of nonlinear PDEs. These equations are discretized first in the spatial domain by using finite volume scheme. Domain discretization will result in a system of ordinary differential equations which is also nonlinear.

Domain discretization

Let N_x and N_r are the numbers of discretization points along the x and r-coordinates. The computational domain is taken as $[0,1]\times [0,1]$ which is covered by cells ${\Omega }_{lm}\equiv [x_{l\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}},x_{l\hspace{0.17em}+\hspace{0.17em}\frac{1}{2}}]\times [r_{pm\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}},r_{pm\hspace{0.17em}+\hspace{0.17em}\frac{1}{2}}]$ for $1\le l\le N_x,$ and $1\le m\le N_{r_p}.$ The characteristic coordinate points in the cell Ω_lm are represented by $(x_l,r_{pm}).$ Here, (21) $$(x_{\frac{1}{2}},x_{\frac{1}{2}})=(0,0),\hspace{1em}{x}_l=\frac{x_{l\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}}+x_{l\hspace{0.17em}+\hspace{0.17em}\frac{1}{2}}}{2},\hspace{1em}{r}_{pm}=\frac{r_{pm\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}}+r_{pm\hspace{0.17em}+\hspace{0.17em}\frac{1}{2}}}{2}$$(21) and for the uniform mesh (22) $$\Delta x=x_{l\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}}-x_{l\hspace{0.17em}+\hspace{0.17em}\frac{1}{2}},\hspace{1em}\Delta r_p=r_p{}_{m\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}}-r_p{}_{m\hspace{0.17em}+\hspace{0.17em}\frac{1}{2}}$$(22) Since (23) $$c_b:=c_b(\tau ,x)\hspace{1em}\text{and}\hspace{1em}{c}_p:=c_p(\tau ,x,r_p)$$(23)

Therefore, for $I_l:=[x_{l\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}},x_{l\hspace{0.17em}+\hspace{0.17em}\frac{1}{2}}]$ and ${\Omega }_{lm}:=[x_{l\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}},x_{l\hspace{0.17em}+\hspace{0.17em}\frac{1}{2}}]\times [r_p{}_{m\hspace{0.17em}-\hspace{0.17em}\frac{1}{2}},r_p{}_{m\hspace{0.17em}+\hspace{0.17em}\frac{1}{2}}],$ the averaged values of the cell $c_{b,l}(\tau )$ and $c_{p,l,m}(\tau ),$ expressed at any time τ, are given as (24) $$c_{b,l}=c_{b,l}(\tau )=\frac{1}{\Delta x_l}{\int }_{I_l}{c}_b(\tau ,x)dx$$(24) (25) $$c_{p,l,m}=c_{p,l,m}(\tau )=\frac{1}{\Delta x_l\Delta r_p{}_m}{\int }_{{\Omega }_{lm}}{c}_p(\tau ,x,r_p)d{r}_{p}dx$$(25)

By integrating Equation (18)(18) $$\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial \tau }+\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial x}=\mathbf{P}\frac{{\partial }^2{\mathbf{c}}_{\mathbf{b}}}{\partial x^2}-\mathbf{\xi }({\mathbf{c}}_{\mathbf{b}}-{\mathbf{c}}_p{|}_{r=1})$$(18) over the interval I_l and utilizing Equations (24)(24) $$c_{b,l}=c_{b,l}(\tau )=\frac{1}{\Delta x_l}{\int }_{I_l}{c}_b(\tau ,x)dx$$(24) and (25)(25) $$c_{p,l,m}=c_{p,l,m}(\tau )=\frac{1}{\Delta x_l\Delta r_p{}_m}{\int }_{{\Omega }_{lm}}{c}_p(\tau ,x,r_p)d{r}_{p}dx$$(25) , we obtain (26) $$\frac{d{\mathbf{c}}_{b,l}}{d\tau }=-\frac{{\mathbf{c}}_{b,l+\frac{1}{2}}-{\mathbf{c}}_{b,l-\frac{1}{2}}}{\Delta x}+\frac{\mathbf{P}}{\Delta x}\left[{\left(\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial x}\right)}_{l+\frac{1}{2}}-{\left(\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial x}\right)}_{l-\frac{1}{2}}\right]-\mathbf{\xi }\left({\mathbf{c}}_{b,l}-{\mathbf{c}}_{p,l,N_{r_p}}\right)$$(26) where $l=1,2,\dots ,N_x,$ and the differential term of the axial diffusion part is approximated as (27) $${\left(\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial x}\right)}_{l\pm \frac{1}{2}}=\pm \frac{({\mathbf{c}}_{b,l\pm 1}-{\mathbf{c}}_{b,l})}{\Delta x}$$(27)

Integration of Equation (19)(19) $$\mathbf{J}\frac{\partial {\mathbf{c}}_p}{\partial \tau }=\frac{\mathbf{\eta }}{{ϵ}_p}\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial {\mathbf{c}}_p}{\partial r}\right)+F_p\frac{L}{u}\mathbf{R}{r}^{\text{het}}$$(19) over the interval Ω_lm gives (28) $$\frac{d{\mathbf{c}}_{p,l,m}}{d\tau }\hspace{0.17em}=\hspace{0.17em}{\mathbf{J}}_{l,m}^{-1}\eta \frac{1}{{ϵ}_p{r}_p^2{}_{m+1/2}\Delta r_p}\left[{({\mathbf{c}}_p)}_{l,m+1/2}-{({\mathbf{c}}_p)}_{l,m-1/2}\right]+F_p\frac{L}{u}{\mathbf{J}}_{l,m}^{-1}\mathbf{R}{r}_{l,m}^{\text{het}}$$(28) where (29) $${({\mathbf{c}}_p)}_{l,m\pm 1/2}=\max \left(\frac{{({\mathbf{c}}_p)}_{l,m+1}-{({\mathbf{c}}_p)}_{l,m}}{\Delta r_p},0\right)r_p^2{}_{m+1/2}+\min \left(\frac{{({\mathbf{c}}_p)}_{l,m+1}-{({\mathbf{c}}_p)}_{l,m}}{\Delta r_p},0\right)r_p^2{}_{m-1/2}$$(29)

Moreover, approximated values are required at the cell interfaces $x_{l\pm \frac{1}{2}},$and $r_p{}_{m\pm \frac{1}{2}}$ for Equations (26)(26) $$\frac{d{\mathbf{c}}_{b,l}}{d\tau }=-\frac{{\mathbf{c}}_{b,l+\frac{1}{2}}-{\mathbf{c}}_{b,l-\frac{1}{2}}}{\Delta x}+\frac{\mathbf{P}}{\Delta x}\left[{\left(\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial x}\right)}_{l+\frac{1}{2}}-{\left(\frac{\partial {\mathbf{c}}_{\mathbf{b}}}{\partial x}\right)}_{l-\frac{1}{2}}\right]-\mathbf{\xi }\left({\mathbf{c}}_{b,l}-{\mathbf{c}}_{p,l,N_{r_p}}\right)$$(26) and (28)(28) $$\frac{d{\mathbf{c}}_{p,l,m}}{d\tau }\hspace{0.17em}=\hspace{0.17em}{\mathbf{J}}_{l,m}^{-1}\eta \frac{1}{{ϵ}_p{r}_p^2{}_{m+1/2}\Delta r_p}\left[{({\mathbf{c}}_p)}_{l,m+1/2}-{({\mathbf{c}}_p)}_{l,m-1/2}\right]+F_p\frac{L}{u}{\mathbf{J}}_{l,m}^{-1}\mathbf{R}{r}_{l,m}^{\text{het}}$$(28) . Various methods can be used to approximate them. We are presenting here the first and second-order approximations together with the TVD-RK scheme, used to get a second-order accuracy in time.^[42]

First order scheme

Here, the concentration vector ${\mathbf{c}}_b$ and c_p are approximated at the interfaces of the cell by applying backward difference formula: (30) $${\mathbf{c}}_{b,l+\frac{1}{2}}={\mathbf{c}}_{b,l}\hspace{1em}{\mathbf{c}}_{b,l-\frac{1}{2}}={\mathbf{c}}_{b,l-1}$$(30) (31) $${\mathbf{c}}_{p,l,m+\frac{1}{2}}={\mathbf{c}}_{p,l,m}\hspace{1em}{\mathbf{c}}_{p,l,m-\frac{1}{2}}={\mathbf{c}}_{p,l,m-1}$$(31)

A first-order accurate numerical scheme is obtained in the axial- coordinates by using the aforementioned approximations given by Equations (30)(30) $${\mathbf{c}}_{b,l+\frac{1}{2}}={\mathbf{c}}_{b,l}\hspace{1em}{\mathbf{c}}_{b,l-\frac{1}{2}}={\mathbf{c}}_{b,l-1}$$(30) and (31)(31) $${\mathbf{c}}_{p,l,m+\frac{1}{2}}={\mathbf{c}}_{p,l,m}\hspace{1em}{\mathbf{c}}_{p,l,m-\frac{1}{2}}={\mathbf{c}}_{p,l,m-1}$$(31) .

Second-order scheme

Here, the cell interface concentrations are calculated in the following manner to get a second-order accurate scheme (32) $${\mathbf{c}}_{b,l+\frac{1}{2}}={\mathbf{c}}_{b,l}+\frac{1}{2}\phi \left({\alpha }_{l+\frac{1}{2}}\right)({\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}),\hspace{1em}{\alpha }_{l+\frac{1}{2}}=\frac{{\mathbf{c}}_{b,l+1}-{\mathbf{c}}_{b,l}+\gamma }{{\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}+\gamma }$$(32) (33) $${\mathbf{c}}_{p,l,m+\frac{1}{2}}={\mathbf{c}}_{p,l,m}+\frac{1}{2}\varphi \left({\beta }_{l,m+\frac{1}{2}}\right)({\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}),\hspace{1em}{\beta }_{l,m+\frac{1}{2}}=\frac{{\mathbf{c}}_{p,l,m+1}-{\mathbf{c}}_{p,l,m}+\gamma }{{\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}+\gamma }$$(33)

A flux-limiting high-resolution scheme is produced by Equations (32)(32) $${\mathbf{c}}_{b,l+\frac{1}{2}}={\mathbf{c}}_{b,l}+\frac{1}{2}\phi \left({\alpha }_{l+\frac{1}{2}}\right)({\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}),\hspace{1em}{\alpha }_{l+\frac{1}{2}}=\frac{{\mathbf{c}}_{b,l+1}-{\mathbf{c}}_{b,l}+\gamma }{{\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}+\gamma }$$(32) and (33)(33) $${\mathbf{c}}_{p,l,m+\frac{1}{2}}={\mathbf{c}}_{p,l,m}+\frac{1}{2}\varphi \left({\beta }_{l,m+\frac{1}{2}}\right)({\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}),\hspace{1em}{\beta }_{l,m+\frac{1}{2}}=\frac{{\mathbf{c}}_{p,l,m+1}-{\mathbf{c}}_{p,l,m}+\gamma }{{\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}+\gamma }$$(33) . Here, $\gamma ={10}^{-10}$ is used to avoid division by zero. The flux limiting functions $\phi$ and $\varphi$ are used to preserve the local monotonicity (positivity) of the numerical scheme.^[42] They are defined as (34) $$\phi \left({\alpha }_{l+\frac{1}{2}}\right)=\max \left(0,\min \left(2{\alpha }_{l+\frac{1}{2}},\min \left(\frac{1}{3}+\frac{2}{3}{\alpha }_{l+\frac{1}{2}},2\right)\right)\right)$$(34) (35) $$\varphi \left({\beta }_{l,m+\frac{1}{2}}\right)=\max \left(0,\min \left(2{\beta }_{l,m+\frac{1}{2}},\min \left(\frac{1}{3}+\frac{2}{3}{\beta }_{l,m+\frac{1}{2}},2\right)\right)\right)$$(35)

The high-resolution scheme is given by Equations (32)(32) $${\mathbf{c}}_{b,l+\frac{1}{2}}={\mathbf{c}}_{b,l}+\frac{1}{2}\phi \left({\alpha }_{l+\frac{1}{2}}\right)({\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}),\hspace{1em}{\alpha }_{l+\frac{1}{2}}=\frac{{\mathbf{c}}_{b,l+1}-{\mathbf{c}}_{b,l}+\gamma }{{\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}+\gamma }$$(32) and (33)(33) $${\mathbf{c}}_{p,l,m+\frac{1}{2}}={\mathbf{c}}_{p,l,m}+\frac{1}{2}\varphi \left({\beta }_{l,m+\frac{1}{2}}\right)({\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}),\hspace{1em}{\beta }_{l,m+\frac{1}{2}}=\frac{{\mathbf{c}}_{p,l,m+1}-{\mathbf{c}}_{p,l,m}+\gamma }{{\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}+\gamma }$$(33) cannot be applied at the boundary intervals. Therefore, the first order backward approximations are applied to the boundary intervals. The fluxes at all other interior interval are computed by using Equations (32)(32) $${\mathbf{c}}_{b,l+\frac{1}{2}}={\mathbf{c}}_{b,l}+\frac{1}{2}\phi \left({\alpha }_{l+\frac{1}{2}}\right)({\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}),\hspace{1em}{\alpha }_{l+\frac{1}{2}}=\frac{{\mathbf{c}}_{b,l+1}-{\mathbf{c}}_{b,l}+\gamma }{{\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}+\gamma }$$(32) and (33)(33) $${\mathbf{c}}_{p,l,m+\frac{1}{2}}={\mathbf{c}}_{p,l,m}+\frac{1}{2}\varphi \left({\beta }_{l,m+\frac{1}{2}}\right)({\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}),\hspace{1em}{\beta }_{l,m+\frac{1}{2}}=\frac{{\mathbf{c}}_{p,l,m+1}-{\mathbf{c}}_{p,l,m}+\gamma }{{\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}+\gamma }$$(33) . It is to be noted that this first-order scheme at the boundary cells does not reduce the overall accuracy of this method.^[42]

To reassure the same second-order accuracy in the time coordinate, a second-order accurate TVD-RK scheme is used for solving Equations (32)–(35). Denoting the right-hand-side of Equations (32)(32) $${\mathbf{c}}_{b,l+\frac{1}{2}}={\mathbf{c}}_{b,l}+\frac{1}{2}\phi \left({\alpha }_{l+\frac{1}{2}}\right)({\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}),\hspace{1em}{\alpha }_{l+\frac{1}{2}}=\frac{{\mathbf{c}}_{b,l+1}-{\mathbf{c}}_{b,l}+\gamma }{{\mathbf{c}}_{b,l}-{\mathbf{c}}_{b,l-1}+\gamma }$$(32) and (33)(33) $${\mathbf{c}}_{p,l,m+\frac{1}{2}}={\mathbf{c}}_{p,l,m}+\frac{1}{2}\varphi \left({\beta }_{l,m+\frac{1}{2}}\right)({\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}),\hspace{1em}{\beta }_{l,m+\frac{1}{2}}=\frac{{\mathbf{c}}_{p,l,m+1}-{\mathbf{c}}_{p,l,m}+\gamma }{{\mathbf{c}}_{p,l,m}-{\mathbf{c}}_{p,l,m-1}+\gamma }$$(33) by $\mathcal{L}({\mathbf{c}}_b,{\mathbf{c}}_p{|}_{r_p=1})$ and $\mathcal{M}({\mathbf{c}}_p),$ a second-order TVD RK scheme updates ${\mathbf{c}}_b$ and c_p through the following two stages: (36) $${\mathbf{c}}_b^{(1)}= {\mathbf{c}}_b^n+\Delta \tau \mathcal{L}({\mathbf{c}}_b^n,{\mathbf{c}}_p^n{|}_{r_p=1}),\hspace{1em}{\mathbf{c}}_b^{n+1}=\frac{1}{2}[{\mathbf{c}}_b^n+{\mathbf{c}}_b^{(1)}+\Delta \tau \mathcal{L}({\mathbf{c}}_b^{(1)},{\mathbf{c}}_p^{(1)}{|}_{r_p=1})]$$(36) (37) $$\hspace{1em}{\mathbf{c}}_p^{(1)}={\mathbf{c}}_p^n+\Delta \tau \mathcal{M}({\mathbf{c}}_p^n),\hspace{1em}{\mathbf{c}}_p^{n+1}=\frac{1}{2}[{\mathbf{c}}_p^n+{\mathbf{c}}_p^{(1)}+\Delta \tau \mathcal{M}({\mathbf{c}}_p^{(1)})]$$(37) where ${\mathbf{c}}_b^n$ and ${\mathbf{c}}_p^n$ denote the solutions at the previous time step τⁿ and ${\mathbf{c}}_b^{n+1},$ ${\mathbf{c}}_p^{n+1}$ are updated solutions at the next time step ${\tau }^{n+1}.$ Moreover, $\Delta \tau$ represents the time step which is calculated under the following Courant-Friedrichs-Lewy (CFL) condition: (38) $$\Delta \tau \le 0.5\min \left(\Delta x,\Delta x^2\min (P{e}_{b,i},P{e}_T),\frac{\Delta r}{2\max (J_{l,m}^{-1}\eta )},\frac{\Delta r^2}{\max (J_{l,m}^{-1}\eta )}\right)$$(38)

The aforementioned numerical algorithm was programmed in C language with 100 × 50 grid points. An Intel processor of core-i5 laptop computer Intel(R) Core(TM) $i5-3210MB$ with a random access memory (RAM) size of 4096 MB was used to execute the program.

Performance consistency criteria

Optimization of performance for non-isothermal reactive chromatographic processes requires suitable performance criteria. Here, we propose the following integral consistency tests for mass and energy balance equations to check the accuracy of the numerical scheme and the formulated model equations. In this section, these tests are performed to analyze the conservation of mass and energy balances for three-component nonlinear reactive model considering the reaction of type, $A\hspace{0.17em}\leftrightarrows \hspace{0.17em}B+C.$

Identity of integrated extents of reaction

A general approach for solving a system with multiple reactions is to calculate the extent to which reactions proceeds. It can be applicable whenever we know: (i) the complete composition of either the inlet or the existing stream from a reactor and (ii) the one constraint (conservation, selectivity, second composition, equilibrium constant, and etc.) for each reaction. Every chemical reaction adheres to the law of conservation, as mass cannot be created or destroyed during a chemical reaction. Thus, the total amount of concentration injected at the column inlet remains conserved during a chemical reaction. Conservation constraint is required to see the extent at which the reactants are converted into products during a reversible reaction, as for such reactions, reactants never attain a complete conversion. A change in the mole numbers n_i of reactants and products participating in a chemical reaction depends on the stoichiometry of the reaction. For the considered chemical reaction A $\rightleftarrows$ B + C, the integrated extent of reaction ξ is defined by the following relation: (39) $$\xi =n_A^{\text{inj}}-n_A^{\text{out}}=n_B^{\text{out}}+n_C^{\text{out}}$$(39) here ξ denotes the total changes occurred in the number of moles due to the chemical reaction and $n_i^{\text{inj}}=C_i^{\text{inj}}{V}^{\text{inj}}$ quantifies the number of moles injected into the column. Here, $V^{\text{inj}}$ represents the injected volume at the inlet of the column and $t_{\text{inj}}$ is the time of injection. For the considered test problems, concentrations of the products $c_B^{\text{inj}}$ and $c_C^{\text{inj}}$ at the column inlet are taken to be zero.

The mole numbers of reactant and products at the column outlet can be calculated by using the following integral formula: (40) $$n_i^{\text{out}}=\dot{V}{\int }_0^{t^{*}}{c}_i(t,z=L)\mathrm{d}t,\hspace{1em}i=A,B,C$$(40) here $\dot{V}$ represents the volumetric flow rate. In this work, the trapezoidal rule is used to approximate the above integral numerically. Moreover, simulations are considered for a longer time in each case study in order to bring back the system in its initial equilibrium state (time $t^{*}$).

The three values of ξ_k obtained from Equation (39)(39) $$\xi =n_A^{\text{inj}}-n_A^{\text{out}}=n_B^{\text{out}}+n_C^{\text{out}}$$(39) can be utilized to calculate the standard deviation as follows: (41) $${\sigma }_{\xi ,i}\left[\%\right]=100\times \sqrt{\frac{\sum _{i=1}^i{({\xi }_i-\overline{\xi })}^2}{3}} ,\hspace{1em}i=A,B,C$$(41)

In the above expression, the average of three ξ_i for $i=A,B,C$ is represented by $\overline{\xi }.$ This standard deviation tends to zero, if the mass balances respect reaction stoichiometry.

Integrated energy balance considering extent of reaction

Energy balance is a useful quantity for analyzing the performance of non-isothermal reactors because the temperature of a chemical reactor is determined by the energy balance for the reactor. Here, the computation of energy balance is performed by comparing the enthalpies entering ($\Delta H^{\text{inj}}$) and leaving ($\Delta H^{\text{out}}$) the system. These enthalpies are defined as (42) $$\begin{array}{c}\Delta H^{\text{inj}}={\rho }^{\mathrm{L}}{c}_p^{\mathrm{L}}\dot{V}{\int }_0^{t^{*}}(T^{\text{inj}}-T^{\text{ref}}) \mathrm{d}t \\ \Delta H^{\text{out}}={\rho }^{\mathrm{L}}{c}_p^{\mathrm{L}}\dot{V}{\int }_0^{t^{*}}(T(t,z=L)-T^{\text{ref}}) \mathrm{d}t \end{array}$$(42) where $\Delta H^{\text{inj}}=0,$ for $T^{\text{inj}}=T^{\text{ref}}.$ Moreover, there will be no overall sorption effect in the case of a complete adsorption and desorption cycle (for sufficiently long $t^{*}$). On the basis of reaction’s effect quantified by heat of reaction $\Delta H_{\mathrm{R}}$ and the extent of reaction ξ, we can derive the following balanced equation, e.g. Equation (41)(41) $${\sigma }_{\xi ,i}\left[\%\right]=100\times \sqrt{\frac{\sum _{i=1}^i{({\xi }_i-\overline{\xi })}^2}{3}} ,\hspace{1em}i=A,B,C$$(41) : (43) $$\Delta H^{\text{out}}+(\Delta H_{\mathrm{R}})\overline{\xi }=0$$(43)

To achieve the accurate numerical simulation, the aforementioned Equation (43)(43) $$\Delta H^{\text{out}}+(\Delta H_{\mathrm{R}})\overline{\xi }=0$$(43) is required to be fulfilled as proof. Let the right-hand-side of Equation (43)(43) $$\Delta H^{\text{out}}+(\Delta H_{\mathrm{R}})\overline{\xi }=0$$(43) , which represents the error in numerical simulation is denoted as $\Delta H_{\text{err}}.$ There are several sources of numerical errors which should be taken into account, such as round off errors, discretization errors, errors in the numerical integrations of the outlet profiles, and etc. Due to these errors, the value $\Delta H_{\text{err}}$ might not be exactly zero, i.e. (44) $$\Delta H^{\text{out}}+(\Delta H_{\mathrm{R}})\overline{\xi }=\Delta H_{\text{err}}$$(44)

Thus, to achieve better fulfillment of the joint mass and energy balances, the value of $\Delta H_{\text{err}},$ should be as small as possible. Due to the accumulation of all possible errors in this more critical consistency check one could expect large errors in $\Delta H_{\text{err}}$ as compared to errors in ξ. A relative percentage error in this energy can be defined as (45) $$E_{\mathrm{H}}\left[\%\right]=100\times \left|\frac{\Delta H_{\text{err}}}{\Delta H_{\mathrm{R}}\overline{\xi }}\right|$$(45)

Numerical case studies

In this section, some test problems are considered to figure out the effects of different parameters that influence separation and conversion in non-isothermal reactive liquid chromatography. The considered case studies also explain coupling between elution profiles of concentration and temperature in the non-isothermal chromatographic reactor. The axial dispersion $D_{b,i}$ is assumed to be the same for all components. Moreover, the values of some kinetic parameters, such as effective mass transfer coefficient $k_{\text{eff},\mathrm{i}},$ effective internal pore diffusivity $D_{\text{eff}},$ effective heat conductivity coefficient ${\lambda }_{\text{eff},\mathrm{z}},$ effective heat transfer coefficient $h_{\text{eff}},$ and the enthalpy of adsorption $\Delta H_{A,j}=\Delta H_A,$ are assumed to be the same for all components. However, the current model equations and the proposed numerical solution technique allow different values for these quantities with respect to components. All the parameters used in these test problems are listed in Table 1 which are in accordance with ranges of parameters typically encountered in HPLC applications. The parameters used in the test problems are similar to those presented in.^[26,28]

Table 1. Parameters used in the test problems.

Description	Symbols	Value	Unit
Bed void volume fraction	${ϵ}_b$	0.4	–
Particle porosity	${ϵ}_p$	0.333	–
Axial Peclet number for concentration	$P{e}_{z,i}$	1500	–
Axial Peclet number for temperature	$P{e}_{z,T}$	1500	–
Dimensionless constant	ηi	2.7142	–
Dimensionless constant	ηT	0.6785	–
Dimensionless constant	βi	40	–
Dimensionless constant	βT	40	–
Henry’s constant for component A	$a_1^{\text{ref}}$	1.5	–
Henry’s constant for component B	$a_2^{\text{ref}}$	0.3	–
Henry’s constant for component C	$a_3^{\text{ref}}$	2.5	–
Interstitial velocity	u	$57.5683\times {10}^{-4}$	m/min
Column length	L	0.25	m
Injection time (dimensionless)	${\tau }_{\text{inj}}$	1.0	–
Injected concentration of component A	$c_A^{\text{inj}}$	3.0	mol/L
Injected concentration of component B	$c_B^{\text{inj}}$	0.0	mol/L
Injected concentration of component C	$c_C^{\text{inj}}$	0.0	mol/L
Reaction equilibrium constant	$K_{eq}^{\text{het}}$	2.0	mol/L
Heterogeneous reaction rate constant	$k^{\text{het}}$	$6.0\times {10}^{-3}$	$mi{n}^{-1}$
Reference temperature	$T^{\text{ref}}$	298	K
Activation energy	$E_A^{\text{het}}$	70	kJ/mol
General gas constant	Rg	0.008314	kJ/molK
Density times heat capacity in liquid phase	${\rho }^L{c}_p^L$	4.0	kJ/Kl
Density times heat capacity in solid phase	${\rho }^S{c}_p^S$	8.0	kJ/Kl

Isothermal case ($\Delta H_A=\Delta H_R=0kJ/\mathrm{mol}$)

Figure 1 demonstrates the isothermal behavior of the process for the considered reversible reaction $A\hspace{0.17em}\leftrightarrows \hspace{0.17em}B+C.$ In this case, the reactant A is injected as a pulse using Danckwerts BCs. Here, $\Delta H_A=\Delta H_R=0 kJ/\mathit{mol},$ thus, no variations in the temperature profile could be expected, i.e., it remains constant. The result indicates that reactant A is continuously converting into products B and C during its propagation through the column. At the same time, the separation between the components A, B, and C is also visible due to their different affinities with the solid bed. This case study can be taken as a reference to inspect non-isothermal behavior.

Figure 1. Isothermal case: $\Delta H_A=0 kJ/\mathit{mol}$ and $\Delta H_R=0 kJ/\mathit{mol}.$ Moreover, $b_j^{\text{ref}}=0$ for j = 1, 2, 3.

Effects of enthalpy of reaction ($\Delta H_A=0 kJ/\mathrm{mol},\Delta H_R\ne 0 kJ/\mathrm{mol}$)

The influence of enthalpy of reaction $\Delta H_R$ on concentration and temperature profiles is investigated in this test problem while keeping the enthalpy of adsorption as $\Delta H_A=0 kJ/\mathit{mol}.$ The results are displayed in Figure 2 for two different values of $\Delta H_R=-20$ and $\Delta H_R=-40.$ Now increasing the value of exothermic reaction has a visible effect on concentration and temperature profiles. Significant rise in the temperature profile can be observed for $\Delta H_R=-20$ and $\Delta H_R=-40.$ It can also be observed that an increase in the magnitude of heat of reaction $\Delta H_R$ enhances the conversion of reactant A into products B and C. For $\Delta H_R=-20$ the rate of conversion is $36(\%)$ which improves further to $41(\%)$ for $\Delta H_R=-40.$ The same trend can also be seen in Table 2.

Figure 2. Effect of enthalpy of reaction: (a,b) is for $\Delta H_A=0 kJ/\mathit{mol}$ and $\Delta H_R=-20 kJ/\mathit{mol},$ (c,d) is for $\Delta H_A=0 kJ/\mathit{mol},$ and $\Delta H_R=-40 kJ/\mathit{mol}.$ Moreover, $b_j^{\text{ref}}=0$ for j = 1, 2, 3.

Table 2. Here, $X_A(\%)=100\times (n_A^{\text{inj}}-n_A^{\text{out}})/n_A^{\text{inj}}.$

Parameters (kJ/mol)	ξA (mol)	ξB (mol)	ξC (mol)	${\sigma }_{\xi ,k}$(%)	$X_A(\%)$	$\Delta H^{\mathit{out}}(kJ)$	$\Delta H_{\mathit{err}}(kJ)$	EH (%)
$\Delta H_A=0,\Delta H_R=0$	0.0049	0.0064	0.0073	0.1006	32.6457	0	0	–
$\Delta H_A=0,\Delta H_R=-20,E_A=70$	0.0054	0.0071	0.0082	0.1118	36.2816	0.2744	0.1365	0.9904
$\Delta H_A=0,\Delta H_R=-40,E_A=70$	0.0062	0.0080	0.0092	0.1265	41.0616	0.2749	−0.0372	0.1192
$\Delta H_A=-20,\Delta H_R=0,E_A=70$	0.0063	0.0081	0.0094	0.1282	41.6865	0.2740	0.2740	–
$\Delta H_A=-40,\Delta H_R=0,E_A=70$	0.0061	0.0082	0.0091	0.1269	40.6351	0.2740	0.2740	–
$\Delta H_A=-60,\Delta H_R=0,E_A=70$	0.0053	0.0078	0.0080	0.1218	35.5920	0.2741	0.2741	–
$\Delta H_A=-60,\Delta H_R=-20,E_A=70$	0.0054	0.0077	0.0081	0.1193	35.9597	0.2744	0.1331	0.9417
$\Delta H_A=-60,\Delta H_R=-20,E_A=100$	0.0072	0.0100	0.0108	0.1538	47.8497	0.2745	0.0884	0.4753
$\Delta H_A=-60,\Delta H_R=-20,E_A=120$	0.0085	0.0116	0.0128	0.1801	56.9529	0.2745	0.0546	0.2481

Effects of enthalpy of adsorption ($\Delta H_A\ne 0 kJ/\mathrm{mol},\Delta H_R=0 kJ/\mathrm{mol}$)

In Figures 3 and 4, the effects of enthalpy of adsorption $\Delta H_A$ are analyzed, while the enthalpy of reaction is kept as $\Delta H_R=0.$ The results are obtained by considering three different values of enthalpy of adsorption, i.e. $\Delta H_A=-20 kJ/\mathit{mol},\Delta H_A=-40 kJ/\mathit{mol}$ (c.f Figure 3) and $\Delta H_A=-60 kJ/\mathit{mol}$ in (c.f Figure 4). The separation of products from reactant in the chromatographic reactor is affected by differences in the adsorption ability of components and by temperature through the individual enthalpies of adsorption, $\Delta H_{A,i}.$ Figures 3 and 4 reveal completely different behavior of concentration and temperature profiles in comparison to those obtained in the previously considered cases. Moreover, it can be observed from quantitative data presented in Table 2 that an increase in the magnitude of $\Delta H_{A,i}$ leads to a significant reduction in conversion and separation. The conversion rate is 41(%), when $\Delta H_A$ is $-20 kJ/\mathit{mol}$ and reduces to 35%, when $\Delta H_A$ is $-60 kJ/\mathit{mol}.$

Figure 3. Effects of enthalpy of adsorption: (a,b) is for $\Delta H_A=-20 kJ/\mathit{mol},$ (c,d) are for $\Delta H_A=-40 kJ/\mathit{mol}.$ In all, Figure 3a–d, $\Delta H_R=0 kJ/\mathit{mol}.$ Moreover, $b_j^{\text{ref}}=0$ for j = 1, 2, 3.

Figure 4. Effect of enthalpy of adsorption: $\Delta H_A=-60 kJ/\mathit{mol}$ and $\Delta H_R=0 kJ/\mathit{mol}.$ Moreover, $b_j^{\text{ref}}=0$ for j= 1, 2, 3.

Joint effects of enthalpies of reaction and adsorption ($\Delta H_A\ne 0 kJ/\mathrm{mol},\Delta H_R\ne 0 kJ/\mathrm{mol}$)

In this test problem, combine effects of enthalpies of adsorption and reaction are evaluated to completely investigate the behavior of non-isothermal chromatographic reactors. The enthalpy of adsorption $\Delta H_A=-60 kJ/\mathit{mol}$ and the enthalpy of reaction $\Delta H_R=-20 kJ/\mathit{mol}$ are considered to be the same for all components. The results are shown in Figure 5. Comparison of Figures 4 and 5 affirmed that concentration and temperature profiles follow similar trends. As reactant A is eluted in the middle of products, the reaction yields more amount of pure products. In Figure 5c one can observe a noticeable decrease in the peak height of reactant A and a significant rise in the peak heights of products B and C. Further, it can be observed from the simulated results shown in Table 2 that larger values of activation energy, i.e. $E_A^{\text{het}}=100 kJ/\mathit{mol}$ and $E_A^{\text{het}}=120 kJ/\mathit{mol},$ enhance the reaction rate. Resultantly, more amount of reactant is converted into products. The rate of conversion increases up to 47 (%) and 56 (%) for the parameters considered. The collective errors in the integral mass and energy balances, denoted by E_H (c.f. Equation (45)(45) $$E_{\mathrm{H}}\left[\%\right]=100\times \left|\frac{\Delta H_{\text{err}}}{\Delta H_{\mathrm{R}}\overline{\xi }}\right|$$(45) ), is less than 1(%) in all the test problems. The corresponding quantitative values given in Table 2 demonstrate the precision of numerical solutions.

Figure 5. Effects of both enthalpies of adsorption and reaction: Here, $\Delta H_A\ne 0 kJ/\mathit{mol}$ and $\Delta H_R\ne 0 kJ/\mathit{mol}.$ Further, $b_j^{\text{ref}}=0$ for j = 1, 2, 3.

Effects of the model parameters η, η_T, β, and β_T

In this subsection, the effects of the intraparticle diffusion coefficients for mass and energy (η and η_T) and mass transfer coefficients for mass and energy (β and β_T) are investigated. Figure 6 gives the results showing the effects of η and η_T. It is important to mention that for the cases where the value of η is altered, the value of η_T remains the same as given in Table 1 and vice-versa. The plots show that small values of η (or η_T) give broadened profiles of low peak heights for both the concentration and temperature. Very similar results are witnessed in Figure 7 when the value of β (or β_T) is reduced.

Figure 6. Effects of η and η_T which represents intraparticle diffusion coefficients. Here, $b_j^{\text{ref}}=0$ (j = 1, 2, 3).

Figure 7. Effects of β and β_T which represents the mass transfer coefficients. Here, $b_j^{\text{ref}}=0$ for j = 1, 2, 3.

Effects of the density times heat capacity ratio of solid to liquid phases ($\frac{{\rho }^S{c}^S}{{\rho }^L{c}^L}$)

In this problem, the effects of varying density time heat capacity ratio $\frac{{\rho }^S{c}^S}{{\rho }^L{c}^L}$ are investigated under the influence of both enthalpies of adsorption and reaction. Retention times and the speed of concentration and temperature profiles inside the reactor are described by this ratio. The results are shown in Figure 8. For $\frac{{\rho }^S{c}^S}{{\rho }^L{c}^L}=0.2$ obtained by taking ${\rho }^s{c}_P^s=8KJ/1K$ and ${\rho }^L{c}_P^L=40KJ/1K,$ see Figure 8a,b, the adsorption related to positive peak of thermal wave is moving slightly faster than concentration pulses and the desorption related negative peak is coupled with concentration pulses. For $\frac{{\rho }^S{c}^S}{{\rho }^L{c}^L}=1,$ obtained by taking ${\rho }^s{c}_P^s=8KJ/1K$ and ${\rho }^L{c}_P^L=8KJ/1K$ (Figure 8c,d), the mean retention times for concentration and temperature profiles are the same and, thus, travel with the same velocity. The figure depicts the coupling between the predicted profiles. For $\frac{{\rho }^S{c}^S}{{\rho }^L{c}^L}=5.0$ obtained by taking ${\rho }^s{c}_P^s=40KJ/1K$ and ${\rho }^L{c}_P^L=8KJ/1K,$ see (Figure 8e,f), the adsorption related peak of temperature is coupled with concentration profiles while the desorption related peak elutes later from the reactor. It can also be observed that more reactant is converted into product for the case $\frac{{\rho }^S{c}^S}{{\rho }^L{c}^L}=1.$

Figure 8. Effects of the ratio $\frac{{\rho }^S{c}_p^S}{{\rho }^L{c}_p^L}$ on concentration and temperature profiles.

Isotherm nonlinearities with respect to concentration

Figures 9 and 10 show the results of numerical calculations for fully nonlinear isotherm given by Equation (7)(7) $$q_{p,i}=\frac{a_i^{\text{ref}}{c}_{p,i}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-\Delta H_{A,i}}{R_g}\left(\frac{1}{T_p}-\frac{1}{T^{\text{ref}}}\right)\right)}{1+\sum _{j=1}^{N_c}{b}_j^{\text{ref}}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-\Delta H_{A,j}}{R_g}\left(\frac{1}{T_p}-\frac{1}{T^{\text{ref}}}\right)\right)c_{p,j}}$$(7) with $b_j^{\text{ref}}=1$ for j = 1, 2, and 3 and $b_1^{\text{ref}}=2,b_2^{\text{ref}}=1,b_3^{\text{ref}}=3,$ respectively. The remaining parameters are exactly the same as given in Table 1. Figures 9a,b and 10a,b give the results for isothermal condition, while Figures 9c,d and 10c,d display the results for non-isothermal condition. A typical Langmuir behavior can be recognized from the right peak tailings of the concentration profiles. Due to non-linearity in the isotherm and in the reaction term, less amount of reactant is converted into the products. The conversion of reactant into the product is improved in the non-isothermal case as compared to the isothermal case. It can also be observed that conversion of reactant decreases further as the values of nonlinearity coefficient $b_j^{\text{ref}}$ increase. Moreover, in the non-isothermal case, significant deviation in the temperature profile could be seen from the reference temperature of 298 K. It should be finally mentioned that the numerical method applied is capable to easily handle other types of non-linearities.

Figure 9. Fully nonlinear isotherm (c.f. Equations (7)(7) $$q_{p,i}=\frac{a_i^{\text{ref}}{c}_{p,i}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-\Delta H_{A,i}}{R_g}\left(\frac{1}{T_p}-\frac{1}{T^{\text{ref}}}\right)\right)}{1+\sum _{j=1}^{N_c}{b}_j^{\text{ref}}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-\Delta H_{A,j}}{R_g}\left(\frac{1}{T_p}-\frac{1}{T^{\text{ref}}}\right)\right)c_{p,j}}$$(7) and (8)(8) $$r^{\text{het}}=k^{\text{het}}(T_p)\left(q_{p,A}-\frac{q_{p,B}{q}_{p,C}}{K_{eq}^{\text{het}}}\right)$$(8) ): Results of numerical calculations for isothermal and non-isothermal conditions. Here, $b_1^{\text{ref}}=1,b_2^{\text{ref}}=1,$ and $b_3^{\text{ref}}=1.$

Figure 10. Fully nonlinear isotherm (c.f. Equation (7)(7) $$q_{p,i}=\frac{a_i^{\text{ref}}{c}_{p,i}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-\Delta H_{A,i}}{R_g}\left(\frac{1}{T_p}-\frac{1}{T^{\text{ref}}}\right)\right)}{1+\sum _{j=1}^{N_c}{b}_j^{\text{ref}}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-\Delta H_{A,j}}{R_g}\left(\frac{1}{T_p}-\frac{1}{T^{\text{ref}}}\right)\right)c_{p,j}}$$(7) and (8)(8) $$r^{\text{het}}=k^{\text{het}}(T_p)\left(q_{p,A}-\frac{q_{p,B}{q}_{p,C}}{K_{eq}^{\text{het}}}\right)$$(8) ): Results of numerical calculations for isothermal and non-isothermal conditions. Here, $b_1^{\text{ref}}=2,b_2^{\text{ref}}=1$ and $b_3^{\text{ref}}=3.$

Conclusion

A multi-component non-isothermal GRM of nonlinear reactive chromatography was formulated and numerically approximated. The mass and energy balances describing non-isothermal reactive processes form a system of convection-diffusion reaction PDEs coupled with an algebraic equation for adsorption isotherms and reaction. An HR-FVS was implemented to numerically approximate the model equations. Several case studies of three-component reactions were considered and analyzed to demonstrate the influence of different parameters on the process performance, such as rate of reaction, retention time, injected volume of reactant, adsorption enthalpies, density time heat capacity ratio of solid to liquid phases, and etc. Consistency tests related to mass and energy conservation were carried out to verify accuracy of the suggested numerical algorithm. It was observed that non-isothermal reactors are more useful to achieve a high conversion rate as compared to isothermal reactors. The computed results are very useful tools for understanding the transport mechanisms in non-isothermal reactors to scale up physiochemical parameters that affect the reactor performance and to optimize experimental conditions.