Modelling and Simulation of a Packed Bed of Pulp Fibers Using Mixed Collocation Method

Abstract

A convenient computational approach for solving mathematical model related to diffusion dispersion during flow through packed bed is presented. The algorithm is based on the mixed collocation method. The method is particularly useful for solving stiff system arising in chemical and process engineering. The convergence of the method is found to be of order 2 using the roots of shifted Chebyshev polynomial. Model is verified using the literature data. This method has provided a convenient check on the accuracy of the results for wide range of parameters, namely, Peclet numbers. Breakthrough curves are plotted to check the effect of Peclet number on average and exit solute concentrations.

1. Introduction

The removal of the soluble material, occupying the interstitial space between the particles and in the pores of the particles, from a saturated packed bed of particles is carried out by the introduction of a solvent, for example, water or weak wash liquor flowing through the bed. Solute removal is associated with diffusion like dispersion of the solvent in the direction of flow, known as longitudinal dispersion. The mechanics involved are the sum of displacement of the fluid containing solute by movement of water plug controlled by fluid mechanics, dispersion due to back mixing, diffusion due to concentration gradient, and adsorption-desorption due to relative affinity of various solutes towards the particle surface. The core problem is the prediction of the behavior of initially sharp interface between the liquids having identical dynamical and kinematical properties. This problem is of considerable practical importance, for example, in determining the efficiency of solvent utilization or filtrate recovery in washing of filter cakes. Mathematically, such chemical engineering processes can be best described by two point boundary value problems as follows: 𝑎₁𝜕²𝑐𝜕𝜉²+𝑎₂𝜕𝑐𝜕𝜉=𝑎₃𝜕𝑐𝜕𝜏+𝑎₄𝜕𝑞,()𝜕𝜏𝜉,𝜏∈𝐷, (1) where 𝐷=Ω×(0,𝑇]=(𝑎,𝑏)×(0,𝑇], along with the boundary conditions as follows: 𝑏₁𝑐+𝑏₂𝜕𝑐𝜕𝜉=𝑘₁𝑏at𝜉=𝑎,₃𝑐+𝑏₄𝜕𝑐𝜕𝜉=𝑘₂at𝜉=𝑏. (2)

A great deal of effort has been applied to compute efficiently the solution of transient partial differential equations (1) and (2) analytically [1–7] and numerical algorithms such as Hermite radial basis function interpolation numerical scheme [8], two-stage Lie-group shooting method [9], finite difference method [10–12], spectral collocation method [13], Sinc differential quadrature method [14], orthogonal collocation method [10, 15–17], fitted mesh collocation method [18], a novel numerical scheme [19], Galerkin/Petrov Galerkin method [20–22], orthogonal collocation on finite elements [23–25], factorized diagonal Padé approximation [26], spline collocation methods [27, 28], variational iteration method (VIM) [29, 30], homotopy analysis method [31, 32], homotopy perturbation method (HPM) [33–35], and energy balance method [36].

In finite-difference method, however, the solution of the system is very unstable and requires strict selection of step size. Nevertheless, the accuracy of numerical solution is not so high. The discretization of even a few PDEs by the method of lines can lead to an extremely large system of ODEs; the numerical solution of which may have severe cost and storage implications. The Galerkin method is unsuitable for large systems because the integration process becomes very tedious. The method of OCM here is preferred over HPM and VIM, due to its simplicity and easy adaptability to computer codes. In OCM the choice of trial function is comparatively simple over other weighted residual methods. Due to its easy adaptability to computer codes, Finlayson [37] has also preferred OCM over other weighted residual methods and also over HPM and VIM.

The present paper emphasizes the use of mixed collocation method. It is a variation of the weighted residual method. It was developed to solve the models related to the transport phenomenon by Villadsen and Stewart [15]. The collocation points are the zeros of the orthogonal polynomials taken from the expansion of the trial function, enforcing the residual to vanish at the collocation points. It is a very simple, elegant, and particularly useful technique for chemical engineering problems. It gained popularity due to less computational time and easy adaptability to the computer programs. Its accuracy is comparable with the analytical one.

2. Theoretical Formulation

Consider a thin slice of a packed bed consisting of homogeneous symmetrical porous particles, as shown in Figure 1, through which filtrate or wash water flows. The mass balance equation for the packed bed can be described by an axial dispersion model as 𝐷_𝐿𝜕²𝑐𝜕𝜉²=𝑢𝜕𝑐+𝜕𝜉𝜕𝑐+𝜕𝜏1−𝜀𝜀𝜕𝑞.𝜕𝜏 (3)

Figure 1 

A simple shell balance.

The adsorption equilibrium relationship, describing the intraparticle solute concentration as a function of the external solute concentration, is linear; that is, 𝑞=𝑘𝑐. The inlet boundary condition is assumed to be Danckwerts as 𝑢𝑐−𝐷_𝐿𝜕𝑐𝜕𝜉=𝑢𝐶_𝑆,at𝜉=0,∀𝜏≥0. (4) At the outlet, no flux is assumed; that is, there is no loss of solute from the bed through the plane at which the displacing liquid is introduced as follows: 𝜕𝑐𝜕𝜉=0,at𝜉=𝐿,∀𝜏0. (5) The initial condition is 𝑐=𝑞=𝐶_𝑖,at𝜏=0,∀𝜉. (6)

The equations (3)–(6) are converted into the dimensionless form using dimensionless parameters mentioned in nomenclature. The following equations are obtained: 1𝜕Pe²𝐶𝜕𝜉^∗2=𝜕𝐶𝜕𝜉^∗+𝜀𝜀+(1−𝜀)𝑘𝜕𝐶𝜕𝜏^∗+1−𝜀𝜀+(1−𝜀)𝑘𝜕𝑄𝜕𝜏^∗;1𝐶−Pe𝜕𝐶𝜕𝜉^∗=0,at𝜉^∗=0,∀𝜏^∗≥0;𝜕𝐶𝜕𝜉^∗=0,at𝜉^∗=1,∀𝜏^∗≥0;𝐶=𝑄=1,at𝜏^∗=0,∀𝜉^∗. (7) The exit solute concentration (𝐶_𝑒) is calculated at 𝜉^∗=1. The average solute concentration (𝐶_av) is calculated by the integration over the bed cross section as 𝐶_av=¹₀𝐶𝑑𝜉^∗. (8)

3. Numerical Procedure

Details about the collocation method are available elsewhere [15, 24, 25, 37]. In this method, the trial function simply converts the system of partial differential equations into a set of differential algebraic equations (DAEs). The system obtained can be solved by using any computer subroutine.

After applying the mixed orthogonal collocation on (7), following set of DAEs is obtained: 1Pe^𝑁_𝑖=0𝐵_𝑗𝑖𝐶_𝑖−^𝑁_𝑖=0𝐴_𝑗𝑖𝐶_𝑖=𝑑𝐶_𝑗𝑑𝜏^∗𝐶,𝑗=1,2,…,𝑁−1;₀−1Pe^𝑁_𝑖=0𝐴_𝑗𝑖𝐶_𝑖=0,𝑗=0;_𝑁𝑖=0𝐴_𝑗𝑖𝐶_𝑖𝐶=0,𝑗=𝑁;_av=^𝑁_𝑖=0𝑊_𝑖𝐶_𝑖. (9) The first- and second-order derivatives (𝐴_𝑗𝑖s and 𝐵_𝑗𝑖s) are calculated by the method given by [24, 37]. The set of DAE’s is solved using MATLAB with ode15s system solver. This program uses the backward differentiation formula to solve the stiff system of differential algebraic equations.

3.1. Collocation Points

In the present study, zeros of shifted Chebyshev polynomials have been used as collocation points, which tend to minimize the error as proposed by Fan et al. [38] and yield satisfactory results for the unsymmetrical boundary value problems. The equidistant spacing in the collocation points is usually not preferred due to Runge divergence phenomenon (Villadsen and Stewart[15]). The 𝑁+1 interpolation points are chosen to be the extreme values of an 𝑁th order shifted Chebyshev polynomial. The zeros of shifted Chebyshev polynomials are calculated using Gauss-Lobatto quadrature formula as cos⁻¹2𝜉^∗_𝑗=𝜋−1𝑁−𝑗𝑁,𝑗=0,1,…,𝑁; (10) five to nineteen zeros are used as collocation points to solve the model.

3.2. Convergence Criterion

The convergence of the method for steady state is checked by applying the following formula: 𝐿^∗=Δ𝑥_ℓ²𝑀⁻¹𝐶,with0<𝐿^∗<1, (11) where 𝑀 is the lumped coefficient matrix and 𝐶 is the column matrix of the collocation solutions at any time 𝑡. For large values of the parameter Pe, the values of 𝐶 may fluctuate. This indicates that more elements should be inserted to make 0<𝐿^∗<1. With the increase in number of elements, the method converges asymptotically (Arora et al. [24]).

4. Verification of the Model

The numerical results are verified by comparing with the analytic ones given by Brenner [1]. A good agreement is found between these two results. The comparison is presented for 17 interior collocation points in Figure 2 for Pe 0, 32, and 80. A deviation of 0%–5% is found between the numerical results and analytic results for these Peclet numbers. However, the magnitude of error has decreased with the increase in collocation points, which is discussed in the next section.

Figure 2 

Comparison between analytic and numerical results for different Pe.

4.1. Error Calculation

The relative error for exit solute concentration is calculated using the formulae (𝐶_ex−𝐶_nm)/𝐶_ex, where 𝐶_ex is the analytic value reported by Brenner [1] and 𝐶_nm is the numerical value of the present case for exit solute concentration. The relative error is plotted for different degrees of freedom for Pe=40 in Figure 3. For small degree of freedom, that is, for three interior collocation points, the error increases very sharply, but as the number of collocation points increases up to seventeen, the error becomes negligible. After seventeenth degree of freedom, the concentration profiles overlap with each other and yield identical results. It justifies the fact that with the increase in the number of interior collocation points after 15, no significant change in the relative error is observed. Also for small or medium range of Pe, the relative error is approaching to zero as the analytic and numerical results are matching perfectly in Figure 2.

Figure 3 

Relative error for various collocation points.

5. Results and Discussion

The analytic solutions of similar type of axial dispersion models using different boundary conditions are available in the literature [1, 3]. These analytic solutions have constituted a basis for the comparison with the numerical solutions obtained from the present study. The results are presented in the form of breakthrough curves for various degrees of freedom for a wide range of Peclet numbers varying from very small to very large.

5.1. Effect of Peclet Number on 𝐶_𝑎𝑣

Figure 4 shows the effect of different Pe varying from 0 to 200 on average solute concentration. For Pe=200, the concentration profile is very peaked and the convergence rate is very fast. However, for small Pe, the concentration profile is much curved and the convergence rate is very slow. It is also very clear from Figure 4 that for large Pe, the behavior of concentration profile is more or less the same, which ultimately shows that for higher Pe the effect of axial dispersion coefficient is the least.

Figure 4 

Effect of different Peclet numbers on 𝐶av.

5.2. Effect of Peclet Number on 𝐶_𝑒

Figure 5 shows that for very small Pe, the exit solute concentration profile is very steep, indicating that large time period is required for the solute to diffuse out from the particle pores. As the Pe increases the curve follows a Gaussian shape (Figure 6). For medium range of Pe, however, the time for washing is large resulting in slow convergence of concentration profiles, but it is still in the range of acceptability.

Figure 5 

Effect of very small Peclet numbers on 𝐶𝑒.

Figure 6 

Gaussian behavior of 𝐶𝑒 curves for various Pe.

Figure 7 shows the effect of large Pe on exit solute concentration. It is observed that as Pe increases, the time for washing decreases sharply. It indicates the fact that for large Pe, axial dispersion coefficient becomes smaller, as a result, more solute will diffuse out from the particle pores.

Figure 7 

Effect of large Peclet numbers on 𝐶𝑒.

5.3. Limiting Cases

Perfect mixing, that is, when Pe→0. In such situation each differential element of the solvent introduced into the bed instantaneously mixes with the contents of the bed, and an equal volume of the fluid is displaced from the bed. The bed behaves like a perfect mixing chamber. For very small Pe, diffusion plays a dominant role since interstitial velocity is very small, resulting in slow convergence of the concentration profiles as shown in Figure 8. This type of situation is not considered ideal for industrial practice.

Figure 8 

The case of perfect mixing.

Perfect displacement, that is, when, Pe→∞. In such a situation, (3) reduces into a partial differential equation of order one. The role of interstitial velocity becomes more important than that of the diffusion. The solution profile becomes very broaden, and time for washing falls very rapidly as shown in Figure 9. In such a situation the initial contents of the bed are pushed out in a piston-like fashion by the displacing fluid. But this situation is not practical as the diffusion can never be zero or the interstitial velocity can never be infinite.

Figure 9 

The case of perfect displacement.

Ideally, it is not possible to remove all the soluble impurities from the packed bed. Therefore an optimum range of Pe should be followed in the industry to keep a balance between the removal of impurities and time of washing. Al-Jabari et al. [20] have also shown that flow characterization of homogeneous packed beds for bed void fraction of about 0.9 can best be represented by Pe=40. This flow is intermediate between the cases of perfect displacement (Pe=∞) and perfect mixing (Pe=0).

5.4. Effect of Peclet Number on Packed Bed

Effect of Pe on exit solute concentration at different locations in the mass transfer zone (MTZ) for constant 𝜏^∗ is presented in Figure 10. For small Pe, more back mixing effect is obtained due to increasing effect of axial dispersion coefficient keeping the interstitial velocity and cake thickness constant. For higher Pe, steeper mass transfer zone is observed. This figure indicates that for small Pe, large time is required for the solute to diffuse out of the particle pores.

Figure 10 

Effect of different Pe on mass transfer zone for fixed value of 𝜏∗.

6. Conclusions

Numerical solution of a packed bed comprising of porous semisolid cylindrical particles is presented. A robust and convenient technique of mixed collocation method shows good agreement with the analytic results. The method is computationally efficient as it consumes less computation time. It is also observed that the efficiency of the results increases with the increase in collocation points. For efficient washing operations, the Pe should lie within the range of 20 to 40. A similar type of nonlinear problem can also be solved in exactly the same manner, as the condition of linearity is of no special use.

Nomenclature

𝑐:	Concentration of solute in the liquor, kg/m3
C:	Dimensionless concentration, (𝑐−𝐶𝑆)/(𝐶𝑖−𝐶𝑆)
𝐶𝑖:	Inlet solute concentration, kg/m3
𝐶𝑆:	Solute concentration in the wash liquor, kg/m3
𝐷𝐿:	Longitudinal dispersion coefficient, m2/s
𝑘:	Mass transfer coefficient, m/s
𝐿:	Thickness of the cake, m
Pe:	Peclet number, dimensionless Pe=𝑢𝐿/𝐷𝐿
𝑞:	Concentration of solute on the particles, kg/m3
𝑄:	Dimensionless concentration, (𝑞−𝐶𝑆)/(𝐶𝑖−𝐶𝑆)
𝑢:	Interstitial velocity through bed, m/s
𝜀:	Porosity of cake, dimensionless
𝜉:	Distance from point of introduction of solvent, m
𝜉∗:	Dimensionless distance, 𝜉/𝐿
𝜏:	Time, s
𝜏∗:	Dimensionless time, 𝜏∗=𝑢𝜏𝜀/𝐿(1−𝜀)𝑘.

Acknowledgments

Ishfaq Ahmad Ganaie is thankful to NBHM, Mumbai, for JRF under the Research Project (no. 2/48(14)/2009/R&D-II/2806). Shelly Arora is thankful to UGC for providing financial assistance in the form of Research Project no. 41-786/2012 (SR).