Finite-volume method for industrial-scale temperature-swing adsorption simulations

We formulate a mathematical model for temperature-swing adsorption systems. A finite-volume method is derived for the numerical solution of the model equations. We specifically investigate the influence of the choice of spatial discretization scheme for the convective terms on the accuracy, convergence rate and general computational performance of the proposed method. The analysis is performed with the nonlinear Dubinin-Radushkevich isotherm representing benzene adsorption onto activated carbon, relevant for gas cleaning in biomass gasification. The large differences in accuracy and convergence between lowerand higher-order schemes for pure scalar advection are significantly reduced when using a non-linear isotherm. However, some of these differences re-emerge when simulating adsorption/desorption cycling. We show that the proposed model can be applied to industrial-scale systems at moderate spatial resolution and at an acceptable computational cost, provided that higher-order discretization is employed for the convective terms. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Separation processes account for roughly 10-15% of the world's energy consumption ( Sholl and Lively, 2016 ). One common method for removing gas phase components (such as volatile organic compounds) from a gas mixture is adsorption onto an adsorbent bed, usually activated carbon ( Yang et al., 2019 ). When the bed becomes saturated with the adsorbing species (denoted the adsorbent ), it must be regenerated before it can be reused for adsorption. This regeneration can be performed by increasing the temperature of the gas passing through the bed, which causes desorption of the adsorbent. This process scheme is referred to as temperature-swing adsorption . Steam is often used as the heat transfer medium in the desorption step. The steam consumption in an unoptimized system can be orders of magnitude higher than in an optimized one ( Schweiger and LeVan, 1993 ), indicating a need for numerical tools that allow for more efficient design, control and operation of these industrial processes.
The optimization of adsorption processes cycling between high and low temperatures requires robust numerical simulation models ( Haghpanah et al., 2013 ). Most significantly, the physicochemical nature of the adsorption process tends to produce sharp concentration and temperature fronts that propagate through the system. The sharpness of these fronts, in conjunction with the large size of the equipment in many industrial applications where adsorption columns can be many meters wide and tall, poses great challenges for many traditional numerical modeling techniques. More specifically, the spatial resolution must be extremely fine to accurately resolve a sharp front, leading to excessive computational costs for many traditional spatial discretization schemes.
In finite volume methods, several methods have been proposed to handle the sharp fronts involved in the adsorption dynamics. High-resolution total variation diminishing (TVD) methods have been particularly successful in reducing nonphysical oscillations around discontinuities while still capturing the smooth part of the solution ( Haghpanah et al., 2013 ). The TVD family of schemes thus manages to combine high-order accuracy with low-order stability, and it does so by introducing flux limiters. For example, the van Leer scheme has been shown to offer second-order accuracy in smooth regions while avoiding oversharpening effects in analyses of chromatographic problems ( Medi and Amanullah, 2011 ).
However, most previous studies have focused primarily on single-step adsorption and desorption dynamics of lab-scale systems. In addition, many models in the literature deal with pressure-swing rather than temperature-swing adsorption, so that the adsorber regeneration is accomplished by lowering the pressure rather than by increasing the temperature. Temperatureswing adsorption is in many ways fundamentally different to pressure-swing adorption, and poses its own challenges to numerical schemes, such as possible velocity oscillations due to gas expansion ( Schweiger and LeVan, 1993 ). Our main interest here is in the construction and analysis of a methodology that can robustly be applied to large adsorber beds, while maintaining enough computational efficiency to allow cycling until industrial steady operational conditions are attained (the so-called cyclic steady state ).
In the present work, we thus derive a mathematical model for temperature-swing adsorption. The model is general in its formulation, but in our own work primarily intended for application in gas cleaning of producer gas from biomass gasification, where benzene is the main adsorbate species, the bed consists of activated carbon and the process is carried out at ambient pressure ( Thunman et al., 2018 ). The emphasis is on the formulation of an accurate and robust finite-volume method for the solution of the proposed model system. We investigate the performance, convergence, computational cost and temporal advancement towards cyclic steady state for the proposed method. Our objective is in particular to assess the effect of spatial resolution on the ability of different discretization schemes to accurately resolve this approach towards cyclic steady state in temperature-swing adsorption with a non-linear isotherm. Fig. 1 illustrates a fixed-bed adsorber system, where either a process gas mixture or a hot inert purge gas flows through a packed-bed reactor system. The model equations for our temperature-swing adsorption system constitute a set of five coupled partial differential equations describing the evolution of five state variables: the gas phase temperature, the gas phase adsorbate concentration, the solid phase temperature, the amount (loading) of the adsorbed species on the solid phase and the gas phase velocity.

Main model -partial differential equations
The gas phase total continuity equation is: where ε is the packed-bed porosity (-), ρ g is the gas phase density (assumed to follow the ideal gas law) (kg/m 3 ), u is the superficial axial velocity (m/s), and q b is the loading of adsorbate species on the solid phase (kg adsorbate/kg solid phase). t denotes time and x is the axial spatial direction. The transport equation for the adsorbate species in the gas phase is: where ω b is the mass fraction of adsorbate species (-).
The adsorption of the adsorbate onto the solid adsorbent material is described by: where a is the specific surface of the adsorbent (m 2 /m 3 ), k is the mass transfer coefficient between the gas and the solid (m/s), and c b,eq is the equilibrium mass concentration given by the adsorption isotherm (kg/m 3 ). The gas phase energy balance equation is: where C p,g is the specific heat of the gas (J/kg,K), T g is the gas phase temperature (K), h is the heat transfer coefficient between the gas and the solid (W/m 2 , K) and T s is the solid phase temperature (K). The solid phase energy balance equation is: where H ad is the heat of adsorption (J/kg). It may be noted that no dispersion terms are included in Eqs.
(1) -(5) , as the Péclet numbers ( Pe = Ud /D eff , where U is the average superficial velocity, d is the diameter of the adsorption column and D eff is the effective in-bed dispersion coefficient) attained in industrial adsorption systems are typically much larger than unity, implying that advective transport dominates over dispersive transport. We will return in more detail to the rationale behind this approach in Section 5 of the paper.

Ancillary expressions -heat and mass transfer correlations
For the bed heat and mass transfer coefficients, h and k , the correlation from Gnielinski (1978) is used, where: and Nu lam = 0 . 664 P r 1 / 3 Re ε 1 / 2 (9) The convective heat transfer coefficient h can be calculated directly from Nu = hd p /k g , where d p is the equivalent diameter of the adsorbent pellet material (m) and k g is the thermal conductivity of the gas phase (W/m,K). The mass transfer coefficient k is then obtained via the Chilton-Colburn mass-heat transport analogy: or, equvialently: where D ab is the molecular diffusivity of the adsorbate in the gas phase mixture, Sc = μ/ρ g D ab is the Schmidt number, μ is the dynamic viscosity of the gas (m 2 /s), and P r = C p,g μ/k g is the Prandtl number. Schematic illustration of an industrial temperature-swing adsorption system. a) A large packed-bed reactor (zoom-in) operates as an adsorber for one or more components in the gas phase. The gas phase mixture to be treated is introduced at the bottom during the adsorption phase. b) Before the bed is completely saturated with adsorbent, the inlet flow is directed into another bed and a hot inert purge gas flow (for example steam) is flushed in reverse mode. After the bed has been regerenerated, it may again act as an adsorber.

Ancillary expressions -adsorption isotherm
A mathematical model for adsorption relies on an adsorption isotherm to predict the amount of adsorbate on the adsorbent as a function of its pressure or concentration at constant temperature. The simplest possible isotherm is the linear Henry adsorption isotherm, which assumes that the amount of adsorbate on the surface is proportional to the partial pressure in the gas. This assumption is known to deteriorate with surface inhomogeneity and adsorbate interaction beyond the low-coverage limit. Instead, the isotherm used in the current work is therefore the well-established Dubinin-Radushkevich isotherm ( Dubinin and Radushkevich, 1947;Foo and Hameed, 2010;Hu and Zhang, 2019 ): where W is the amount of adsorbate on the adsorbent at equilibrium (mg/g), β (-), E 0 (kJ/mol) and W 0 (mg/g) are adsorbentspecific constants, and where R is the universal gas constant (J/mol,K), and P 0 and P are saturation pressure (Pa) and partial pressure (Pa) of the adsorbent, respectively. This isotherm has been widely used for adsorption in microporous carbonaceous materials, such as activated carbon ( Nguyen and Do, 2001 ). It was originally developed as a semiempirical expression for adsorption isotherms based on the Polyani potential theory of adsorption ( Polyani, 1932 ). The model parameters in the isotherm expression reflect the physical and chemical adsorbate-adsorbent interaction ( Chiang et al., 2001;Terzyk et al., 2005 ). Here, we choose to employ values representative of benzene adsorption onto an industrially relevant activated carbon, resulting in β = 1 . 0 , W 0 = 500 mg/g, E 0 = 14 . 0 kJ/mol and d p = 1 cm. The saturation pressure P 0 is evaluated from the Antoine equation: where A = 4 . 01814 , B = 1203 . 835 and C = −53 . 226 are used (for T in K and P 0 in bar).

Initial and boundary conditions
At t = 0 , the state variables are assigned initial values, e.g. prescribed solid and gas phase temperatures, q b = 0 kg/kg and ω b = 0 .
Dirichlet boundary conditions are prescribed for the state variables at the inlet boundary. Since there are no dispersion terms included in the model, the often-used Neumann boundary condition of zero gradient of the state variables at the bed outlet is not needed. More advanced boundary conditions employed in adsorption/desorption cycling are described later.

Finite volume method
The finite volume method is suitable for solving the conservation laws of Eqs. (1) -(5) . The method is inherently conservative and allows for the flexibility to freely include source terms to the equations. Generally, the domain of interest is discretized into cells of (possibly) different volumes and the simplification that any quantity is constant within the cell is used. In the current implementation, the values are stored at the cell center. The equations are then integrated both in time and space over the volumes. Here, the equations are spatially one-dimensional and the cell volumes are kept equal with a constant normal area A , thus also the spatial spacing is constant and reduced to a general x .
The general conservation equation for a general variable φ is: where S represents a source term. The temporal term is discretized as: where the superscript n refers to the current time step, n + 1 to the next timestep and V i the volume of cell i .
The advective term is discretized as: with t being the time step and y denoting the time where the term is evaluated. Setting y = n gives an explicit time discretization while y = n + 1 an implicit one. In the current work, all terms are taken implicitly when possible. For example, the upwind scheme allows for an implicit discretization, hence when this scheme is used the variable is taken implicitly. The van Leer scheme does not allow an implicit discretization, hence the variable has to be taken explicitly. As both u and φ are solved for in a linear matrix system, they cannot both be implicitly taken into the matrix. Instead, for Eqs. (2) and (4) , u is taken explicitly. When performing the discretization of the advective term, the term [ A ρu φ] must be evaluated at the two faces i ± 1/2 separating the volumes at i ± 1. As the variables are stored in the cell center, they must be interpolated to the face. For the explicit term [ ρu ] i +1 / 2 that is done with linear interpolation as: The evaluation for the i − 1 / 2 face is done in the same fashion and thus only a generic face f will be discussed in the following. The procedure is the same for any face and, in fact, when assembling the matrix to be solved, it is the faces that are being iterated over, rather than the cells.
Four different interpolation schemes for the variable φ are evaluated in this work. The first-order upwind scheme is compared to three common high order/high resolution schemes: the ULTIMATE QUICKEST, van Leer and MUSCL schemes.

First-order upwind
The first-order upwind scheme simply uses the value for φ from the cell in the upwind direction of the flow, denoted C , i.e.,

TVD Schemes
The van Leer and MUSCL schemes used are TVD limiters ( Sweby, 1984 ) applied to a central difference scheme ( van Leer, 1974 ). The van Leer scheme ( van Leer, 1974 ) is : while the MUSCL scheme is: where ˜ φ x is the variable φ at a position x normalized as: and the indices C, U, D refer to the upwind, second upwind and downwind cells with respect to the face. From Eq. (17) and (19) the face value φ f can be calculated after the upwind and downwind cells have been identified.

ULTIMATE QUICKEST
The ULTIMATE QUICKEST scheme is the ULTIMATE limiter ( Leonard, 1991 ) applied to the QUICKEST scheme. As for the van Leer and MUSCL schemes, information from the upwind, second upwind and downwind nodes ( C, U, D ) is used. The scheme only limits the QUICKEST discretization under certain circumstances. More specifically, if and holds, the unconstrained QUICKEST scheme is used, and: Where CFL is the local Courant number, ( CF L = u t/ε x when u is the linear velocity). It is important to note that it is the linear velocity through the column, and not the velocity of the adsorption front, which should be used to evaluate the Courant number ( Medi and Amanullah, 2011 ). If the criterion in Eq. (20) does not hold, but instead the first-order upwind scheme is used. Otherwise, the ULTIMATE limiter is applied to the QUICKEST scheme, that is, φ f is calculated from Eq. (23) but the following limits are enforced:

Source terms
Any source term S in an equation is simply integrated as: and, as with the advective terms, the source term is taken implictly when convenient enough, otherwise it is explicitly evaluated.

Discretization of each equation
In the current work, the different schemes are only applied to the advection term of the adsorbate transport Eq. (2) . The advection term in the fluid energy equation is always discretized with the upwind scheme, while special attention is paid to the discretization of the total continuity Eq. (1) . Namely, to ensure conservation, the interpolation of ρ g u must be consistent between the different equations, and since the term is linearly interpolated in the adsorbate and energy equation, it is also discretized the same way in the continuity equation. The procedure resembles a central difference scheme, with the exception that the density and velocity are jointly interpolated as:

Solution process
The discretization is done for all conservation Eq. (1) -(5) . However, the conservation equations are coupled and several variables that are solved for appear in multiple equations, for example the term ∂ ∂t ( q b ρ b ) . To solve for all variables simultaneously, all equations are assembled in the same matrix system, enabling a coupled solution. The code for simulating the system is written in C++ with a fully object-oriented layout. The properties depending on the solution variables ( ρ g through the ideal gas law, h, k and c b,eq ) are all evaluated with values from time n .

Simulations of transport and adsorption
Numerical simulations of adsorption problems are known to be challenging and time consuming due to the non-linearity of the governing equations and the presence of strong temporal and spatial gradients ( Webley and He, 20 0 0 ). More specifically, the solutions are known to exhibit sharp fronts. These sharp fronts appear due to the advective nature of the system (high Péclet numbers) and are further significantly amplified by the nature of the adsorption process, which gives rise to a large source term on the right hand side of Eq. (3) that dominates the process. It is known that the unconditionally stable first-order upwind discretization of the convective term (physically equivalent to the tanks-in-series approach Fogler, 2016 ) results in low accuracy due to numerical smearing of the fronts ( Webley and He, 20 0 0; Haghpanah et al., 2013 ), whereas higher-order schemes such as QUICK are prone to produce oscillatory results in the presence of sharp gradients ( Webley and He, 20 0 0 ).
For this reason, we first establish the performance of the four different discretization schemes for the convective term for a pure advection case, in which a step change in the inlet concentration from c| z=0 = 0 to a constant value c in takes place at t = 0 with c = 0 ∀ x as the initial condition. No interaction with the bed material (adsorption/desorption) is activated in these simulations, i.e. ∂ ∂t ( q b ρ b ) = 0 . In the absence of physical and numerical diffusion, the system response to this step change should constitute the pure advection of a sharp and distinct front in the c -profile from the inlet to the outlet of the bed over a time equal to the mean residence time, τ = L/U, in the system. We choose to compare the concentration profiles obtained with the different discretization schemes at t = 0 . 5 τ, as illustrated in Fig. 2 .
Several important observations can be made in Fig. 2 . First of all, it is clear that the first-order upwind scheme introduces significant numerical dispersion, resulting in a failure to reproduce a sharp profile even at a resolution of 300 mesh points. ULTIMATE QUICKEST exhibits a better performance, but the MUSCL and van Leer schemes are even better, with very sharp fronts already at a resolution of approximately 100 mesh points.
We next turn our attention to the simulation of an adsorption step, in which an initially empty bed ( q b = 0 ∀ x ) is exposed to adsorbate at an inlet concentration of c in starting from time t = 0 .
Here, adsorption is activated, leading to a gradual build-up of adsorbate in the adsorbent over time. There are two main effects on the adsorbate profile in the bed in comparison to the pure advection cases: firstly, the front is typically sharpened due to the action of the adsorption source term ( Webley and He, 20 0 0 ), and secondly, the propagation velocity of the step in the profile through the bed is now governed by the adsorption dynamics, and thus significantly slower than the advection time scale τ . A comparison of the concentration profile characteristics obtained with the different discretization schemes is displayed in Fig. 3 . It can be seen from Fig. 3 that the previously identified trend from the advection test remains: the first-order upwind scheme exhibits most numerical dispersion, whereas the van Leer and  MUSCL schemes again perform similarly and superior to the firstorder scheme, and ULTIMATE QUICKEST falls somewhere in between. However, a more striking conclusion from the adsorption test is the fact that the differences between the discretization schemes are dramatically reduced. This situation may be more easily investigated by plotting the differences between the first-order upwind and the van Leer schemes on top of each other, as in Fig. 4 . The performance of the two schemes is dramatically different at both coarse and fine mesh resolutions for pure advection, whereas the performance is similar already at the coarse resolution for the adsorption case, and close to identical at the fine resolution. This very important observation is in contrast to previous findings that suggest that at least 10-15 mesh points are needed within the front itself, and that over 200 node points are needed for the first-order upwind scheme to approach the prediction of a higher-order scheme (the QUICK scheme with an adaptive smoothing algorithm) at 30 node points ( Webley and He, 20 0 0 ). The main difference between the current work and ( Webley and He, 20 0 0 ) is that we use a non-linear isotherm and physical mass transfer rates, whereas the previous investigation was performed for a linear isotherm system with artificially increased mass transfer rates to mimic the effect of the sharper fronts observed in non-linear isotherm systems. The present results suggest that numerical modeling of non-linear adsorption systems is less sensitive to numerical dispersion from lower-order spatial discretization than previously believed.
It is also of interest to study systematically the effect of the mesh resolution on the performance of the different schemes. Here, this is performed using the L 1 norm for the concentration variable, defined as: 27) where N is the number of mesh points, c x is the concentration at the node in position z = x, x is the (uniform) mesh spacing, and the subscript ref refers to the reference solution, which here is taken to be the solution obtained with the van Leer scheme with N = 10 0 0 . For mesh spacings different to that of the reference solution, linear interpolation is used to obtain the corresponding values of c ref,x . Convergence plots for the advection and adsorption test cases are shown in Figs. 5 and 6 , respectively. Due to the very close similarities between MUSCL and van Leer, the MUSCL results have been omitted from these figures.
Overall, the reported L 1 -values are up to one order of magnitude larger for pure advection than for simultaneous adsorption and advection. The slope of the L 1 norm for the first-order upwind scheme for pure advection is approximately −0.5, whereas that of ULTIMATE QUICKEST is −0.8 and that of van Leer −1.1. The convergence is thus slowest for first-order upwind and fastest for van Leer, as expected. For the cases with adsorption, the slopes change to approximately −0.9 for all schemes, indicating that also the convergence rate of the tested schemes becomes much more similar for the adsorption process.
We have also quantified the computational cost of the simulations to find out how they scale with the mesh resolution, N . The computational cost does not depend much on the discretization scheme and scales as N 2.3 for pure advection and N 2.1 for adsorption. This performance is better than the one reported by Webley and He (20 0 0) for their pressure-swing adsorption model  ( N 2.5 ). It also indicates that the matrix solver, which uses LU decomposition and scales as N 3 , contributes a significant part of the computational cost. Nevertheless, in the investigated range of N , the ancillary computations in the time evolution, being proportional to N , still contribute a non-negligible amount to the overall cost.

Numerical dispersion
In addition to quantifying the convergence rate, it is also of interest to quantify the numerical dispersion introduced by the more diffusive schemes. This numerical dispersion, D num , can be obtained by fitting the step response result to a convection-diffusion equation with the axial dispersion D a , by equating D num = D a : ∂c ∂t The fitting is performed with the non-linear least squares trust region reflective algorithm lsqnonlin wrapped around the 1-D PDE solver pdepe in MATLAB R2017b. Here, pdepe uses 10 0 0 node points and a time step corresponding to CF L = 10 −3 . The initial condition is c = 0 ∀ x and the boundary conditions are c| z=0 = c in and ∂c ∂z | z= L = 0 . The equivalent number of tanks in series, n tanks , for a tanks-in-series model with the same axial dispersion can then be estimated from the fitted D a via ( Fogler, 2016 ): It is also of interest to compare the numerical dispersion with the actual physical dispersion expected for the flow through the adsorber bed. An effective packed-bed dispersion coefficient, D eff , can be estimated from Wakao and Kaguei (1982) : For the current conditions, D AB ≈ 8 · 10 −6 m 2 /s ( Gustafson and Dickhut, 1994 ) and d p = 1 cm, resulting in D e f f ≈ 10 −3 m 2 /s. For industrial applications where smaller particle diameters are sometimes employed, D eff can go down towards 10 −4 m 2 /s. Table 1 shows the obtained numerical dispersion coefficients from the step response test simulations. Results are given for firstorder upwind and ULTIMATE QUICKEST, as the low dispersion coefficients obtained for the van Leer and MUSCL schemes cause oscillations in the solution procedure for pdepe (for example, our results for the van Leer scheme at a mesh resolution of 300 mesh points become even lower than the molecular diffusivity).
It is clear that the numerical dispersion attained at a low to moderate spatial resolution (10-30 mesh points) with a highresolution scheme is of the same order of magnitude as the expected effective physical dispersion in the bed. This observation indicates that the numerics contribute a diffusive effect that, for the relevant industrial conditions investigated here, is of the same order of magnitude as the physical effect that would have resulted from adding dispersive terms to the model equations. We thus conclude that the in-bed dispersion may to a first approximation be indirectly accounted for via the numerical dispersion, in the spirit of e.g. implicit Large Eddy Simulation for turbulent flow ( Grinstein et al., 2007 ). In addition, as the performance of different discretization schemes becomes more similar at adsorption conditions, this methodology is likely to be acceptable independent of the actual scheme chosen.

Adsorption/desorption cycles and cyclic steady state (CSS)
Industrial gas cleaning by adsorption proceeds by alternation between two modes of operation: an adsorption phase followed by a desorption phase. During the adsorption phase, the incoming gas mixture may contain non-zero concentrations of the adsorbate(s), whereby the bed is slowly saturated. At some point before breakthrough occurs (that is, before the outlet adsorbate concentration increases above some threshold value), the bed is purged by an inert gas mixture, which may enter either from the inlet or from what was the outlet during the adsorption phase. In temperatureswing adsorption, as studied here, the inert purge gas during this second phase has a higher temperature, which shifts the equilibrium adsorbate concentration that the adsorbent can hold, causing desorption to occur. The desorption process lasts until enough ad- sorbate has been removed from the bed, and the bed has been regenerated and is able to adsorb again.
If the inlet boundary conditions are altered repeatedly according to an adequately chosen cycling program in this way, the adsorber column may approach a state where the state vector (concentrations, loadings and temperature) returns identically at the end of a cycle to its value at the start of the same cycle ( Todd et al., 2001 ). This is a mathematical definition of a cyclic steady state (CSS). In adsorption modeling, studying the process from an initially clean bed until cyclic steady state is important for many reasons. One is that the concentration profiles and column dynamics may operate in entirely different mass transfer regimes as compared to the simple adsorption step onto a clean bed during cycling ( Webley and He, 20 0 0 ), which implies that CSS assessment is a necessary aspect of any investigation into numerical adsorption modeling. It is the functionality and behavior at CSS that is industrially relevant, hence the appropriate initial conditions for the simulations to be industrially relevant must reflect the CSS ( Gautier et al., 2018 ). Additionally, the cyclic steady state condition is the one that must be known for the actual process to be evaluated in a quantitative manner ( Haghpanah et al., 2013 ). Finally, the cyclic steady state can be computationally expensive to predictespecially if the process is nonisothermal, where the slow thermal dynamics of the system may disguise large effects that accumulate over many cycles -leading to that the final column profiles are appreciably different from the initially evolving ones ( Todd et al., 2001 ).
For these reasons, we study the performance of the proposed model during adsorption/desorption cycling according to the cycling program given in Table 2 . The velocity during both phases is 0.1 m/s, and the bed is regenerated from the outlet toward the inlet (reversed flow), as in common industrial practice.
In Fig. 7 , the approach towards cyclic steady state is illustrated for the bed loading in a simulation with van Leer discretization and 100 mesh points. It is seen that more than a half of the bed is completely regenerated (free from adsorbate) starting from the rightmost boundary, whereas the first 15% is completely saturated. The intermediate loading profile changes with every cycle, but it seems to converge towards a final steady condition as the number of cycles increases. In the industrial setting, knowledge about this steady condition is crucial for reaching optimal performance and minimized energy consumption, but also to assess bed ageing. In the process of cleaning the gas produced in biomass gasification from benzene, for example, low concentrations of impurities such as napthalene may adsorb irreversibly to the carbon and lower the capacity of the beds over time. The specifics of such phenomena will depend on the cyclic steady state of the beds, highlighting the importance of being able to accurately and robustly predict the industrial cycling behaviour already in the design process.
To arrive at an evaluation of the performance of a given adsorption system setup and operation, we are normally interested in a global performance indicator. One such indicator is the total integrated bed loading after regeneration, defined as: In an industrial application, one could for example employ multiobjective optimization to identify a cycling program that minimizes Q b together with the energy consumption associated with the regeneration step ( Haghpanah et al., 2013;Xu et al., 2019 ).
For the remaining discussion, we choose to contrast the performance of the first-order upwind scheme (as a representative of traditional tanks-in-series modelling for adsorption processes) with that of the van Leer scheme (as a representative of a highresolution scheme). Fig. 8 illustrates how Q b is changing as a function of the number of cycles in the adsorption/desorption process for these two different discretization schemes at two different mesh resolutions (30 and 100 mesh points). Interestingly, it can be observed that the convergence towards the cyclic steady state is slower in the first-order upwind simulation with 30 mesh points than in the other three cases. Moreover, the two van Leer simulations converge towards the same value of Q b at cyclic steady state, whereas this global prediction from the first-order upwind simulations changes with spatial resolution. In other words, despite the surprisingly good performance of first-order upwinding in the adsorption single step test cases, the performance in this more realistic and comprehensive test case identifies significant advantages of using higher-order schemes. However, even if the van Leer scheme allows for a fast and efficient capturing of the global performance indicator already at a resolution of 30 node points, a detailed look into the final bed loading profiles at cyclic steady state reveals that the distribution of the adsorbate is still changing somewhat with the spatial resolution (cf. Fig. 9 ). For both the van Leer and the first-order upwind schemes, the bed loading profile attains a sharper corner at the end of the fully loaded section, along with a decrease in the loading in the region closest to the fully regenerated part of the bed.
The computational cost for an adsorption/desorption cycle scales approximately as N 2.2 , indicating that the performance is on par with that observed for the advection and adsorption step test cases. At a resolution of 30 mesh points, it takes approximately one second of real time to advance the simulation time ten seconds.
The overall analysis points to that approximately 30 mesh points are sufficient to resolve the convergence towards a cyclic steady state with an accuracy similar to that obtained with much finer spatial resolution. This observation is in line with the conclusion from Haghpanah et al. (2013) for investigations of CSS attainment in pressure-swing adsorption simulations. Moreover, at this level of resolution, the numerical dispersion present in the numerical model is of the same order of magnitude as the physical axial dispersion in the bed, further underlining the applicability of this mesh resolution for simulations of industrial adsorption units.

Conclusions
We present a finite-volume method for numerical simulations of industrial temperature-swing adsorption processes. It is shown that for an adsorption process characterized by a nonlinear Dubinin-Radushkevich isotherm, valid for benzene adsorption onto activated carbon, the effect of the spatial discretization scheme on the attainable accuracy in the resolution of sharp concentration fronts is much less pronounced than for pure advection. Both the accuracy and convergence rate of the first-order upwind scheme (or, equivalently, a tanks-in-series approach) become similar to that of higher-order TVD schemes (e.g. van Leer) for simple adsorption step test cases. However, during adsorption/desorption cycling towards a cyclic steady state, the higher-order scheme outperforms the first-order upwind one in both the rate of convergence of a global performance indicator, as well as in the sensitivity of the spatially resolved profiles to the mesh resolution em-ployed. Robust and accurate predictions of the cyclic steady state are required in the design of industrial temperature-swing adsorption systems if the long-term behaviour is also to be properly assessed.
The numerical dispersion introduced by the ULTIMATE QUICK-EST scheme (as a representative of a typical higher-order scheme) at a moderate mesh resolution (30 node points) is shown to be of the same order of magnitude as the expected axial dispersion in the flow through the packed bed, so that the proposed model formulation without dispersion terms may produce solutions that implicitly account for physical in-bed dispersion. As the use of as few as 30 node points is shown to be acceptable also for adsorption/desorption cycling with a higher-order scheme, this level of spatial resolution can be recommended for numerical simulations of industrial adsorption processes with the proposed model.
In conclusion, the presented model is shown to be able to reproduce the characteristics of packed-bed adsorption systems during industrial adsorption/desorption cycling at moderate spatial resolution and at good computational efficiency. The model therefore fulfills the requirements necessary for robust and accurate numerical simulations of non-linear large-scale industrial temperature-swing adsorption systems, such as commonly used gas cleaning equipment in biomass gasification.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.