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


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][3][4][5][6][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 À B, A À B þ C and A þ B À B þ C are extensively investigated in the literature for isothermal case. [14][15][16][17][18][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 À 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][22][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][25][26][27][28][29] However, thermal effects have been discussed extensively for gas chromatography. [30][31][32][33][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][36][37][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 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 eff, i is the effective mass transfer coefficient, R p is the radius of spherical particles, the factor 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 j 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 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 eff ¼ p D p is the effective internal pore diffusivity, r 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 In the above equations, T b is temperature of the bulk fluid, T p is fluid temperature inside the particle pores, k eff, z represents effective axial heat conductivity coefficient, q L is density per unit volume in the liquid phase, c L P is heat capacity for liquid phase and h 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 here DH A, j represents the enthalpy of adsorption of j-th component, DH R is the enthalpy of reaction, r het is the heterogeneous reaction rate, k eff, e denotes the internal heat diffusivity coefficient. Moreover, where q S is the density per unit volume of solid phase and c S p is the corresponding heat capacity. The considered density and heat capacity q L , q S , c L p , and c S p 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 k c and of the thermal wave u T inside the column can be approximated from Equations (1)-(4) as [29] u k c ffi It is clear from Equation (6) that the speed u k c 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 @q k @c k and for a larger ratio of , the speeds u k c 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 where a ref i denotes the Henry's constant and b ref i represents the nonlinearity coefficient of the ith component at the reference temperature, R g is the gas constant, and T 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 À B þ C) is given as here k het represent the heterogeneous forward reaction rate constant and K het eq represent reaction equilibrium constant, respectively. The Arrhenius equation uses activation energy E het A to describe the temperature effect on the chemical reaction rate k het ðT s Þ as an exponential function of absolute temperature: To simplify the notations and reduce the number of parameters appearing in the model equations, the following dimensionless quantities are introduced: here L denotes the column length, while Pe b, i and Pe T is the peclet numbers, B c, i and B T is the Biot numbers for mass and energy, while g c, i , g T , n c, i , and n T are the dimensionless constants. After using the above dimensionless parameters in the mass and energy balances (c.f. Equations (1)-(4)), we obtain (14)) are solved using the appropriate initial and boundary conditions (BCs).
The initial conditions are given as and T init p 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 is temperature of the injected sample. In this work, we have taken T inj b , T init b , and T ref the same. At the column outlet (x ¼ 1), the zero Neumann BCs are utilized: For Equations (12) and (14), radial BCs at r ¼ 0 and r ¼ 1 are expressed as 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][41][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 where The model Equations (18) and (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 Â ½0, 1 which is covered by cells X lm ½x l À 1 2 , x l þ 1 2 Â ½r pm À 1 2 , r pm þ 1 2 for 1 l N x , and 1 m N r p : The characteristic coordinate points in the cell X lm are represented by ðx l , r pm Þ: Here, and for the uniform mesh , Since c b :¼ c b ðs, xÞ and c p :¼ c p ðs, x, r p Þ Therefore, for I l :¼ ½x l À 1 2 , x l þ 1 2 and X lm :¼ ½x l À 1 2 , x l þ 1 2 Â ½r pm À 1 2 , r pm þ 1 2 , the averaged values of the cell c b, l ðsÞ and c p, l, m ðsÞ, expressed at any time s, are given as By integrating Equation (18) over the interval I l and utilizing Equations (24) and (25), we obtain where l ¼ 1, 2, :::, N x , and the differential term of the axial diffusion part is approximated as Integration of Equation (19) over the interval X lm gives where Moreover, approximated values are required at the cell interfaces x l6 1 2 , and r pm6 1 2 for Equations (26) and (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 secondorder accuracy in time. [42] First order scheme Here, the concentration vector c b and c p are approximated at the interfaces of the cell by applying backward difference formula: c p, l, mþ 1 2 ¼ c p, l, m c p, l, mÀ 1 2 ¼ c p, l, mÀ1 A first-order accurate numerical scheme is obtained in the axial-coordinates by using the aforementioned approximations given by Equations (30) and (31).

Second-order scheme
Here, the cell interface concentrations are calculated in the following manner to get a second-order accurate scheme A flux-limiting high-resolution scheme is produced by Equations (32) and (33). Here, c ¼ 10 À10 is used to avoid division by zero. The flux limiting functions u and / are used to preserve the local monotonicity (positivity) of the numerical scheme. [42] They are defined as The high-resolution scheme is given by Equations (32) and (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) and (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 righthand-side of Equations (32) and (33) by Lðc b , c p r p ¼1 Þ and Mðc p Þ, a second-order TVD RK scheme updates c b and c p through the following two stages: 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 ! 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 À B þ C, the integrated extent of reaction n is defined by the following relation: here n denotes the total changes occurred in the number of moles due to the chemical reaction and n inj i ¼ C inj i V inj quantifies the number of moles injected into the column. Here, V inj represents the injected volume at the inlet of the column and t inj is the time of injection. For the considered test problems, concentrations of the products c inj B and c inj C 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: here _ 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 n k obtained from Equation (39) can be utilized to calculate the standard deviation as follows: In the above expression, the average of three n i for i ¼ A, B, C is represented by n: 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 (DH inj ) and leaving (DH out ) the system. These enthalpies are defined as 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 DH R and the extent of reaction n, we can derive the following balanced equation, e.g. Equation (41): To achieve the accurate numerical simulation, the aforementioned Equation (43) is required to be fulfilled as proof. Let the right-hand-side of Equation (43), which represents   DH out þ ðDH R Þ n ¼ DH err Thus, to achieve better fulfillment of the joint mass and energy balances, the value of DH 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 DH err as compared to errors in n. A relative percentage error in this energy can be defined as Numerical case studies In this section, some test problems are considered to figure out the effects of different parameters that influence     parameters used in the test problems are similar to those presented in. [26,28] Isothermal case (DH A 5DH R 50kJ=mol) Figure 1 demonstrates the isothermal behavior of the process for the considered reversible reaction A ! B þ C: In this case, the reactant A is injected as a pulse using Danckwerts BCs. Here, DH A ¼ DH R ¼ 0 kJ=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.
Effects of enthalpy of reaction (DH A 50 kJ=mol, DH R 6 ¼ 0 kJ=mol) The influence of enthalpy of reaction DH R on concentration and temperature profiles is investigated in this test problem while keeping the enthalpy of adsorption as DH A ¼ 0 kJ=mol: The results are displayed in Figure 2 for two different values of DH R ¼ À20 and DH 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 DH R ¼ À20 and DH R ¼ À40: It can also be observed that an increase in the magnitude of heat of reaction DH R enhances the conversion of reactant A into products B and C. For DH R ¼ À20 the rate of conversion is 36ð%Þ which improves further to 41ð%Þ for DH R ¼ À40: The same trend can also be seen in Table 2.
Effects of enthalpy of adsorption (DH A 6 ¼ 0 kJ=mol, DH R 50 kJ=mol) In Figures 3 and 4, the effects of enthalpy of adsorption DH A are analyzed, while the enthalpy of reaction is kept as DH R ¼ 0: The results are obtained by considering three different values of enthalpy of adsorption, i.e. DH A ¼ À20 kJ=mol, DH A ¼ À40 kJ=mol (c.f Figure 3) and DH A ¼ À60 kJ=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 Joint effects of enthalpies of reaction and adsorption (DH A 6 ¼ 0 kJ=mol, DH R 6 ¼ 0 kJ=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 DH A ¼ À60 kJ=mol and the enthalpy of reaction DH R ¼ À20 kJ=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 het A ¼ 100 kJ=mol and E het A ¼ 120 kJ=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)), is less than 1(%) in all the test problems. The corresponding quantitative values given in Table 2 demonstrate the precision of numerical solutions. Effects of the model parameters g, g T , b, and b T In this subsection, the effects of the intraparticle diffusion coefficients for mass and energy (g and g T ) and mass transfer coefficients for mass and energy (b and b T ) are investigated. Figure 6 gives the results showing the effects of g and g T . It is important to mention that for the cases where the value of g is altered, the value of g T remains the same as given in Table 1 and vice-versa. The plots show that small values of g (or g 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 b (or b T ) is reduced.   Figure 9. Fully nonlinear isotherm (c.f. Equations (7) and (8)): Results of numerical calculations for isothermal and non-isothermal conditions. Here, b ref Effects of the density times heat capacity ratio of solid to liquid phases ( q S c S q L c L ) In this problem, the effects of varying density time heat capacity ratio q S c S q 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 q S c S q L c L ¼ 0:2 obtained by taking q s c s P ¼ 8KJ=1K and q L c L P ¼ 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 q S c S q L c L ¼ 1, obtained by taking q s c s P ¼ 8KJ=1K and q L c L P ¼ 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 q S c S q L c L ¼ 5:0 obtained by taking q s c s P ¼ 40KJ=1K and q L c L P ¼ 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 q S c S q L c L ¼ 1: Isotherm nonlinearities with respect to concentration Figures 9 and 10 show the results of numerical calculations for fully nonlinear isotherm given by Equation (7) with b ref j ¼ 1 for j ¼ 1, 2, and 3 and b ref 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   Figure 10. Fully nonlinear isotherm (c.f. Equation (7) and (8)): Results of numerical calculations for isothermal and non-isothermal conditions. Here, b ref as compared to the isothermal case. It can also be observed that conversion of reactant decreases further as the values of nonlinearity coefficient b ref j increase. Moreover, in the nonisothermal 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 nonlinearities.

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. ORCID Shamsul Qamar http://orcid.org/0000-0002-7358-6669