Exact equivalence at cyclic steady state between isothermal diffusion and linear driving force models for linear adsorption systems

To match the dynamics of a linear driving force model and the diffusion equation is of great practical importance in the design and optimization of adsorption separation processes. A frequency response analysis is applied to show that it is not possible to arrive at an equivalence based on a single parameter. Using this as the basis, a universal equivalence for the linear problem is constructed and closed form analytical expressions for the two parameters are derived for the sphere and slab geometries. The two parameters represent the increased effective mass transfer coefficient and a reduction in the active volume of the particle, both corresponding to the internal concentration profiles of the diffusion equation at cyclic steady state.


Notation a
Surface to volume ratio of solid, m -1 A Integration constants b Inverse of dimensionless penetration depth, Greek letters Ratio of half-cycle time and diffusion time constant Column void fraction Ratio of the effective volume of LDF equivalent model and the volume of the solid Dimensionless spatial coordinate, r

Introduction
One of Shivaji Sircar's highly cited papers (Sircar and Hufton 2000) is on "Why does the Linear Driving Force for adsorption kinetics work?" and describes how different methods can be used to arrive at the equivalence between the LDF model and the diffusion equation. This is a topic to which he returned most recently in 2018 (Sircar 2018), showing that this equivalence is a problem of significant practical impact. What may seem an academic question and a technique that may be superseded by ever increasing computational power, is in fact a very important topic even today. This is due to the fact that the design of adsorption separation processes requires the simulation to cyclic steady state (CSS) of multiple columns with a complex schedule. If an optimiser is to be used to obtain Pareto fronts in order to find conditions that minimize energy consumption and maximize productivity, 1000s of calculations to CCS are required, each in turn consisting in the simulation of 100s of cycles. Computationally efficient methods are needed. By lumping the mass transfer resistance into an LDF coefficient, the model of an adsorption column is reduced by at least one dimension when the process is either micro-or macro-pore diffusion limited, and the resulting computations are typically one order of magnitude faster. That the topic of the LDF equivalence is so important is also evidenced by the number of contributions that have been focussed on this topic since the pioneering work of (Glueckauf 1955). Glueckauf produced the equivalence from a series expansion valid for long times of the problem of diffusion in a single particle subject to a perfect step change in external concentration. This led to two expressions and the simpler one for spherical particles is what is so widely used: He then showed that this was a very good approximation to be used in the simulation of breakthrough curves and discussed also the quadratic driving force model, QDF (Sircar and Hufton 2000), of Vermeulen (1953) as an improved approximation for nearly rectangular isotherms. The QDF model can be seen as a concentration dependent LDF coefficient and therefore will not be discussed in detail here. Sircar and Hufton (2000) showed that in fact for a long column where constant pattern behaviour was established the LDF model was a good approximation for a system that could be described using the Langmuir isotherm.
The equivalence can also be derived considering that the concentration profile resulting from the diffusion equation has to be symmetric with respect to the center of the particle and is therefore a sum of even functions or a polynomial with only even powers (Villadsen and Michelsen 1978). Equation 1 is recovered using the first order approximation that the profile is a parabola, which corresponds to the most important term for long times. Higher order approximations are possible (Villadsen and Michelsen 1978; Do and Rice 1986), but these increase complexity and are not easy to generalise to multicomponent systems. When considering the dynamics of an adsorption column, the solution to a pulse experiment leads to the definition of the moments of the response. The first moment will depend only on equilibrium properties, while the second moment, i.e. the dispersion of the peak, will depend on both equilibrium and kinetic parameters (Ruthven 1984). Therefore from the first and second moment, which combined can be used to calculate the height equivalent to a theoretical plate, it is possible to arrive at general expressions that lump multiple mass transfer resistances into an equivalent LDF parameter. Also in this case, for a single diffusion resistance Eq. 1 is obtained.
Both the original derivation, based on a long-time approximation, and the parabolic concentration profile give an insight into when the simple LDF equivalence fails. When cycling of the external concentration is fast, compared to the diffusion time constant deviations are observed and in this case Eq. 1 is generalised to with Ω a function of the half-cycle time. Nakao and Suzuki (1983) were the first to investigate this numerically simulating a single particle subject to a sinusoidal change in the external concentration. Matching the amplitudes of the oscillations from the two models, they found that Ω increases rapidly for faster cycle times. Before discussing developments that followed, it is important to note that this problem has an analogue in heat transfer, which is mathematically identical. In 1942 Hausen (1942) solved the same problem and arrived at effectively the same trend, which was later rediscovered by Nakao and Suzuki 40 years later. In Hausen's more easily accessible book (Hausen 1983) all geometries are discussed as well as the limiting solutions valid for very fast cycles and the parabolic profile for slow cycles. The only minor difference is that Nakao and Suzuki give a limit value for slow cycles of 2 , corresponding to the first eigenvalue of the diffusion equation, while Hausen starts from the value corresponding to a parabolic profile, i.e. Glueckauf's 15 for the sphere. Figure 1 shows the comparison of the two approximations. Note that Hausen also gave analytical expressions that can be found in the heat transfer literature (Hausen 1983;Schmidt and Willmott 1981), which are similar to recent formulations for the spherical geometry (Hossain et al. 2016). The values of Ω for the spherical geometry that can be found in the literature show a strong enhancement in the effective mass transfer coefficient for fast cycles.
In the mass transfer literature, Alpay and Scott (1992) appear to be the first to have rediscovered Hausen's analytical limit valid for very fast cycle times. They used the solution to the diffusion equation in a semi-infinite medium confirming the numerical results of Nakao and Suzuki (1983). Additional studies have considered the effect of nonlinearity (Carta and Cincotti (1998); Zhang and Ritter 1997;Mendes et al. 1994), extension to multicomponent mixtures (Serbezov and Sotirchos 2001), and the shape of the forcing function (Sheng and Costa 1997), but any effort that is too complex has to be balanced against the direct solution of the diffusion equation.
The extension of the LDF equivalence to full processes was investigated by Raghavan et al. (1986) and is also discussed by Sircar and Hufton (2000). In both cases the correlations are qualitatively similar to the trend of Nakao and Suzuki (1983), but the values of Ω do not overlap and are lower. In the case of the heatless dryer (Raghavan et al. 1986), the coefficient Ω was found to appear to level-off to a value between 30 and 40 for very fast cycles. This was later found to be the result of an insufficient number of collocation points in the solution of the diffusion equation (Ahn and Brandani 2005). The difference between single particle dynamics and full process simulations indicates that the approach is only an approximation and an improved method can be found.
To understand the limitations of trying to match the dynamics of the single particle one has to consider that in a general flow system subject to a sinusoidal input the response will have two important properties: the amplitude of the resulting sinusoidal function at the outlet and the phase lag between the input and output functions (Stephanopoulos 1983). The approach of (Hausen 1942;Nakao and Suzuki 1983) matches only the amplitude of the oscillations, which can therefore exit the system at a time that is not the same as the full diffusion model. Dias and Rodrigues (1998) used a frequency response analysis to calculate Ω and showed that different results are obtained when matching the amplitude or the phase lag. Mendes et al. (1994) used a numerical approach and first calculated Ω to match the amplitudes and then added a phase shift to achieve a full match. Their conclusion was that there exists no universal equivalence criterion.
The aim of this contribution is to derive the analytical expressions that allow to calculate for a linear isothermal system, i.e. a system with a linear adsorption isotherm and a constant diffusion coefficient, the exact equivalence between the diffusion equation and the LDF model at cyclic steady state. This will be done introducing a second physically meaningful parameter and will be generalized to process simulations.

From column mass balance to individual particle
The heart of an adsorption process simulator is the column model. In the standard dispersed plug-flow one dimensional model (Ruthven 1984) the mass balance is written as With the associated Danckwerts boundary conditions. Here q is the average adsorbed phase concentration; c is the fluid phase concentration; u is the interstitial velocity; J is the dispersive flux; and is the void fraction of the column.
This model reduces to the well mixed adsorber when the length of the column becomes very short and the axial dispersion term dominates (Aris 1991). In this limit, the overall mass balance becomes Here F is the volumetric flowrate; V S is the volume of the solid in the column; and V F is the volume of the fluid in the column.
The dynamics of the adsorbed phase reduces in both cases to that of the individual particle when the volumetric flowrate is very high. In this limit all the particles will see the same external concentration. This shows how single particle dynamics are equivalent to column dynamics at least in one limiting condition, therefore an exact equivalence has to apply also to single particle dynamics.

Single particle exact equivalence from frequency response
We assume a sinusoidal fluctuation, Δc sin t t c , with halfcycle time, t c , of the fluid concentration and develop the solution in terms of dimensionless deviation functions (Stephanopoulos 1983), starting from Eq. 1 and The "Appendix" shows the derivation of the Laplace domain transfer functions between the forcing function, oscillating external concentration, and the average adsorbed phase concentration for the sphere and the slab.
where = t c D R 2 is the ratio of the half-cycle time and the diffusion time constant. This parameter determines if a cycle is fast. When > 1 Glueckauf's constant equivalence is a good approximation. For fast cycles, when < 1 , the equivalence becomes a function of the parameter .
With a sinusoidal forcing function, the CSS amplitude ratio and phase lag can be determined from the Laplace domain solution of the transfer function substituting s = i . This is equivalent to the traditional use of i (Stephanopoulos 1983), with the frequency due to use of the half-cycle time in the definition of the dimensionless time. The resulting amplitude ratio is the modulus of the transfer function, while the phase lag is the argument.
As discussed in Dias and Rodrigues (1998) it is not possible to match both amplitude ratio and phase lag with a single empirical constant. Figure 2 shows the single parameters obtained by matching either the phase lag or the amplitude ratios for the spherical and slab geometries. The plots show also existing correlations (Nakao and Suzuki 1983;Hausen 1983), and demonstrate that indeed the frequency response analysis can provide in an elegant manner the approximations obtained previously. Figure 3 shows the difference in phase lags that results from using the value of Ω that matches the amplitude ratio. Note that the diffusion model lags behind and this is important to understand what happens in column and process simulations. If there is a phase shift, the driving force will not be correct and this is not apparent when matching only the single particle. Process based numerical approaches (Sircar and Hufton 2000;Raghavan et al. 1986) show correlations for Ω with values that are smaller than those of Nakao and Suzuki (1983). This is most likely due to the effect of the incorrect match of the phase lag. We note that a match of the phase lag is consistent with Glueckauf's limiting value of Ω = 15 . The other important consideration is the fact that with the single parameter approach we obtain a coefficient that increases monotonically, giving the impression that faster cycle times are always beneficial and result in a significant increase of the equivalent mass transfer coefficient. The values of Ω from the phase lag are larger than the corresponding values that match the amplitudes, with an asymptote for fast cycles that is inversely proportional to , not √ . The values are diverging as the cycles become faster. Figure 4 shows snapshots of internal concentration profiles for the spherical and slab geometries, generated from the amplitude ratio and phase lag calculated from the Laplace domain transfer function of the internal profile: for the spherical geometry see Eq. 27. These profiles show two important features for fast cycles: the profiles deviate significantly from the parabolic shape, with very sharp gradients near the surface; and the internal core does not see that the external concentration is fluctuating. This is consistent with the approach of assuming a semi-infinite medium, but more importantly it shows that the true equivalence should also include the fact that the increased mass transfer coefficient corresponding to steeper internal gradients near the surface has as a counterpart a decreasing fraction of solid that is active.
To arrive at an exact equivalence between the two models, there is the need to reduce the volume of the solid and increase the mass transfer coefficient. With two equations, matching both amplitude and phase lag, the problem would appear to be numerical (Rouse and Brandani 2001), but in fact with the appropriate understanding of the physical meaning of the parameters closed analytical expressions can be obtained.
The two-parameter equivalence starts by matching the phase lag to fix the value of Ω . Then the ratio of the resulting moduli is the volume correction. This can be understood considering that when the adsorbed amounts are matched, this ratio becomes the ratio of the volumes of the diffusion and LDF models, . Note that this second parameter will not affect the phase lag, because the expression of the argument contains the ratio of the real and imaginary parts and therefore it is also equivalent to the ratio of the real and imaginary parts of the adsorbed amounts, i.e. the volumes in the calculation of the phase lag cancel.
Therefore the two-parameter LDF equivalence is given by The argument of a complex number z = x + iy is given by the inverse tangent of the ratio of the imaginary and real parts of the complex number, y x . Therefore the transfer function of the LDF model can be written as The parameters are found from Equation 10 can be used directly with any mathematical software or language that handles complex numbers and functions. For completeness we report the full analytical expressions for the sphere and slab, which can be used when only real numbers and functions are available.
The equivalence parameters reduce to simple expressions for slow and very fast cycles. For the sphere.
While for the slab.
From the physical meaning of the parameter it is possible to see that 1 b is the dimensionless penetration depth of the semi-infinite medium (Carslaw and Jaeger 1959). Figures 5  and 6 show the two parameters and their limiting values/ trends as a function of for the sphere and slab, respectively.

From single particle to full column
Now that the exact equivalence is established for the single particle we proceed to show that this is in fact the equivalence needed for any process. Consider the well mixed adsorber equations and take the Laplace transform of Eq. 4. For deviation variables the initial condition is 0.
For the system to be linear, a constant flowrate was assumed, i.e. a dilute system consistent with a linear isotherm.
This can be rearranged to obtain the transfer function for the well mixed adsorber . Therefore, if the exact equivalence is applied to the single particle, also the well mixed adsorber will be exactly equivalent.
The same applies to the full column, because each differential slice of a full column is a well mixed adsorber. Locally, at each position along the column, the mass balance in the Laplace domain can be written in terms of the fluid phase concentration and the transfer function We have therefore proven that a universal equivalence at CSS does exist for the linear LDF and diffusion models for linear adsorption systems. Figure 7 shows the concentration at 1/10 the column length for an adsorption column with the diffusion model, the exact equivalence and single parameter equivalences. The purpose of this figure is primarily to show that single parameter corrections do not reproduce correctly how the concentration waves are dampened in the column. The match of the full diffusion model and the exact equivalence serves only as a verification of the proof above. The two cannot be distinguished and the figure shows the equivalent LDF results as symbols.
To clarify how the equivalence is applied, the mass balance equation in a column is written as with Ω = 23.2 corresponds to matching the amplitude of the transfer function and this means that ΔQ ΔC is the same as the ratio obtained from the diffusion equation, but the shift changes the driving force which gives an incorrect overall result. Ω = 16 was obtained by matching the amplitude ratio at this position. It is a lower value than Nakao and Suzuki (1983) and a different value would be obtained at a different position. This is consistent with what is reported in Sircar and Hufton (2000) that single parameter correlations are dependent of the configuration.
The exact equivalence provides a very useful means of testing the numerical discretization used inside the particle for the diffusion equation, at least for the linear model. For very fast cycles Ahn and Brandani (2005) have shown how this can be optimised to achieve high accuracy without having to increase the number of discretization points by scaling two regions appropriately. Figure 8 shows the effect of reducing the number of internal collocation intervals in a cubic orthogonal collocation on finite elements solution of the diffusion equation within the column model. Figure 7 was obtained with 50 elements, which is over specified. What has not been discussed so far is the effect of the shape of the forcing function, but the exact equivalence is extended easily to this case. Any function can be represented by a Fourier series and each harmonic will have its corresponding equivalence. As a result an exact transfer function can be obtained for any shape of the inlet concentration constructing the Fourier series of the outlet concentration. What is really important though is that the higher harmonics affect progressively smaller effective volumes therefore only the primary harmonics need to be considered. This is particularly true when considering an adsorption column instead of the single particle dynamics. As the concentration wave travels over a long column the higher frequencies will be progressively dampened. An adsorption column effectively acts as a filter to high frequencies, and this indicates that for design purposes the use of the first harmonic based on the half-cycle time constant should be sufficient to define the equivalent LDF coefficients. This is also consistent with the arguments of Sircar and Hufton (2000), who pointed out that what is important is not to match an instantaneous snapshot but the average over the entire cycle.

Conclusions
General analytical expressions have been derived that allow to determine the equivalence between the LDF model and the diffusion equation for linear adsorption systems. At fast cycle times two physical parameters are needed: the first corresponds to the enhancement, Ω , of the effective mass transfer coefficient as the cycle time is decreased; the second, , takes into account the fact that only the external portion of the volume is active at fast cycle times. For relatively slow cycles the equivalence reduces to Glueckauf's expression, which in the case of a sphere gives = 1 and Ω = 15 , while for a slab Ω = 3.
The equivalence is exact at cyclic steady state. The transient to CSS is accelerated by the reduction of the effective volume. This is a further indirect advantage of this equivalence, especially in design and optimization studies. The equivalence is applicable also in heat transfer applications such as direct contact recuperators (Hausen 1983;Schmidt and Willmott 1981;Levenspiel 2014).
In systems with multiple adsorbates that are characterised by sufficiently different diffusion coefficients, different equivalence parameters are needed for each adsorbate. A full discussion of this aspect, which depends also on whether the multicomponent adsorption isotherm is linear or not, was beyond the scope of this contribution, but due to the physical nature of the equivalence parameters it is possible to tackle this problem and this will be the focus of a separate investigation. Here we note that this is a likely case in macropore diffusion limited systems, where there is a high equilibrium selectivity and in mixtures containing water in hydrophilic adsorbents.
It is important to emphasise that this equivalence should be not used to estimate the number of cycles needed to reach CSS nor used in process control applications, where the inertia of the system has to be accounted for correctly. In such cases either a separate transient for the portion of the solid that is inert at CSS has to be added or the diffusion model has to be solved, and for this optimized discretization strategies are available (Ahn and Brandani 2005) and should be used.
This contribution is in memoriam of Shivaji Sircar who has been a champion of adsorption engineering throughout his career in both industry and academia. He has inspired many researchers in the field, reminding many that it is important to try and focus research on industrially relevant problems. It is hoped that he would have appreciated the solution of a fundamental problem to which he often returned over the years.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.
Applying the boundary conditions A 1 = −A 2 = A and Therefore Equation 5 can be written in dimensionless form as Transformed in the Laplace domain and substituting the solution (27) Dividing by s the transfer function of the diffusion equation for the single spherical particle is obtained. The same procedure allows to obtain the transfer function for the slab geometry. Given that the mathematical derivation for the two geometries differs only by one boundary condition, the full derivation for the slab geometry is omitted. For the LDF model, Eq. 2 can be rewritten in terms of deviation variables and transformed in the Laplace domain giving Which rearranged gives the transfer function for the LDF model for both geometries.