Numerical approximation of a nonisothermal two-dimensional general rate model of reactive liquid chromatography

Correspondence Shamsul Qamar, Department of Mathematics, COMSATS University Islamabad, Islamabad 45550, Pakistan. Email: shamsul.qamar@comsats.edu.pk Abstract A nonlinear and nonisothermal two-dimensional general rate model is formulated and approximated numerically to allow quantitatively analyzing the effects of temperature variations on the separations and reactions in liquid chromatographic reactors of cylindrical geometry. The model equations form a nonlinear system of convectiondiffusion-reaction partial differential equations coupled with algebraic equations for isotherms and reactions. A semidiscrete high-resolution finite volume method is modified to approximate the system of partial differential equations. The coupling between the thermal waves and concentration fronts is demonstrated through numerical simulations, and important parameters are pointed out that influence the reactor performance. To evaluate the precision of the model predictions, consistency checks are successfully carried out proving the accuracy of the predictions. The results allow to quantify the influence of thermal effects on the performance of the fixed beds for different typical values of enthalpies of adsorption and reaction and axial and radial Peclet numbers for mass and heat transfer. Furthermore, they provide useful insight into the sensitivity of nonisothermal chromatographic reactor operation.


INTRODUCTION
Reactive chromatographic column is a multifunctional reactor, which integrates reaction and separation processes into the same unit. 1 This process has been found capable for separating complex mixtures with desired results. [2][3][4][5][6][7][8][9][10] To achieve reactions and separations simultaneously, the packing materials should behave as a catalyst and should have different affinities toward the reactants and products. For example, consider a reversible reaction process of the type A ⇄ B + C. The reactant A is injected as a rectangular pulse into the chromatographic column. With the presence of the catalyst inside the column, the reactant A reacts and produces products B and C. The situation will be optimum if the reactant A travels in the middle of the products, the forward reaction will be enhanced, and the backward reaction will be suppressed. For a sufficiently long elution time of reactant A inside the column, a complete conversion of component A and production of products B and C can be achieved. 11 Generally, the goal of a reactive separation is to produce a single high-value product. More ideas on the foundations and operations of chromatographic reactors are given by some researchers. 3,5,6,12,13 Studies of thermal reactions in liquid-phase chromatographic reactors have been carried out by only a few groups of researchers. [14][15][16][17][18][19] Especially, authors in Ref. [17], have presented experimental study of exothermic esterification reaction catalyzed by an acidic ion-exchange resin by considering the same A⇄B + C reaction. Modeling of chromatographic operations is very helpful for understanding and analyzing dynamic composition fronts in chromatographic reactors without engaging in extensive experiments. Different one-and two-dimensional (2D) mathematical models, considering diverse levels of complexities, are available in the literature, which provide detailed information about the mass transfer and partitioning processes. 4,[20][21][22][23][24][25] The nonisothermal models of liquid chromatography can be linearized under the assumptions of small changes in the concentrations and temperature. 25,26 Such linear assumptions are valid when the injected volume of the sample is small and is diluted. The Laplace transform method is often applied to solve the resulting linear model equations analytically due to its moment generating property. 25,26 The derived analytical solutions provide fruitful information about the elution profiles propagating through the column and about the influence of kinetic and thermodynamic parameters on the process. Moreover, these solutions could be helpful to validate the numerical solutions of nonlinear models. However, it has been shown in our previous article 26 that such solutions are only valid and meaningful for small values of enthalpy of adsorption and give overpredicted solutions for larger values of enthalpy of adsorption. On the other hand, for concentrated (or large volume) samples and large temperature variations, consideration of nonlinear models, using nonlinear isotherms and reactions, is required. These models are not solvable analytically, and, thus, a numerical solution technique is required to approximate the model equations. Such models and the corresponding numerical solution techniques are more flexible and general because they handle the cases of both diluted and large volumed samples.
In this work, the general rate model of liquid chromatography is utilized to simulate nonlinear 2D nonisothermal reactive chromatography. The goal is to quantify the effects of temperature variations on conversion and separation in reactive nonisothermal liquid chromatography in the presence of axial and radial dispersions. A second-order semidiscrete finite volume method is used to solve the equations of the model. Case studies of chemical reactions involving three-component mixtures are given to show the coupling of thermal waves and concentration fronts. Different significant parameters that largely affect the performance of the reactor are identified. The results show that a heat capacity ratio of solid phase to liquid phase plays an important role in the case of nonisothermal chromatographic reactor operation. Moreover, small values of radial dispersion coefficients generate radial concentration and temperature gradients inside the reactor.
The novelty of this current work includes (a) the extension of our previous linear models (cf Refs. [25] and [26]) to nonlinear and nonisothermal two dimensional-general rate model (2D-GRM) for liquid chromatography, which incorporates two energy balance equations in the model; (b) introduction of nonisothermal reaction process in the 2D-GRM of liquid chromatography; and (c) application of the 2D finite volume scheme to numerically approximate the current nonisothermal 2D-GRM. Furthermore, as mentioned before, linearizations of adsorption isotherms and reaction term are needed to apply an analytical solution approach. Thus, the same nonlinear and nonisothermal reaction behavior presented in this article cannot be accommodated in the analytical framework presented in our previous article. 26 The considered model equations and the corresponding numerical solution technique are more flexible and general, which are capable of handling the cases of both diluted and large volume samples.
The current nonisothermal 2D model could be very useful in various situations, for example, (i) for an imperfect injection at the inlet of the column (thus introducing a radial profile at the inlet of the column), (ii) when radial temperature gradients are present (these temperature gradients are connected with the radial concentration gradients), and (iii) when the column is not packed homogeneously (this situation is more prevalent for larger columns processes). All of these issues are mostly neglected in practical applications; hence the onedimensional (1D) models become acceptable. However, for their relevance and effects, 2D models are required. With our current nonisothermal 2D model, we could study situations (i) and (ii) only. The situation (iii) needs further extensions in the model equations (eg, considering nonconstant column porosities), which is currently under investigation.
This article is further structured in the following manner: In Section 2, a 2D-GRM is formulated to simulate dynamical processes in nonisothermal liquid chromatographic reactors of cylindrical geometry. Section 3 presents the derivation of a high-resolution flux-limiting finite volume method used for solving the equations. Section 4 presents basic formulations for checking consistency of the results. Some case studies are presented in Section 5 followed by the conclusions in Section 6.

T WO-DIMENSIONAL NONISOTHERMAL GENERAL RATE MODEL
The migration of multicomponent concentration bands in the nonisothermal liquid chromatography is a particular case of the reaction-diffusion processes. Such phenomena can be modeled mathematically by applying partial differential equations, obtained by writing the mass and energy balances in a F I G U R E 1 Schematic diagram of a chromatographic column of cylindrical geometry column slice. In the derivation of our model, it was assumed that the column is thermally insulated and homogeneously packed, a volumetric flow rate is constant, and a heterogeneous reaction is taking place in the solid phase. Moreover, mass and energy transfer between the stationary and mobile phases, axial and radial dispersions, and intraparticle pore diffusion are incorporated in the following mass and energy balance equations.
Let , , and represent the time coordinate, axial coordinate that stretches through the column length, and the radial coordinate along the column radius, respectively. Both the reactant and products move through the axis of the column in the -direction by means of convection and axial dispersion, spreads through the radius of the column in the -direction by radial dispersion, and the reactant decays and produces the products due to chemical reactions in the solid phase. To have significant effects of mass transfer in the radial direction, the following injection conditions are assumed. Thus, we introduce a new parameter̃, which divides the column's inlet cross section into an inner cylindrical core and an outer annular ring (see Figure 1). By doing so, the injection of sample mixtures into the column can be done either via the inner core, via the outer ring, or via the whole cross section (which is the case wheñis set equal to the radius of the column ). Thus, the mass balance equations are given for = 1, 2, … , as where , is the concentration of th component in the bulk phase, , denotes the same component concentration in the pores of the particle, and denotes the fluid phase velocity, which, due to the excess of solvent, is assumed to be constant even in the case of component-specific different axial dispersion coefficients , . Furthermore, the phase ratio is denoted by = (1 − )∕ , where denotes the external porosity, , denotes the dispersion coefficient of the th component in the radial direction, and ext, denotes the coefficient of external mass transfer. Lastly, is the radial coordinate of spherical particles of radius and stands for the number components in the mixture.
The corresponding equation of mass for the solute in the particles pores can be given as where , is the solid-phase concentration at local equilibrium for th component, , is the pore diffusivity for the th component, is the corresponding stoichiometric coefficient of th component, and het denotes the heterogeneous reaction rate in the solid phase.
Owing to the nonisothermal nature of the column, the corresponding energy balance of the column, assuming also heat conductivity in the radial direction of the column, is given as where and are, respectively, temperatures of the bulk fluid and fluid inside the particles pores, ef f, represents the effective axial heat conductivity, ef f , denotes the effective radial heat conductivity, and ℎ ef f is the effective particle to fluid heat transfer coefficient. An energy balance including the possible development of radial temperature profiles inside particles pores is expressed as ( Here, denotes the internal heat diffusivity coefficient, and are the densities of liquid and solid phases, and and are the corresponding heat capacities. The , , , and are considered independent of temperature, which is valid in a small range of temperature. Furthermore, Δ , is the enthalpy of adsorption of th component and Δ denotes the enthalpy of reaction. Assuming the case of a nondispersive, nonreactive, and fully equilibrium chromatography, the propagation speeds of concentration profiles and of the thermal wave can be estimated using Equations (1)-(4) as (cf Ref. [14]) It is evident from Equation (5) that the speed of concentrations is dependent on the respective local gradients of adsorption isotherms and that the propagation speed of temperature wave is influenced by the ratio of density times heat capacity. For less retained components or for smaller gradients of isotherms and for sufficiently larger ratio of S S p L L p , the speeds of concentration fronts are higher than the speed of thermal wave. The nonlinear adsorption isotherm is given as where ref and ref are, respectively, denoting the Henry's constant and nonlinearity coefficient of the th component at reference temperature, is the gas constant, and ref represents the reference temperature.
The chemical reactions inside a chromatographic reactor can be catalyzed homogeneously or heterogeneously or in both ways. In other words, the reaction could occur either in the liquid phase (ie, homogenous reaction) or in the particle phase (ie, heterogeneous reaction) or in both phases. In the case of homogeneously catalyzed reaction, the separation of catalyst has to be taken into account. On the other hand, heterogeneously catalyzed reactions usually occur in the case of 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 (solidphase) reaction. The corresponding reaction rate for a threecomponent model reaction ( ⇄ + ) is given as Here, het and het , respectively, represent the forward heterogeneous rate of reaction constant and reaction equilibrium constant. The Arrhenius equation is used to characterize the temperature effects on the chemical reaction rates using the activation energies het : Next, we introduce the following dimensionless variables for a reduction of the equations parameters: where denotes the column length, , and , are the axial Peclet numbers, , and , are the radial Peclet numbers, and and are Biot numbers for mass and energy, respectively. Further, , , , and are the dimensionless constants. After using the above dimensionless parameters in the mass and energy balances (cf Equations 1-4), we obtain for = 1, 2, … , : , The initial conditions are given as , ( , , , 0) = 0, ( , , , 0) where init and init represent the initial bulk and particle temperatures. The following boundary conditions are considered for Equations (10) and (12) at = 0 and = 1: Also, for Equations (10) and (12), the following Danckwerts boundary conditions are used for injections via the inner region: Here, the symbols  . (18e)

NUMERICAL SCHEME
Different numerical schemes can be utilized to approximate the solutions of the above model equations (see Refs. [21,22,27,28]) and the references therein. In this study, a semidiscrete finite volume method is suggested for approximating the model equations in the axial and radial coordinates. The method is simple, compact, easily implemented, and second-order accurate. 27,29 A second-order total variation diminishing Runge-Kutta (TVD-RK) scheme is utilized to solve the system of ordinary differential equation in the time coordinate. 27,30 The order of accuracy of our suggested scheme has already been verified analytically and numerically in our previous article. 27 It was proved there that this scheme is second-to-third order accurate. This selected second-order Runge-Kutta method for the time discretization and the considered flux-limiting function for calculating fluxes at the cell interfaces guarantee the second-order accuracy of the scheme. 27,29 Furthermore, the scheme is capable of capturing (resolve) accurately sharp fronts of the concentration and temperature profiles. We get the following system of equations from Equations (10)-(13) by using the adsorption isotherm in Equation (6) and considering a mixture of three components: If all the nonlinearity coefficients ref are zero then the above Jacobian matrix becomes very simple. To derive the scheme, we first discretize the domain of computations.
Therefore, for , , , = , , , ( ) On integrating Equation (19) over the interval and utilizing Equations (25) and (26), we arrive to the expression of the form: where = 1, 2, … , and = 1, 2, … , . The derivatives appearing in the above equations are approximated as ( ) Integration of Equation (20) where the interface flux is expressed as  (27) and (29) are approximated by using the following schemes along with the TVD-RK scheme to get a second-order accuracy in time. 27,30

First-order scheme
The axial flux vectors and at the interfaces of the cell are approximated as follows: The above approximations above provide a first-order accuracy of the scheme along the axial and radial coordinates.

Second-order scheme
Here, the axial flux vectors are approximated by following expressions: A flux-limiting high-resolution scheme is produced by Equations (33) and (34). A small value of , for example, = 10 −10 , is selected to avoid the situation of division by zero. The local monotonicity of the scheme is preserved by utilizing the flux limiting functions and as defined below (cf Ref. [27]): It is impossible to apply the high-resolution scheme up to the boundary intervals. To overcome this difficulty, the aforementioned first-order approximation of fluxes is used at the interfaces of boundary cells, whereas the suggested secondorder approximation of fluxes is well applicable at the interfaces of interior cells. In our case, the interstitial velocity is positive. Thus, we need to take first-order approximation in the left boundary cell only, whereas such an approximation is not needed in the right boundary cell. Further, it has been found that global accuracy of the suggested algorithm is not diminished by such first-order approximations in the boundary cells. 27 Lastly, to solve Equations (27)-(36), we apply a secondorder TVD-RK scheme to get the second-order accuracy in time. 30 The right-hand side of Equations (27) and (29) is represented by  ( , | =1 ) and ( ). To update the stages of and , the following second-order TVD-RK scheme is utilized 30 : (1) = + Δ ( ), +1 = 1 2 [  (27) and (29), defined below, is used to choose the time-step Δ : ) . (39) The above method was programmed in C language for grid points of 60 × 30 × 10. obtained from our numerical scheme. 14,15 Application of these integral consistency checks for mass and energy balance equations are very useful for the validation of numerical results and for checking the accuracy of the formulated model. The consistency check is required to verify the correctness of the applied numerical scheme and to evaluate the conservativity of mass and energy balances.

Considering a chemical reaction
Let describes a change in the number of moles due to the chemical reaction, that is .
The total moles injected in the column is described as where inj is the volume injected and inj is the injected concentration. Furthermore, at the column outlet, we have wherėis the volumetric flow rate. Standard deviations for the three values of are calculated as , (%) = 100 × √ ∑ =1 ( −̄) 2 3 , = 1, 2, 3.
Here,̄represents the average of for = 1, 2, 3. For the energy balance, we have where Δ inj is the enthalpy entering and Δ out is the enthalpy leaving the system. Lastly, a relative percentage error H (%) is defined as where In the above equation, Δ err = 0 if the equation is satisfied exactly. However, several sources of numerical errors are involved in the computations, such as discretization errors, round off errors, and errors in the numerical integrations of the outlet profiles. Because of these errors, the right-hand side of Equation (46) is not zero, and the smaller is the value of Δ err the better is the fulfillment of the coupled mass and energy balances.

CASE STUDIES ON THERMAL EFFECTS
In this section, the effects of different parameters that affect separation and conversion in nonisothermal reactive liquid chromatography are analyzed. The considered test problems also explain coupling between thermal fronts and concentration in the nonisothermal chromatographic reactor. The value of enthalpy of adsorption is taken the same for all components, that is Δ , = Δ . Moreover, the values of the kinetic parameters , and , are also assumed the same for all the components, although they may vary according to components in practice. All the parameters used in the test problems are listed in Table 1, which were used by Sainio et al. 17 and Vu and Seidel-Morgenstern 19 in their experimental studies for exothermic esterification reaction catalyzed by an acidic ionexchange resin. In all the test problems, except the last test problem for the fully nonlinear case, the reference nonlinearity coefficients ref ( = 1, 2, 3) in Equation (6) are set equal to zero. Further, a grid of mesh points 60 × 30 × 10 is used in all calculations.

Isothermal case ( = = kJ/mol)
Here, we start with an isothermal case which is an ideal reference case for the study of nonisothermal behavior. The reactant was injected in the annular region as a Danckwerts boundary condition (BC). The 1D and three-dimensional (3D) plots of Figure 2 show the behavior of concentration and temperature profiles under isothermal conditions for the considered three components reversible reaction problem in which only the reactant ,1 is injected to the reactor. The results indicate the presence of products ,2 and , 3 . The values of the Henry's constants show strong but not complete separation and conversion. With the current values of Δ = Δ = 0 kJ/mol, no changes are seen in the temperature profile. For the selected value of radial Peclet number = 37.5, which corresponds to the small radial dispersion , radial variations in the concentration profiles can be easily seen in Figure 2C. Because of injection through the inner region, the value of concentration is larger at the column center, that is at = 0.

Effects of enthalpy of reaction (
= kJ/mol, ≠ kJ/mol) Now, we evaluate the effect of the enthalpy of reaction (Δ ) on concentration and temperature profiles, while keeping Δ = 0 kJ/mol and other parameters are the same as given in Table 1. Once again, injection is done through the inner region of the inlet cross section. The results of  Figures 3C and 3D). On the temperature profile, there is a smooth temperature variation and a single peak, indicating that heat release has occurred as compared with the isothermal case of Figure 2. An increase in the concentration profiles of both products can be observed for Δ = −20 kJ/mol and even more clearly for Δ = −40 kJ/mol. This caused an increase in the conversion 1 (%) of ,1 (reactant) to products ,2 and ,3 , which is evident from the listed results in Table 2. The 3D plots of the results in Figure 3 are given in Figure 4. In this case, radial variations can be observed in both concentration and temperature profiles for the selected = 37.5.

Effects of enthalpy of adsorption ( ≠ kJ/mol, = kJ/mol)
This case evaluates the effects of Δ , whereas the enthalpy of reaction is neglected (ie, Δ = 0 kJ/mol). The separation of the products ,2 and ,3 is based on the differences in the adsorption strength, which depends on the temperature through the enthalpy of adsorption. The effects of the enthalpy of adsorption on the concentration profiles are clearly seen in Figure 5. This also reflects the changes seen in the conversion rates given in Table 2, which shows the conversion of products as 33% when Δ = −20 kJ/mol, 32% when Δ = −40 kJ/mol and back to 28% when Δ = −60 kJ/mol. Again, the 3D plots of the results in Figure 5 are presented in Figure 6 to see the occurrence of the radial changes.

Joint effects of enthalpies of adsorption and reaction ( ≠ kJ/mol, ≠ kJ/mol)
For an entire illustration of the nonisothermal behavior of chromatographic reactors, we investigate the joint effects of enthalpies of adsorption and reaction on both the temperature and concentration profiles. We see in Figure 7 A, which is obtained for Δ = −20 kJ/mol and Δ = −60 kJ/mol, that there is a noticeable increase in the peak height of reactant A and a slight decrease in the peak heights of the products ,2 and ,3 as compared to the results in rate is increased from 28% to 41% which is evident from

Effects of the model parameters , , , and
The effects of the intraparticle diffusion coefficients for mass and energy ( and ) and the mass transfer coefficients for mass and energy ( and ), are investigated under this subsection. Figure 8 gives the results showing the effects of and . It is important to mention here that for the cases where the value of is altered, the value of remains as given in Table 1 and vice versa. The plots show that reducing the value of (or ) results in a reduction of the peak heights and a more diffusive profile for both the concentration and temperature. A very similar result can be witnessed in Figure 9 for the cases when the value of (or ) is reduced.

Effects of the ratio of solid to liquid phases density multiplied by heat capacity ( )
Here, we consider the effects of the ratio along with both enthalpies of adsorption and reaction on the concentration and temperature profiles (cf Equation 5), while other parameters remain the same as given in Table 1. The results shown so far have been obtained for the case where both the waves of the temperature and concentration profiles were moving at closer speeds, that is, = 1 implies that both concentration and temperature profiles are almost coupled. This is clearly illustrated in Figure 10 A. For the case = 0.2 (Figure 10 B), the temperature wave is seen to be moving slightly faster than the concentration profile. For = 5 (Figure 10 C), it is evident that the concentration profile moves quicker than that of the temperature profile. Therefore, the concentration and temperature profiles move at different velocities for the cases = 0.2 and = 5.

Effects of the radial and axial Peclet numbers (
The effects of kinetic parameters , , , , , and , , which, respectively, denote radial dispersion coefficients for the mass and energy, and axial dispersion coefficients for the mass and energy, are considered for nonzero enthalpies of reaction and adsorption. It should be noted that for each of the case considered, except the radial and axial Peclet numbers, all other parameters remained the same as given in Table 1. Figure 11 displays the effects of = 1.5 and = 150, which corresponds to the larger and smaller radial dispersion coefficient for mass . We can conclude from the results that slow radial dispersion (large radial Peclet number) causes a more diffusive profile, which can be clearly seen on the plots of Figures 11 E and 11 F. Figure 12 shows a similar result for the effects of , = 1.5 and , = 150, which corresponds to larger and smaller radial dispersion coefficient for energy ef f,r . It should also be noted that for such small values of and , , the results of current 2D-GRM reduces to the 1D-GRM. Figure 13 displays the results showing the effects of varying , ( , ) and , ( ef f,z ) for the components. It is evident that small values of , and , (large values of , and ef f,z ) generate more diffusive and tailed profiles of the components.  Figures 14 and 15 show the 1D and 3D plots of the results of numerical calculations for the nonlinear isotherm given by Equation (6) with = 1 for = 1, 2, 3. All other parameters are the same as given in Table 1. Figures 14 A and 14 B (and Figures 15 A and 15 B) give the results for isothermal condition, whereas Figures 14 C and 14 D (and Figures 15 C and 15 D) display the results for nonisothermal condition. A typical Langmuir behavior can be recognized from the peak tailings of the concentration profiles. The temperature deviates in the nonisothermal case significantly from the reference temperature of 300 K [−6 K,+11 K]. This leads to small but visible deviations in the elution profiles compared to the isothermal case. It should be finally mentioned that the suggested numer-ical method is capable of handling other types of nonlinear isotherms than that quantified by Equation (6).

CONCLUSION
The effect of temperature variation arising from the exothermic reaction on 2D reactive chromatography was investigated by formulating and applying a 2D general rate model of cylindrical geometry. The model forms a system of convection-diffusion reaction partial differential equations. A high-resolution finite volume method was applied to numerically approximate the solutions of the model equations. Several case studies of reversible reactions were carried out. The consistency tests were carried out to verify accuracy of