- Email: [email protected]
Numerical solution of hyperbolic models of transport in bidisperse solids
Computers and Chemical Engineering 24 (2000) 1981 – 1995 www.elsevier.com/locate/compchemeng
Numerical solution of hyperbolic models of transport in bidisperse solids F. Liu a, S.K. Bhatia a,*, I.I. Abarzhi b b
a Department of Chemical Engineering, The Uni6ersity of Queensland, Brisbane, Qld 4072, Australia Engineering Thermophysics Institute, Ukranian Academy of Sciences, Zhelyabo6a 2a, 252057 Kie6, Ukraine
Received 17 August 1999; received in revised form 17 April 2000; accepted 22 May 2000
Abstract This article considers the numerical solution of generalized transport models for adsorption that include the presence of finite local relaxation times. Such models arise in adsorption problems where the small diffusivity in micropores can lead to local nonequilibrium even at the macroscopic process time scale. We consider the important case of a bidisperse-structured adsorbent, the transport in which is a problem of long-standing interest to chemical engineers. In the formulation considered, the microparticle model involves a non-linear hyperbolic differential equation with source terms, and in the macroparticle model the differential equation involves a singularly perturbed parabolic diffusion problem. There are some difficulties in solving these equations due to the hyperbolic nature of the microparticle transport, conventionally assumed to obey the Fickian model leading to a parabolic diffusion equation. In this paper, an efficient numerical technique is presented to simulate these processes. A combination of the upwind method and fitted mesh collocation method is applied to solve the problem. Numerical solutions were found to match the analytical solution of the traditional model when the adsorption isotherm is linear, for macropore diffusion control or micropore diffusion control. Simulations of the adsorption and desorption dynamics are also presented for a Langmuir isotherm. The numerical scheme offers a more generalized alternative that can be used for both the model forms, with and without consideration of finite local relaxation times. © 2000 Elsevier Science Ltd. All rights reserved. Keywords: Adsorption; Bidispersed particles; Upwind method; Fitted mesh; Collection method; Hyperbolic
1. Introduction The description of adsorbate transport in microporous solids is an important component of the design and optimization of the adsorption processes. Generally microporous adsorbents, such as activated carbon or pelletized zeolites, are bidisperse and have both micropores as well as macropores. The description of transport processes in such bidisperse solids is difficult due to the complex and largely unknown nature of the pore network, and a variety of structural idealizations have been proposed, based on which approximate methods for diffusion and adsorption have been developed. The simplest and the most popular of these idealizations considers the solid particle as comprising an aggregate of microporous spherical grains, with the intergrain macropores providing the channels for intraparticle transport. The modeling of transport in such a solid was first presented by Ruckenstein, Vaidyanathan and Youngquist (1971), and considered the diffusion processes over two length scales: the particle scale (macropore diffusion) and the local/microscopic scale (micropore diffusion). Since this development, numerous studies have used this model for interpreting adsorption/desorption dynamics in porous adsorbents (Bhatia, 1987); however, some workers (Gray & Do, 1990) have added further arbitrary mechanisms representing particle scale transport of
* Corresponding author. Tel.: +61-7-33654263; fax: + 61-7-33654199. E-mail address: [email protected] (S.K. Bhatia) 0098-1354/00/$ - see front matter © 2000 Elsevier Science Ltd. All rights reserved. PII: S0098-1354(00)00596-2
1982
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
adsorbed species in order to improve the fitting of the experimental data. A sounder basis for particle scale adsorbate transport has been subsequently provided (Bhatia, 1997; Bhatia, Gray & Do, 1991) by considering the asymmetric macropore concentration profile over the grain surface, and the resulting two-dimensional transport in the grain. This leads to a through flux across the grains, and facilitates particle scale adsorbate transport, though at the cost of considerable additional complexity. Most recently Abarzhi (2000) has provided a further development in the analysis of adsorbate transport by considering the effect of local relaxation times in the macropore and micropore networks. The classical model for micropore diffusion in the grains (Ruthven, 1984) involves use of the Fickian flux model Nm = − Dm(Cm)9Cm
(1)
with a concentration dependent diffusivity following the Darken relation. As indicated by Luikov (1980) and Abarzhi (2000) such a model neglects the effect of local relaxation times in the diffusion process, and is not exact. While negligible at the single pore level the local relaxation time has been considered important for the intragrain pore network as a whole in some cases. When finite local relaxation time is considered, the resulting microparticle flux model is non-Fickian, and has the form a
(Nm +Nm = − Dm(Cm)9Cm (t
(2)
in which a represents a local relaxation time scale, that may be expected to be inversely proportional to the microparticle diffusivity. In the presence of a concentration-dependent diffusivity, as is the case for adsorbate transport, a would be inversely proportional to Dm(Cm). In general it may be anticipated that the local relaxation time is negligible in the intergrain macropores, where the diffusivities are three to four orders of magnitude larger that those in the micropores. Such a model would be most relevant to molecular sieving carbons and zeolites where the extremely small micropore diffusivities may lead to finite relaxation times. Abarzhi (2000) has shown that the local non-equilibrium is important at times comparable with the local relaxation time, although near equilibrium the flux law approaches the Fickian formulation in Eq. (1). To enable the application of adsorption models involving finite local relaxation times to actual experimental data it is pertinent that the numerical solution of such models be explored, as analytical solutions are extremely complex and not possible for the case of concentration-dependent diffusivities. The solution of such models must clearly be a precursor to their application, and studies in this direction should facilitate further research along these lines and experimental validation of the models. To this end the present article is devoted to the numerical solution of the transport model for adsorption in bidisperse solids, with the intragrain flux represented by Eq. (2) in place of Eq. (1). In the present case, combination of Eq. (2) with the balance equation (Cm = − 9 · Nm (t
(3)
for the grain micropores, provides for a hyperbolic conservation law as opposed to the parabolic model arising from the use of the classical Fickian approach. This requires specialized and complex solution techniques, a problem that is compounded by the coupling of the grain balance equation with the macropore balance equation, which is parabolic in nature. In addition, the latter often involves a singularly perturbed problem (Liu & Bhatia, 1999). In the present article, we propose a combination of an upwind method (Bermudez & Vazquez, 1994) for the intragrain hyperbolic transport and the fitted mesh collocation technique (Liu & Bhatia, 1999) for the macropore transport, along with the use of the differential/algebraic system solver package (DASSL) developed by Petzold (1982). The latter considerably improves the speed of the upwind methods, which conventionally rely on finite difference approximations. Although upwind methods are not new, their use in hyperbolic mass-transport problems has not received significant attention. The present article presents an initial study in this direction, demonstrating the success of these techniques for such problems. For sufficiently small relaxation times it is demonstrated that the solution converges to that of the parabolic model as expected, so that the technique proposed is a generalization valid for both the model forms.
2. Mathematical model We consider simultaneous diffusion and adsorption in a spherical sorbent pellet comprising an aggregate of small
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1983
spherical microparticles of uniform size. The mass balance equation for the distribution of sorbate in the macrosphere is (Liu & Bhatia, 1999): s
n
(C*M 1 ( (C*M = 2 h2 +3(1 −s)g[N*(t, h, j)]j = 1 m (t h (h (h
(4)
with boundary conditions h= 1:
(C*M = Bi(C*e − C*M ) (h
(5)
h= 0:
(C*M = 0. (h
(6)
A differential mass balance equation for the distribution of sorbate in the microsphere is
n
(C*m 2N*m (N*m = −g + , (t j (j
(N*m (C*m *) = − N*m +D*(C m m (t (j
b*a*(C*) m
(7)
n
(8)
with boundary conditions j= 1: C*= I*(C *M ), m m j =0:
(9)
(C*m = 0, N*m = 0 (j
(10)
where I *m is a dimensionless adsorption isotherm. I *(C *M)= C *M when the sorption isotherm is linear, and I *(C *M) = m m C *M(1+K *)/(1 +K *C *M), when the sorption isotherm follows the Langmuir form. Here b*= Dpsab/R 2p. m m For the adsorption cycle, the initial condition is t= 0: C*M =C*=0, C*e = 1, N*m = 0 m
(11)
while for the desorption cycle, the initial condition is t =0: C*M =C*=1, C*e = 0, N*m = 0. m
(12)
The parameter s is the fraction of equilibrium uptake in macro-voids and the quantity (1− s) is that in the microsphere. Parameter g is called the dynamic parameter because it is a measure of the relative rates of diffusion in the macropores and diffusion in the micropores. We now rewrite the above non-dimensional model in the form:
n
1 ( (C*M (C*M =a1 2 h2 +a4[N*] m j = 1, (ts h (h (h 1
(13)
h=1:
(C*M = Bi(C*e − C*M ), (h
(14)
h= 0:
(C*M =0 (h
(15)
and (C*m (N*m N* +a2 = − 2a2 m, (ts 1 (j j
(16)
(N*m a3D*(C *) a3 m m (C* m + =− N*, m (ts 1 b*a*(C*) (j b*a*(C*) m m
(17)
j =1: C*= I*(C *M ), m m
(18)
j =0:
(C*m =0, N*= 0 m (j
with the model parameters having two alternate definitions: Type A: a1 =1/s1, a2 =sg/s1, a3 =s/s1, a4 =1, ts 1 = tss1.
(19)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1984
Type B: a1 =1, a2 =sg, a3 =s, a4 =s1, ts 1 =ts where s1 =3(1−s)g, ts =t/s. The fractional uptake is evaluated by determining the amount of sorbate that has entered the particle after changing the external concentration. The fractional uptake (Liu & Bhatia, 1999) is given by FA =
M( t M(
with
(20)
&
M( t =3 s
1
& &
h 2C*M dh+ 3(1 −s)
0
1
h2
0
1
n
j 2C*djdh , m
0
M( =s + (1− s)I*(1) m
(21) (22)
3. Upwind method for hyperbolic conservation laws with source terms in microparticles The differential equations in the microparticle involve nonlinear hyperbolic Eqs. (16) and (17) with source terms. The section deals with the extension of some upwind schemes to hyperbolic systems of conservation laws with source terms. Upwind methods for hyperbolic conservation laws have undergone considerable development in recent years, leading to a variety of implementation techniques such as flux-difference and flux-splitting methods (Bermudez & Vazquez, 1994). Moreover, the foundations of the theoretical work are well established. A review of these has been given by LeVeque (1990) and Godlewski and Raviart (1991). In this section, the upwind methods are only used in spatial descretization. The partial differential equations (Eqs. (16) and (17)) in the microparticle are reduced to an ordinary differential equation system. Let us consider the hyperbolic system of conservation laws with source terms (v(j, ts 1) (F(v(j, ts 1)) =G(j, v(j, ts 1)) + (ts 1 (j where v(j, ts 1)=
(23)
C*(j, ts 1) m , ts 1) N*(j, m
(24)
Á  a2N*m à à F(v)= Ãa D*(C*) Ã, à 3 m m C*m à Äa*(C*)b* Å m
(25)
Á N*  −2a2 m à à j G(j, v)= à Ã. a3 Ã− N*m Ã Ä a*(C*)b* Å m
(26)
To discretize in space we use the finite volume method. The domain is discretized by taking a uniform mesh {ji =(i− 1)Dj, Dj = ji −ji − 1, (i =1, …, N + 1)} and each interval is divided into two to form a cell Ci = (ji −Dj/2, ji + Dj/2. Let us consider a discretization in which the notes do not lie on the boundary. The approximate solution W(j, ts 1) to v(j, ts 1) at ts 1 can be considered to be a piecewise constant function given by
W(ji, ts 1)=Wi (ts 1) for j Ci = ji −
Dj Dj , ji + . 2 2
(27)
To obtain the spatial approximation of Eq. (23), we integrate this over the cell Ci. This leads to the following scheme: dW Fi + (1/2) −Fi − (1/2) + =G(ji, W(ji, ts 1)) = Gi dts 1 Dj
(28)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1985
where Fi 9 (1/2) and Gi denote approximations of F(W(ji 9 (1/2))) and (1/Dj) Ci G(j, W)dj, respectively. Upwind schemes are obtained by replacing Fi 9 (1/2) by a numerical flux function f. In a similar way we introduce a numerical source function c to upwind the source term. More precisely, we take for each i: Fi + (1/2) = fi (Wi, Wi + 1),
(29)
Fi − (1/2) = fi (Wi − 1, Wi ),
(30)
Gi = c(ji − 1, ji, ji + 1, Wi − 1, Wi, Wi + 1).
(31)
This represents an upwind scheme corresponding to numerical fluxes F(Wi − 1) +F(Wi ) 1 − Q(Wi ) (Wi − Wi − 1), 2 2 F(Wi ) +F(Wi + 1) 1 − Q(Wi ) (Wi + 1 − Wi ), fi (Wi, Wi + 1)= 2 2
fi (Wi − 1, Wi )=
(32) (33)
where Q is a matrix characteristic of each Q-scheme having a continuous dependence on Wi. Some notations in Bermudez and Vazquez (1994) will be used in this paper. For any diagonalizable matrix A, let L be the diagonal form of A given by its eigenvalues, then A = XLX − 1
(34)
where X is the matrix of the eigenvetors of A. We recall that A = X L X − 1
(35)
and L ij = li dij.
(36)
Q-scheme equals to the Jacobian matrix A of the flux valued at the value of Wi : Q(Wi ) = A(Wi ) =X L X − 1 =diag( l1(Wi ) , l2(Wi ) )
(37)
where diag ( l1(Wi ) , l2(Wi ) ) is a diagonal matrix and l1(Wi )=
'
a2a3D*m,i , a*b* i
l2(Wi ) = −
'
a2a3D*m,i a*b* i
(38)
*(j Here D*mi =D*(C m m i, tsi ))and a* i = a*(C*(j m i, ts 1)) The discretization of the flux term is Fi + (1/2) − Fi − (1/2) 1 1 = [fi (i, i +1) −fi (i −1, i )]= Dj Dj Dj
' '
Á N*(i+ 1, j )−N*(i 1 a2a3D*m,i  m m −1, j ) − [C*(i+ 1, j )−2C*(i, j )+ C*(i− 1, j )] à m m m à a2 2 2 a*b* à à a D * C *(i+ 1, j ) −C *(i −1, j ) 1 a2a3D*m,i 3 m,i m m à − [N*(i+ 1, j )− 2N*(i, j )+ N*(i− 1, j )] à m m m Ä a*b* 2 2 a*b* Å i i
(39)
The numerical source function c will be taken of the form c(ji − 1, ji, ji + 1, Wi − 1, Wi, Wi + 1) = c L(ji − 1, ji, Wi − 1, Wi ) + c R(ji, ji + 1, Wi, Wi + 1)
(40)
where left and right numerical source functions are given by 1 c L(ji − 1, ji, Wi − 1, Wi ) = [I+ Q(Wi ) Q − 1(Wi )] 2 G. (ji − 1, ji, Wi − 1, Wi ) =
'
Á N*(i− 1, j ) +N*(i, j) 1 a2a3 N*(i− 1, j )+N*(i, j)  m m m m − à à − a2 ji − 1 +ji 2 a*b*D * 2 i m,i à à j) a3 N*(i− 1, j )+ N*(i, j )à m +1, j ) +N*(i, m m m à − a2a3D*m,i − N*(i − Ä a*b* ji − 1 +ji 2a*b* 2 Å i i
'
(41)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1986
and 1 c R(ji, ji + 1, Wi, Wi + 1) = [I− Q(Wi ) Q − 1(Wi )]G. (ji, ji + 1, Wi, Wi + 1)= 2
'
Á N*(i, j )+ N*(i 1 a2a3 N*(i, j )+ N*(i+ 1, j )  m m +1, j ) m m + à − a2 à ji + ji + 1 2 a*b*D 2 * i m,i à à j ) + N*(i j )+ N*(i+ 1, j ) à a3 N*(i, m m +1, j ) m m à a2a3D*m,i N*(i, − a*b* ji +ji + 1 2a*b* 2 Å Ä i i
'
(42)
where G. (x, y, U, V)=G((x +y)/2, (U + V)/2) The above development leads to the upwind schemes:
!
' n '
dC*(i, j) 1 N*(i 1 m m +1, j ) − N*(i m −1, j ) − =− a2 dts 1 Dj 2 2 + a2
a2a3D*m,i [C*(i+ 1, j )− 2C*(i, j )+C*(i− 1, j )] m m m a*b* i
N*(i− 1, j )+N*(i, j ) N*(i, j ) +N*(i 1 m m m m +1, j ) + − ji − 1 + ji ji +ji + 1 2
(i = 2, 3, …, N; j= 1, 2, …, m2 +2)
!
j) 1 a3D*m,i C*(i 1 dN*(i, m m +1, j ) − C*(i m −1, j ) =− − dts 1 Dj a*b* 2 2 i +
'
'
"
n
a2a3 N*(i+ 1, j )− N*(i− 1, j ) m m 2 a*b*D *m,i i
n
a2a3D*m,i [N*(i+ 1, j )− 2N*(i, j )+N*(i− 1, j )] m m m a*b* i
n
(43)
"
a3 N*(i− a2a3D*m,i N*(i−1, j ) +N*(i, j ) N*(i, j ) +N*(i+ 1, j ) 1, j )+ 2N*(i, j )+ N*(i− 1, j ) m m m m m m − m + 2a*b* a*b* ji − 1 +ji ji +ji + 1 2 i i
(i = 2, 3, …, N; j= 1, 2, …, m2 +2)
(44)
4. Fitted mesh collocation method for parabolic equation in macroparticles The differential equation in macroparticles involves a parabolic equation. Often Crank–Nicolson or other finite-difference and finite-element methods have been used to solve the problem (Sohn & Szekely, 1972). Useful though such techniques are, their application to practical problems can present difficulties, particularly when steep gradients exist in the grains. Indeed when g 1 or g 1,the differential equations involve singularly perturbed problems. The stability of classical finite-difference and finite-element methods for the time-dependent differential equations then depends on the associated parameter g (Liu & Zheng, 1986; Liu & McElwain, 1997; Liu, McElwain & Donskoi, 1998). A combination of a collocation technique and a suitable integration method has been traditionally used for providing the approximate solution. However, the numerical integration of these equations is only effective when g is of the order of unity. When the parameter g is either very large or very small, the numerical integration can often present difficulty, and in order to overcome this problem a fitted mesh collocation method has been presented by Liu and Bhatia (1999). Using a combination of a fitted mesh technique and the collocation method, the PDE system in macroparticles is reduced to an ordinary differential equation. In this section, this technique is used. In the macroparticles, the unit interval lying on the h-axis is divided into the two subintervals VM1 = [0, h¯ *] c and VM2 = [h¯ *, ¯ *c =1 −h¯ c =1 −min{1/2, a1 ln(m + 1)} c 1], where h In the subinterval VM1, we define h1 =h/h¯ *, c and Eqs. (13) and (15) then become
n
a ( (C*M1 (C*M1 = 1 h2 +a4[N*(t m s 1, j, h1)]j = 1, (ts 1 h¯ *c 2h 21 (h1 1 (h1 h1 = 0:
(C*M1 =0. (h1
(45) (46)
At boundary h = h¯ *, c continuity of C * M and its derivative must be imposed: C*M1 h 1 = 1 =C*M2 h 2 = 0,
)
)
1 (C*M2 1 (C*M1 = . (h2 h 2 = 0 h¯ *c (h1 h 1 = 1 (1− h¯ *) c
(47) (48)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1987
In subinterval VM2 we define h2 =(h−h¯ *)/(1 − h¯ *) and Eqs. (13) and (14) become c c
n
(C*M2 2 (C*M2 1 ( 2C*M2 =a2 + + a4[N*(t m s 1, j, h2)]j = 1, 2 2 [h2(1 − h¯ *) ¯ *](1 − h¯ *) (ts 1 (1− h¯ *) (h 2 (h2 c c +h c c h2 = 1:
(49)
1 (C*M2 =Bi(C*e −C*M2 ) (1−h¯ *) (h2 c
(50)
Eqs. (45) and (49) are transformed into a system of m2 ordinary differential equations of the initial-value type by the orthogonal collocation technique of Villadsen and Stewart (1967). Here m2 = m1 + m1, m1 = (m− 1)/2. Roots of Jacobi polynomials are used to space the collocations points in the two subintervals V*M1 = [0, 1] and V*M2 =[0, 1]. We note that the problems are symmetric about h1 = 0, but the problems have no special symmetry properties about h2 =0. Orthogonal collocation can be applied directly to Eqs. (45) and (49) to obtain the working equations: dC*M1 ( j ) a1 m1 + 1 = 2 % B (1) C* (l) +a4N*(N + 1, j ), m h¯ *c l = 1 j,i M1 dts 1
( j=1, …, m1)
(51)
n
m1 + 2 m1 + 2 dC*M2 ( j ) 1 2 = a1 % B (2) % A (2) 1, j +m1), j,i C* M2 (l) + j,i C* M2 (l) + a4N*(N+ m 2 (2) dts 1 (1− h¯ *) [x j (1 − h¯ *)+h ¯ *](1− h¯ *) l=1 l=1 c c c c
( j = 2, …, m1 +1)
(52)
where the matrices A (1) and B (1) correspond to spherical geometry with symmetry properties, while the points x (2) and the matrices A (2) and B (2) correspond to slab geometry with no special symmetry properties (Villadsen & Stewart, 1967; Finlayson, 1972). The boundary conditions, (Eqs. (47) and (48)), are also transformed, resulting in the following forms: C*M1 (m1 + 1)= C*M2 (1) =
n
m1 + 2 1 1 1 m1 (1) (2) % A C * (l)− % A C* (l) . 1,l M (2) 2 ((A (1) ¯ *) ¯ *)) h¯ *c l = 1 m 1 + 1,l M1 (1− h¯ *) l=2 m 1 + 1,m 1 + 1/h c −A 1,1/(1− h c c
(53)
5. Numerical treatment In order to solve the ordinary differential equations (Eqs. (40), (41), (48) and (49)), the boundary conditions are numerically treated. At boundary j = 0. Eq. (19) for C *m can be approximated by a second order finite difference scheme, and we obtain 1 C*(i, j )= [4C*(2, j )− C*(3, j )], m m m 3
( j =1, …, m2 + 1)
(54)
and N*(1, j )= 0, m
( j =1, …, m2 +1)
(55)
At boundary j=1, Linear isotherm: C*(N + 1, j ) =C*M ( j ), m N*(N+ 1, j )= m
n
1 2Dj dC*M ( j ) 4N*(N, j ) − N*(N − 1, j )− , m m dts 1 4Dj +3 a2
Langmuir isotherm: C*(N + 1, j ) = m N*(N+ 1, j )= m
( j =1, …, m2 + 1)
C*M ( j )*(1 + K*) m , 1 + K*C m * M( j )
(56) ( j =1, …, m2 + 1)
( j= 1, …, m2 + 1)
(57) (58)
n
1 2Dj(1+K*) dC*M ( j ) m 4N*(N, j ) − N*(N − 1, j )− , m m 2 4Dj +3 a2[1+ K*C * ( j )] dts 1 m M
( j =1, …, m2 + 1)
(59)
where C*M ( j )=C*M1 ( j ), when j 5m1 +1 and C*M ( j )= C*M2 ( j −m1 − 1), when m1 + 1Bj. At boundary h2 =1, orthogonal collocation can be applied directly to Eq. (14) to obtain the following equation: 1 m1 + 2 (2) C* (l) = Bi(C*e − C*M2 (m1 +2)), % A 1− h¯ *c l = 1 m 1 + 2,l M2
(60)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1988
In this paper, this new numerical method (USFMCMD) uses a combination of upwind schemes, fitted mesh technique, collocation method and the DASSL integration method. DASSL uses the backward differentiation formulas of order one through five to solve a system of the differential/algebraic equations (Petzold, 1982). Once the values of C *(i, j ) and N *(i, j ) at time C*M1 ( j ) and C*M2 ( j ) at time ts 1 have been computed, the integrals m m in Eq. (21) for the fractional uptake are evaluated by the following expression:
!
M( t =3 s
m1 + 1
m1 + 2
l=1
l=1
n
2 (2) % w (1) ¯ *c 3C*M1 (l) + % w (2) ¯ *)+h ¯ *) ¯ *)C *M2 (l) l h l (x l (1 − h c c (1−h c
Fig. 1. Variation of fractional adsorption with time under macropore control, showing the effect of N with m =7, g = 10 + 4, s= 10 − 4 and b*=10 − 3.
Fig. 2. Variation of fractional adsorption with time under macropore control, showing the effect of m with N =16, g = 10 + 4, s=10 − 4 and b*=10 − 3.
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1989
Fig. 3. Convergence of hyperbolic model in microparticle control regime, showing the effect of N on the dynamics with m =7, g = 10 − 5, s =10 − 5 and b* = 10 − 3.
Fig. 4. Convergence of macropore transport equations in microparticle control regime, showing the effect of m on the dynamics with N =16, g=10 − 5, s = 10 − 5 and b*= 10 − 3.
+ 3(1−s)
m1 + 1
m1 + 2
l=1
l=1 2 m
n"
2 (2) % w (1) ¯ *c 3S(l) + % w (2) ¯ *)+h ¯ *) ¯ *)S(l +m1) l h l (x l (1 − h c c (1−h c
(61)
in which S(l) approximates 10 j C*(j, hl )dj by Simpson’s rule, and w (1) and w (2) are the integration weights for the l l collocation scheme. 6. Numerical results and discussion In this section, we present some numerical results for adsorption problems in bidisperse solids using the above method.
1990
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
6.1. Case 1: Comparison with analytical solution for b* = 0 (Ruckenstein et al., 1971). We assume 1. linear isotherm, i.e. I *(C *M)=C *M; m 2. external mass transfer offers negligible resistance, i.e. Bi ; 3. diffusivities in microparticle are constant, i.e. a*= 1, D *= 1. m
6.1.1. Macropore diffusion control When g1, i.e. when sorbent particle is large and microparticle is small, and the particle has high capacity, micropore diffusion is fast while the macropore transport is relatively slow. The system dynamics are then controlled by macropore diffusion. Figs. 1 and 2 show the effects of the discretization parameters N and m, respectively, for the USFMCMD, when g =10 + 4, s =10 − 4 and b*= 10 − 3. Also shown in the figures is the analytical solution corresponding to b* =0 for the same values of g and s. In this simulation, type A parameters were used. In this case, the parameters a1 and a2 take on the value of 0.3334× 10 − 4 and 0.3334×10 − 4, respectively, and the equations for the macroparticle as well as the microparticle lead to singularly perturbed problems. For this case the microparticle dynamics are not of significance and convergence to the analytical solution for b*= 0 is expected. It is apparent from Figs. 1 and 2 that the numerical solutions (USFMCMD) yields convergence, and the result is independent of N for N ] 4, and of m for m ]5. The USFMCMD took only 0.55 s of CPU time on a Pentium II/266 computer and is highly efficient. 6.1.2. Micropore diffusion control When g 1, i.e. when the sorbent particle is small, microparticle is large, and the particle has low capacity, micropore diffusion is slow, while the macropore is relatively fast. Thus, the system dynamics are controlled by micropore diffusion. Figs. 3 and 4 show the effects of N and m respectively for the USFMCMD, with g=10 − 5, s = 10 − 5 and b*=10 − 3. Also shown in the figures is the analytical solution corresponding to b*= 0 for the same values of g and s. In this simulation, type B parameters are used, and a1 = 1 and a2 = 10 − 10, so the microparticle equation yields a singularly perturbed problem. Due to the small value of gb* the local relaxation time is insignificant and convergence to the analytical solution for b* =0 is anticipated. It is apparent from Figs. 3 and 4 that numerical solutions (USFMCMD) are in good agreement with the analytical solution. It can be seen that the USFMCMD yields convergence and is independent of N for N] 8, and of m for m] 3. Again it was seen that method is highly efficient (it took about 0.55 s of CPU time).
Fig. 5. Variation of fractional adsorption with time under microparticle transport control, showing the effect of N with m =7, g= 10 − 5, s=10 − 5, b* =10 + 3 and Langmuir isotherm.
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1991
Fig. 6. Variation of fractional adsorption with time under microparticle transport control, showing the effect of b* with N = 30, m =7, g=10 − 5, s=10 − 5 and a Langmuir isotherm.
Fig. 7. Simulated micropore concentration profiles for a bidisperse sorbent particle under microparticle transport control, with N =30, m=7, g =10 − 5, s = 10 − 5, b* = 10 + 4 and a Langmuir isotherm.
6.2. Case 2: Simulation of the adsorbate transport with a Langmuir isotherm. We assume that 1. Langmuir isotherm, i.e. I*(C *M )= m
C*M (1+ K*) m ; 1+ K*C m * M
D*(C *) m m =
1 ; 1 + K*(1 − C*) m m
a*(C*)= 1+ K*(1− C*); m m m
2. external mass transfer offers negligible resistance, i.e. Bi .
1992
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
6.2.1. Micropore diffusion control We consider a single particle sorption model based on micropore diffusion control g = 10 − 5, s =10 − 5 with the Langmuir isotherm (K *m =1). Under this condition, the macropore dynamics is not of significance, and convergence with respect to the microparticle discretization parameter is of interest. Fig. 5 shows the effect of N, with m=7 and b*=103 on the USFMCMD fraction adsorption curve, demonstrating convergence for N] 20. It may be noted that in the microparticle control regime (g 1) the relative significance of the relaxation time is indicated by the parameter b*g, as the time can now be rescaled to gt. The convergence of the hyperbolic model with respect to N, even for the low value of b*g of 10 − 2, indicates the robustness of the present scheme. Fig. 6 shows the effect of b*. As indicated above the parameter b*g reflects the significance of the local non-equilibrium in this regime and, as seen in Fig. 6,
Fig. 8. Variation of fractional adsorption with time under macroparticle transparent control, showing the effect of b* with N =30, m =7, g= 10 + 4, s = 10 − 4 and a Langmuir isotherm.
Fig. 9. Variation of fractional adsorption with time, showing the effect of b* with N =30, m= 7, g= 1, s= 10 − 4 and a Langmuir isotherm.
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1993
Fig. 10. Variation of fractional adsorption with time, showing the effect of b* with N =30, m = 7, g=1, s= 10 − 4 and a Langmuir isotherm, during desorption.
with decrease in the value of this quantity below about 10 − 3 convergence to the classical parabolic model (with b*g= 0) occurs. The latter solution was obtained using a collocation scheme (Liu & Bhatia, 1999) in place of the upwind method described here. Even at b*g =10 − 2 the effect of local non-equilibrium is weak, but is very significant for b*g =0.1. Interestingly, with increase in b*g the system is initially sluggish, but has more rapid uptake at larger times compared with the case for b*g =0. Fig. 7 shows the simulated micropore concentration profiles for a bidisperse sorbent particle at various times, obtained using USFMCMD N= 30, m=7, b*= 104, indicating the moving steep profiles under the chosen condition. These steep profiles are akin to a moving step function, consistent with the wave nature of hyperbolic systems.
6.2.2. Macropore diffusion control We consider a single particle sorption model based on macropore diffusion control g = 10 + 4, s = 10 − 4 with the Langmuir isotherm (K *m =1). Fig. 8 shows the effect of b*, with N= 30 and m= 7, on the USFMCMD fraction adsorption curve. In this case since macropore transport controls, the micropore dynamics is not expected to influence the uptake significantly. This is evident from the figure, demonstrating insensitivity to the value of b*, even up to b*= 1. Convergence of the results for the chosen value of m was confirmed by varying this parameter. 6.2.3. Bimodal diffusion control When g= O(0), both macropore diffusion and micropore diffusion are important. The system dynamics therefore depends on the macrosphere as well as microsphere sizes. Figs. 9 and 10 show the effect of b* with N= 30, m=7, and g= 1 and s= 10 − 4 using USFMCMD for adsorption and desorption, respectively. In this case the effect or b* is not as strong as that under microparticle diffusion control; however, as before with increase in b* the initial response is more sluggish, but later on it overtakes the result for small values of b*. This occurs for both adsorption and desorption, and is due to the wave like nature of the transport in the presence of finite local relaxation time.
7. Conclusions Numerical techniques for adsorption problems involving non-linear hyperbolic equations and steep gradients in bidisperse solids have been described in the present article. Upwind schemes were applied to hyperbolic equations with source terms in the microparticle, and fitted mesh collocation technique was applied to the parabolic equation in the macroparticle. A combination of the upwind schemes, fitted mesh collocation method and the differential-al-
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1994
gebraic equation solver DASSL is applied to solve the problem. This method is numerically evaluated for comparison with the analytical solution (b* =0) when the adsorption isotherm is linear, and for convergence of the upwind scheme in the microparticle as well as the fitted mesh collocation in the macroparticle. The method is rapid and highly efficient. This method is also used to examine adsorption and desorption problems when the adsoption isotherm has the Langmuir form. Studies on the effect of the model parameter representing the local relaxation time show this to have significant influence on the dynamics particularly under conditions of microparticle transport control.
8. Notations Bi CM Cm Cmb Cb C0 Cm0 C *M C *m De Dp Dm Dm0 D *m FA Im I *m km Km K *m Nm N *m rM rm Rp Rm T Greek a(Cm) a* ab s1 b* oM g s j h t ts ts 1
Biot number (=kmRp/De) macropore concentration microsphere concentration (=Im(Cb)) bulk concentration during adsorption, and C0 during desorption initial (t= 0) macropore concentration microsphere concentration in equilibrium with Cb (=CM/Cb) for adsorption, (= CM/C0) for desorption (=Cm/Cmb) effective diffusivity particle macropore diffusivity (=De/oM) micropore diffusivity [= Dm(Cb)] dimensionless surface diffusion coefficient fractional adsorption a suitable micropore adsorption isotherm dimensionless isotherm external mass transfer coefficient Langmuirian constant (=Km×Cb) for adsorption, (= Km×C0) for desorption micropore flux (NmRm/[CmbDm(Cb)]) particle radical coordinate microsphere radial coordinate radius of particle radius of microsphere time coordinate letters local relaxation time (= a/ab) [= a(Cmb)] [= 3(1−s)g)] [= Dpsab/R 2p] macroporosity [= Dm0R 2p/(DpR 2ms)] {=oMCb/[oMCb+(1−oM)Cmb]} (= rm/Rm) (=rM/Rp) (= sDpt/R 2p) (=t/s) (=s1ts )
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1995
Acknowledgements This research has been supported by the Australian Research Council through the Large Grant Scheme. References Abarzhi, I. I. (2000). Wave mechanism of mass transfer for kinetics of adsorption in biporous media. Colloids and Surfaces A, 164, 105–113. Bermudez, A., & Vazquez, M. E. (1994). Upwind methods for hyperbolic conservation laws with source terms. Computational Fluid, 23, 1049 – 1071. Bhatia, S. K. (1987). Modelling the pore structure of coal. American Institute of Chemical Engineering Journal, 33, 1707 – 1718. Bhatia, S. K. (1997). Transport in bidisperse adsorbents: significance of the macroscopic adsorbate flux. Chemical Engineering Science, 52, 1377 – 1386. Bhatia, S. K., Gray, P. G., & Do, D. D. (1991). Modelling of sorption of gaseous sorbates in bidispersed structured solids. Gas Separation and Purification, 5, 49 – 55. Finlayson, B. A. (1972). The method of weighted residuals and 6ariational principles. Academic Press. Gray, P. G., & Do, D. D. (1990). Adsorption and desorption of gaseous sorbates on a bidispersed particle with Freundlich isotherm-III contribution of surface diffusion to the sorption dynamics of sulphur dioxide on activated carbon. Gas Separation and Purification, 4, 149–157. Godlewski, E., & Raviart, P. A. (1991). Hyperbolic systems of conservation laws, Mathematiques & Applications, N. 3/4, Ellioses-Edition Marketing. LeVeque, R. (1990). Numerical methods for conser6ation laws. Basel: Birkha¨user. Liu, F., & Bhatia, S. K. (1999). Computationally efficient solution techniques for adsorption problems involving steep gradients in bidisperse particles. Computers and Chemical Engineering, 23, 933–943. Liu, F., & McElwain, D. L. S. (1997). A computationally efficient solution technique for moving-boundary problems in finite media. IMA Journal of Applied Mathematics, 59, 71–84. Liu, F., & Zheng, X. (1986). Uniformly convergent difference scheme for the singular perturbation of a self adjoint elliptic partial differential equation. Applied Mathematical Mechanics, 7, 495–504. Liu, F., McElwain, D. L. S., & Donskoi, E. (1998). The use of a modified Petrov – Galerkin method for gas – solid reaction modelling. IMA Journal of Applied Mathematics, 60, 1–14. Luikov, A. V. (1980). Heat and mass transfer. Moscow: Mir Publishers. Petzold, L. R. (1982). A Description of DASSL: A Differential/Algebraic Equation System Solver, Sandia Technical Report SAND82-8637 Livermore, CA USA. Ruckenstein, E., Vaidyanathan, A. S., & Youngquist, G. R. (1971). Sorption by solids with bidisperse pore structures. Chemical Engineering Science, 26, 1305 – 1318. Ruthven, D. M. (1984). Principles of adsorption and adsorption processes. New York, USA: Wiley. Sohn, H. Y., & Szekely, J. (1972). A structural model for gas– solid reactions with a moving boundary-III. A general dimensionless representation of the irreversible reaction between a porous solid and a reactant gas. Chemical Engineering Science, 27, 763 – 778. Villadsen, J. V., & Stewart, W. E. (1967). Solution of boundary-value problems by orthogonal collocation. Chemical Engineering Science, 22, 1483 – 1501.
.
.
Numerical solution of hyperbolic models of transport in bidisperse solids F. Liu a, S.K. Bhatia a,*, I.I. Abarzhi b b
a Department of Chemical Engineering, The Uni6ersity of Queensland, Brisbane, Qld 4072, Australia Engineering Thermophysics Institute, Ukranian Academy of Sciences, Zhelyabo6a 2a, 252057 Kie6, Ukraine
Received 17 August 1999; received in revised form 17 April 2000; accepted 22 May 2000
Abstract This article considers the numerical solution of generalized transport models for adsorption that include the presence of finite local relaxation times. Such models arise in adsorption problems where the small diffusivity in micropores can lead to local nonequilibrium even at the macroscopic process time scale. We consider the important case of a bidisperse-structured adsorbent, the transport in which is a problem of long-standing interest to chemical engineers. In the formulation considered, the microparticle model involves a non-linear hyperbolic differential equation with source terms, and in the macroparticle model the differential equation involves a singularly perturbed parabolic diffusion problem. There are some difficulties in solving these equations due to the hyperbolic nature of the microparticle transport, conventionally assumed to obey the Fickian model leading to a parabolic diffusion equation. In this paper, an efficient numerical technique is presented to simulate these processes. A combination of the upwind method and fitted mesh collocation method is applied to solve the problem. Numerical solutions were found to match the analytical solution of the traditional model when the adsorption isotherm is linear, for macropore diffusion control or micropore diffusion control. Simulations of the adsorption and desorption dynamics are also presented for a Langmuir isotherm. The numerical scheme offers a more generalized alternative that can be used for both the model forms, with and without consideration of finite local relaxation times. © 2000 Elsevier Science Ltd. All rights reserved. Keywords: Adsorption; Bidispersed particles; Upwind method; Fitted mesh; Collection method; Hyperbolic
1. Introduction The description of adsorbate transport in microporous solids is an important component of the design and optimization of the adsorption processes. Generally microporous adsorbents, such as activated carbon or pelletized zeolites, are bidisperse and have both micropores as well as macropores. The description of transport processes in such bidisperse solids is difficult due to the complex and largely unknown nature of the pore network, and a variety of structural idealizations have been proposed, based on which approximate methods for diffusion and adsorption have been developed. The simplest and the most popular of these idealizations considers the solid particle as comprising an aggregate of microporous spherical grains, with the intergrain macropores providing the channels for intraparticle transport. The modeling of transport in such a solid was first presented by Ruckenstein, Vaidyanathan and Youngquist (1971), and considered the diffusion processes over two length scales: the particle scale (macropore diffusion) and the local/microscopic scale (micropore diffusion). Since this development, numerous studies have used this model for interpreting adsorption/desorption dynamics in porous adsorbents (Bhatia, 1987); however, some workers (Gray & Do, 1990) have added further arbitrary mechanisms representing particle scale transport of
* Corresponding author. Tel.: +61-7-33654263; fax: + 61-7-33654199. E-mail address: [email protected] (S.K. Bhatia) 0098-1354/00/$ - see front matter © 2000 Elsevier Science Ltd. All rights reserved. PII: S0098-1354(00)00596-2
1982
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
adsorbed species in order to improve the fitting of the experimental data. A sounder basis for particle scale adsorbate transport has been subsequently provided (Bhatia, 1997; Bhatia, Gray & Do, 1991) by considering the asymmetric macropore concentration profile over the grain surface, and the resulting two-dimensional transport in the grain. This leads to a through flux across the grains, and facilitates particle scale adsorbate transport, though at the cost of considerable additional complexity. Most recently Abarzhi (2000) has provided a further development in the analysis of adsorbate transport by considering the effect of local relaxation times in the macropore and micropore networks. The classical model for micropore diffusion in the grains (Ruthven, 1984) involves use of the Fickian flux model Nm = − Dm(Cm)9Cm
(1)
with a concentration dependent diffusivity following the Darken relation. As indicated by Luikov (1980) and Abarzhi (2000) such a model neglects the effect of local relaxation times in the diffusion process, and is not exact. While negligible at the single pore level the local relaxation time has been considered important for the intragrain pore network as a whole in some cases. When finite local relaxation time is considered, the resulting microparticle flux model is non-Fickian, and has the form a
(Nm +Nm = − Dm(Cm)9Cm (t
(2)
in which a represents a local relaxation time scale, that may be expected to be inversely proportional to the microparticle diffusivity. In the presence of a concentration-dependent diffusivity, as is the case for adsorbate transport, a would be inversely proportional to Dm(Cm). In general it may be anticipated that the local relaxation time is negligible in the intergrain macropores, where the diffusivities are three to four orders of magnitude larger that those in the micropores. Such a model would be most relevant to molecular sieving carbons and zeolites where the extremely small micropore diffusivities may lead to finite relaxation times. Abarzhi (2000) has shown that the local non-equilibrium is important at times comparable with the local relaxation time, although near equilibrium the flux law approaches the Fickian formulation in Eq. (1). To enable the application of adsorption models involving finite local relaxation times to actual experimental data it is pertinent that the numerical solution of such models be explored, as analytical solutions are extremely complex and not possible for the case of concentration-dependent diffusivities. The solution of such models must clearly be a precursor to their application, and studies in this direction should facilitate further research along these lines and experimental validation of the models. To this end the present article is devoted to the numerical solution of the transport model for adsorption in bidisperse solids, with the intragrain flux represented by Eq. (2) in place of Eq. (1). In the present case, combination of Eq. (2) with the balance equation (Cm = − 9 · Nm (t
(3)
for the grain micropores, provides for a hyperbolic conservation law as opposed to the parabolic model arising from the use of the classical Fickian approach. This requires specialized and complex solution techniques, a problem that is compounded by the coupling of the grain balance equation with the macropore balance equation, which is parabolic in nature. In addition, the latter often involves a singularly perturbed problem (Liu & Bhatia, 1999). In the present article, we propose a combination of an upwind method (Bermudez & Vazquez, 1994) for the intragrain hyperbolic transport and the fitted mesh collocation technique (Liu & Bhatia, 1999) for the macropore transport, along with the use of the differential/algebraic system solver package (DASSL) developed by Petzold (1982). The latter considerably improves the speed of the upwind methods, which conventionally rely on finite difference approximations. Although upwind methods are not new, their use in hyperbolic mass-transport problems has not received significant attention. The present article presents an initial study in this direction, demonstrating the success of these techniques for such problems. For sufficiently small relaxation times it is demonstrated that the solution converges to that of the parabolic model as expected, so that the technique proposed is a generalization valid for both the model forms.
2. Mathematical model We consider simultaneous diffusion and adsorption in a spherical sorbent pellet comprising an aggregate of small
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1983
spherical microparticles of uniform size. The mass balance equation for the distribution of sorbate in the macrosphere is (Liu & Bhatia, 1999): s
n
(C*M 1 ( (C*M = 2 h2 +3(1 −s)g[N*(t, h, j)]j = 1 m (t h (h (h
(4)
with boundary conditions h= 1:
(C*M = Bi(C*e − C*M ) (h
(5)
h= 0:
(C*M = 0. (h
(6)
A differential mass balance equation for the distribution of sorbate in the microsphere is
n
(C*m 2N*m (N*m = −g + , (t j (j
(N*m (C*m *) = − N*m +D*(C m m (t (j
b*a*(C*) m
(7)
n
(8)
with boundary conditions j= 1: C*= I*(C *M ), m m j =0:
(9)
(C*m = 0, N*m = 0 (j
(10)
where I *m is a dimensionless adsorption isotherm. I *(C *M)= C *M when the sorption isotherm is linear, and I *(C *M) = m m C *M(1+K *)/(1 +K *C *M), when the sorption isotherm follows the Langmuir form. Here b*= Dpsab/R 2p. m m For the adsorption cycle, the initial condition is t= 0: C*M =C*=0, C*e = 1, N*m = 0 m
(11)
while for the desorption cycle, the initial condition is t =0: C*M =C*=1, C*e = 0, N*m = 0. m
(12)
The parameter s is the fraction of equilibrium uptake in macro-voids and the quantity (1− s) is that in the microsphere. Parameter g is called the dynamic parameter because it is a measure of the relative rates of diffusion in the macropores and diffusion in the micropores. We now rewrite the above non-dimensional model in the form:
n
1 ( (C*M (C*M =a1 2 h2 +a4[N*] m j = 1, (ts h (h (h 1
(13)
h=1:
(C*M = Bi(C*e − C*M ), (h
(14)
h= 0:
(C*M =0 (h
(15)
and (C*m (N*m N* +a2 = − 2a2 m, (ts 1 (j j
(16)
(N*m a3D*(C *) a3 m m (C* m + =− N*, m (ts 1 b*a*(C*) (j b*a*(C*) m m
(17)
j =1: C*= I*(C *M ), m m
(18)
j =0:
(C*m =0, N*= 0 m (j
with the model parameters having two alternate definitions: Type A: a1 =1/s1, a2 =sg/s1, a3 =s/s1, a4 =1, ts 1 = tss1.
(19)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1984
Type B: a1 =1, a2 =sg, a3 =s, a4 =s1, ts 1 =ts where s1 =3(1−s)g, ts =t/s. The fractional uptake is evaluated by determining the amount of sorbate that has entered the particle after changing the external concentration. The fractional uptake (Liu & Bhatia, 1999) is given by FA =
M( t M(
with
(20)
&
M( t =3 s
1
& &
h 2C*M dh+ 3(1 −s)
0
1
h2
0
1
n
j 2C*djdh , m
0
M( =s + (1− s)I*(1) m
(21) (22)
3. Upwind method for hyperbolic conservation laws with source terms in microparticles The differential equations in the microparticle involve nonlinear hyperbolic Eqs. (16) and (17) with source terms. The section deals with the extension of some upwind schemes to hyperbolic systems of conservation laws with source terms. Upwind methods for hyperbolic conservation laws have undergone considerable development in recent years, leading to a variety of implementation techniques such as flux-difference and flux-splitting methods (Bermudez & Vazquez, 1994). Moreover, the foundations of the theoretical work are well established. A review of these has been given by LeVeque (1990) and Godlewski and Raviart (1991). In this section, the upwind methods are only used in spatial descretization. The partial differential equations (Eqs. (16) and (17)) in the microparticle are reduced to an ordinary differential equation system. Let us consider the hyperbolic system of conservation laws with source terms (v(j, ts 1) (F(v(j, ts 1)) =G(j, v(j, ts 1)) + (ts 1 (j where v(j, ts 1)=
(23)
C*(j, ts 1) m , ts 1) N*(j, m
(24)
Á  a2N*m à à F(v)= Ãa D*(C*) Ã, à 3 m m C*m à Äa*(C*)b* Å m
(25)
Á N*  −2a2 m à à j G(j, v)= à Ã. a3 Ã− N*m Ã Ä a*(C*)b* Å m
(26)
To discretize in space we use the finite volume method. The domain is discretized by taking a uniform mesh {ji =(i− 1)Dj, Dj = ji −ji − 1, (i =1, …, N + 1)} and each interval is divided into two to form a cell Ci = (ji −Dj/2, ji + Dj/2. Let us consider a discretization in which the notes do not lie on the boundary. The approximate solution W(j, ts 1) to v(j, ts 1) at ts 1 can be considered to be a piecewise constant function given by
W(ji, ts 1)=Wi (ts 1) for j Ci = ji −
Dj Dj , ji + . 2 2
(27)
To obtain the spatial approximation of Eq. (23), we integrate this over the cell Ci. This leads to the following scheme: dW Fi + (1/2) −Fi − (1/2) + =G(ji, W(ji, ts 1)) = Gi dts 1 Dj
(28)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1985
where Fi 9 (1/2) and Gi denote approximations of F(W(ji 9 (1/2))) and (1/Dj) Ci G(j, W)dj, respectively. Upwind schemes are obtained by replacing Fi 9 (1/2) by a numerical flux function f. In a similar way we introduce a numerical source function c to upwind the source term. More precisely, we take for each i: Fi + (1/2) = fi (Wi, Wi + 1),
(29)
Fi − (1/2) = fi (Wi − 1, Wi ),
(30)
Gi = c(ji − 1, ji, ji + 1, Wi − 1, Wi, Wi + 1).
(31)
This represents an upwind scheme corresponding to numerical fluxes F(Wi − 1) +F(Wi ) 1 − Q(Wi ) (Wi − Wi − 1), 2 2 F(Wi ) +F(Wi + 1) 1 − Q(Wi ) (Wi + 1 − Wi ), fi (Wi, Wi + 1)= 2 2
fi (Wi − 1, Wi )=
(32) (33)
where Q is a matrix characteristic of each Q-scheme having a continuous dependence on Wi. Some notations in Bermudez and Vazquez (1994) will be used in this paper. For any diagonalizable matrix A, let L be the diagonal form of A given by its eigenvalues, then A = XLX − 1
(34)
where X is the matrix of the eigenvetors of A. We recall that A = X L X − 1
(35)
and L ij = li dij.
(36)
Q-scheme equals to the Jacobian matrix A of the flux valued at the value of Wi : Q(Wi ) = A(Wi ) =X L X − 1 =diag( l1(Wi ) , l2(Wi ) )
(37)
where diag ( l1(Wi ) , l2(Wi ) ) is a diagonal matrix and l1(Wi )=
'
a2a3D*m,i , a*b* i
l2(Wi ) = −
'
a2a3D*m,i a*b* i
(38)
*(j Here D*mi =D*(C m m i, tsi ))and a* i = a*(C*(j m i, ts 1)) The discretization of the flux term is Fi + (1/2) − Fi − (1/2) 1 1 = [fi (i, i +1) −fi (i −1, i )]= Dj Dj Dj
' '
Á N*(i+ 1, j )−N*(i 1 a2a3D*m,i  m m −1, j ) − [C*(i+ 1, j )−2C*(i, j )+ C*(i− 1, j )] à m m m à a2 2 2 a*b* à à a D * C *(i+ 1, j ) −C *(i −1, j ) 1 a2a3D*m,i 3 m,i m m à − [N*(i+ 1, j )− 2N*(i, j )+ N*(i− 1, j )] à m m m Ä a*b* 2 2 a*b* Å i i
(39)
The numerical source function c will be taken of the form c(ji − 1, ji, ji + 1, Wi − 1, Wi, Wi + 1) = c L(ji − 1, ji, Wi − 1, Wi ) + c R(ji, ji + 1, Wi, Wi + 1)
(40)
where left and right numerical source functions are given by 1 c L(ji − 1, ji, Wi − 1, Wi ) = [I+ Q(Wi ) Q − 1(Wi )] 2 G. (ji − 1, ji, Wi − 1, Wi ) =
'
Á N*(i− 1, j ) +N*(i, j) 1 a2a3 N*(i− 1, j )+N*(i, j)  m m m m − à à − a2 ji − 1 +ji 2 a*b*D * 2 i m,i à à j) a3 N*(i− 1, j )+ N*(i, j )à m +1, j ) +N*(i, m m m à − a2a3D*m,i − N*(i − Ä a*b* ji − 1 +ji 2a*b* 2 Å i i
'
(41)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1986
and 1 c R(ji, ji + 1, Wi, Wi + 1) = [I− Q(Wi ) Q − 1(Wi )]G. (ji, ji + 1, Wi, Wi + 1)= 2
'
Á N*(i, j )+ N*(i 1 a2a3 N*(i, j )+ N*(i+ 1, j )  m m +1, j ) m m + à − a2 à ji + ji + 1 2 a*b*D 2 * i m,i à à j ) + N*(i j )+ N*(i+ 1, j ) à a3 N*(i, m m +1, j ) m m à a2a3D*m,i N*(i, − a*b* ji +ji + 1 2a*b* 2 Å Ä i i
'
(42)
where G. (x, y, U, V)=G((x +y)/2, (U + V)/2) The above development leads to the upwind schemes:
!
' n '
dC*(i, j) 1 N*(i 1 m m +1, j ) − N*(i m −1, j ) − =− a2 dts 1 Dj 2 2 + a2
a2a3D*m,i [C*(i+ 1, j )− 2C*(i, j )+C*(i− 1, j )] m m m a*b* i
N*(i− 1, j )+N*(i, j ) N*(i, j ) +N*(i 1 m m m m +1, j ) + − ji − 1 + ji ji +ji + 1 2
(i = 2, 3, …, N; j= 1, 2, …, m2 +2)
!
j) 1 a3D*m,i C*(i 1 dN*(i, m m +1, j ) − C*(i m −1, j ) =− − dts 1 Dj a*b* 2 2 i +
'
'
"
n
a2a3 N*(i+ 1, j )− N*(i− 1, j ) m m 2 a*b*D *m,i i
n
a2a3D*m,i [N*(i+ 1, j )− 2N*(i, j )+N*(i− 1, j )] m m m a*b* i
n
(43)
"
a3 N*(i− a2a3D*m,i N*(i−1, j ) +N*(i, j ) N*(i, j ) +N*(i+ 1, j ) 1, j )+ 2N*(i, j )+ N*(i− 1, j ) m m m m m m − m + 2a*b* a*b* ji − 1 +ji ji +ji + 1 2 i i
(i = 2, 3, …, N; j= 1, 2, …, m2 +2)
(44)
4. Fitted mesh collocation method for parabolic equation in macroparticles The differential equation in macroparticles involves a parabolic equation. Often Crank–Nicolson or other finite-difference and finite-element methods have been used to solve the problem (Sohn & Szekely, 1972). Useful though such techniques are, their application to practical problems can present difficulties, particularly when steep gradients exist in the grains. Indeed when g 1 or g 1,the differential equations involve singularly perturbed problems. The stability of classical finite-difference and finite-element methods for the time-dependent differential equations then depends on the associated parameter g (Liu & Zheng, 1986; Liu & McElwain, 1997; Liu, McElwain & Donskoi, 1998). A combination of a collocation technique and a suitable integration method has been traditionally used for providing the approximate solution. However, the numerical integration of these equations is only effective when g is of the order of unity. When the parameter g is either very large or very small, the numerical integration can often present difficulty, and in order to overcome this problem a fitted mesh collocation method has been presented by Liu and Bhatia (1999). Using a combination of a fitted mesh technique and the collocation method, the PDE system in macroparticles is reduced to an ordinary differential equation. In this section, this technique is used. In the macroparticles, the unit interval lying on the h-axis is divided into the two subintervals VM1 = [0, h¯ *] c and VM2 = [h¯ *, ¯ *c =1 −h¯ c =1 −min{1/2, a1 ln(m + 1)} c 1], where h In the subinterval VM1, we define h1 =h/h¯ *, c and Eqs. (13) and (15) then become
n
a ( (C*M1 (C*M1 = 1 h2 +a4[N*(t m s 1, j, h1)]j = 1, (ts 1 h¯ *c 2h 21 (h1 1 (h1 h1 = 0:
(C*M1 =0. (h1
(45) (46)
At boundary h = h¯ *, c continuity of C * M and its derivative must be imposed: C*M1 h 1 = 1 =C*M2 h 2 = 0,
)
)
1 (C*M2 1 (C*M1 = . (h2 h 2 = 0 h¯ *c (h1 h 1 = 1 (1− h¯ *) c
(47) (48)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1987
In subinterval VM2 we define h2 =(h−h¯ *)/(1 − h¯ *) and Eqs. (13) and (14) become c c
n
(C*M2 2 (C*M2 1 ( 2C*M2 =a2 + + a4[N*(t m s 1, j, h2)]j = 1, 2 2 [h2(1 − h¯ *) ¯ *](1 − h¯ *) (ts 1 (1− h¯ *) (h 2 (h2 c c +h c c h2 = 1:
(49)
1 (C*M2 =Bi(C*e −C*M2 ) (1−h¯ *) (h2 c
(50)
Eqs. (45) and (49) are transformed into a system of m2 ordinary differential equations of the initial-value type by the orthogonal collocation technique of Villadsen and Stewart (1967). Here m2 = m1 + m1, m1 = (m− 1)/2. Roots of Jacobi polynomials are used to space the collocations points in the two subintervals V*M1 = [0, 1] and V*M2 =[0, 1]. We note that the problems are symmetric about h1 = 0, but the problems have no special symmetry properties about h2 =0. Orthogonal collocation can be applied directly to Eqs. (45) and (49) to obtain the working equations: dC*M1 ( j ) a1 m1 + 1 = 2 % B (1) C* (l) +a4N*(N + 1, j ), m h¯ *c l = 1 j,i M1 dts 1
( j=1, …, m1)
(51)
n
m1 + 2 m1 + 2 dC*M2 ( j ) 1 2 = a1 % B (2) % A (2) 1, j +m1), j,i C* M2 (l) + j,i C* M2 (l) + a4N*(N+ m 2 (2) dts 1 (1− h¯ *) [x j (1 − h¯ *)+h ¯ *](1− h¯ *) l=1 l=1 c c c c
( j = 2, …, m1 +1)
(52)
where the matrices A (1) and B (1) correspond to spherical geometry with symmetry properties, while the points x (2) and the matrices A (2) and B (2) correspond to slab geometry with no special symmetry properties (Villadsen & Stewart, 1967; Finlayson, 1972). The boundary conditions, (Eqs. (47) and (48)), are also transformed, resulting in the following forms: C*M1 (m1 + 1)= C*M2 (1) =
n
m1 + 2 1 1 1 m1 (1) (2) % A C * (l)− % A C* (l) . 1,l M (2) 2 ((A (1) ¯ *) ¯ *)) h¯ *c l = 1 m 1 + 1,l M1 (1− h¯ *) l=2 m 1 + 1,m 1 + 1/h c −A 1,1/(1− h c c
(53)
5. Numerical treatment In order to solve the ordinary differential equations (Eqs. (40), (41), (48) and (49)), the boundary conditions are numerically treated. At boundary j = 0. Eq. (19) for C *m can be approximated by a second order finite difference scheme, and we obtain 1 C*(i, j )= [4C*(2, j )− C*(3, j )], m m m 3
( j =1, …, m2 + 1)
(54)
and N*(1, j )= 0, m
( j =1, …, m2 +1)
(55)
At boundary j=1, Linear isotherm: C*(N + 1, j ) =C*M ( j ), m N*(N+ 1, j )= m
n
1 2Dj dC*M ( j ) 4N*(N, j ) − N*(N − 1, j )− , m m dts 1 4Dj +3 a2
Langmuir isotherm: C*(N + 1, j ) = m N*(N+ 1, j )= m
( j =1, …, m2 + 1)
C*M ( j )*(1 + K*) m , 1 + K*C m * M( j )
(56) ( j =1, …, m2 + 1)
( j= 1, …, m2 + 1)
(57) (58)
n
1 2Dj(1+K*) dC*M ( j ) m 4N*(N, j ) − N*(N − 1, j )− , m m 2 4Dj +3 a2[1+ K*C * ( j )] dts 1 m M
( j =1, …, m2 + 1)
(59)
where C*M ( j )=C*M1 ( j ), when j 5m1 +1 and C*M ( j )= C*M2 ( j −m1 − 1), when m1 + 1Bj. At boundary h2 =1, orthogonal collocation can be applied directly to Eq. (14) to obtain the following equation: 1 m1 + 2 (2) C* (l) = Bi(C*e − C*M2 (m1 +2)), % A 1− h¯ *c l = 1 m 1 + 2,l M2
(60)
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1988
In this paper, this new numerical method (USFMCMD) uses a combination of upwind schemes, fitted mesh technique, collocation method and the DASSL integration method. DASSL uses the backward differentiation formulas of order one through five to solve a system of the differential/algebraic equations (Petzold, 1982). Once the values of C *(i, j ) and N *(i, j ) at time C*M1 ( j ) and C*M2 ( j ) at time ts 1 have been computed, the integrals m m in Eq. (21) for the fractional uptake are evaluated by the following expression:
!
M( t =3 s
m1 + 1
m1 + 2
l=1
l=1
n
2 (2) % w (1) ¯ *c 3C*M1 (l) + % w (2) ¯ *)+h ¯ *) ¯ *)C *M2 (l) l h l (x l (1 − h c c (1−h c
Fig. 1. Variation of fractional adsorption with time under macropore control, showing the effect of N with m =7, g = 10 + 4, s= 10 − 4 and b*=10 − 3.
Fig. 2. Variation of fractional adsorption with time under macropore control, showing the effect of m with N =16, g = 10 + 4, s=10 − 4 and b*=10 − 3.
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1989
Fig. 3. Convergence of hyperbolic model in microparticle control regime, showing the effect of N on the dynamics with m =7, g = 10 − 5, s =10 − 5 and b* = 10 − 3.
Fig. 4. Convergence of macropore transport equations in microparticle control regime, showing the effect of m on the dynamics with N =16, g=10 − 5, s = 10 − 5 and b*= 10 − 3.
+ 3(1−s)
m1 + 1
m1 + 2
l=1
l=1 2 m
n"
2 (2) % w (1) ¯ *c 3S(l) + % w (2) ¯ *)+h ¯ *) ¯ *)S(l +m1) l h l (x l (1 − h c c (1−h c
(61)
in which S(l) approximates 10 j C*(j, hl )dj by Simpson’s rule, and w (1) and w (2) are the integration weights for the l l collocation scheme. 6. Numerical results and discussion In this section, we present some numerical results for adsorption problems in bidisperse solids using the above method.
1990
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
6.1. Case 1: Comparison with analytical solution for b* = 0 (Ruckenstein et al., 1971). We assume 1. linear isotherm, i.e. I *(C *M)=C *M; m 2. external mass transfer offers negligible resistance, i.e. Bi ; 3. diffusivities in microparticle are constant, i.e. a*= 1, D *= 1. m
6.1.1. Macropore diffusion control When g1, i.e. when sorbent particle is large and microparticle is small, and the particle has high capacity, micropore diffusion is fast while the macropore transport is relatively slow. The system dynamics are then controlled by macropore diffusion. Figs. 1 and 2 show the effects of the discretization parameters N and m, respectively, for the USFMCMD, when g =10 + 4, s =10 − 4 and b*= 10 − 3. Also shown in the figures is the analytical solution corresponding to b* =0 for the same values of g and s. In this simulation, type A parameters were used. In this case, the parameters a1 and a2 take on the value of 0.3334× 10 − 4 and 0.3334×10 − 4, respectively, and the equations for the macroparticle as well as the microparticle lead to singularly perturbed problems. For this case the microparticle dynamics are not of significance and convergence to the analytical solution for b*= 0 is expected. It is apparent from Figs. 1 and 2 that the numerical solutions (USFMCMD) yields convergence, and the result is independent of N for N ] 4, and of m for m ]5. The USFMCMD took only 0.55 s of CPU time on a Pentium II/266 computer and is highly efficient. 6.1.2. Micropore diffusion control When g 1, i.e. when the sorbent particle is small, microparticle is large, and the particle has low capacity, micropore diffusion is slow, while the macropore is relatively fast. Thus, the system dynamics are controlled by micropore diffusion. Figs. 3 and 4 show the effects of N and m respectively for the USFMCMD, with g=10 − 5, s = 10 − 5 and b*=10 − 3. Also shown in the figures is the analytical solution corresponding to b*= 0 for the same values of g and s. In this simulation, type B parameters are used, and a1 = 1 and a2 = 10 − 10, so the microparticle equation yields a singularly perturbed problem. Due to the small value of gb* the local relaxation time is insignificant and convergence to the analytical solution for b* =0 is anticipated. It is apparent from Figs. 3 and 4 that numerical solutions (USFMCMD) are in good agreement with the analytical solution. It can be seen that the USFMCMD yields convergence and is independent of N for N] 8, and of m for m] 3. Again it was seen that method is highly efficient (it took about 0.55 s of CPU time).
Fig. 5. Variation of fractional adsorption with time under microparticle transport control, showing the effect of N with m =7, g= 10 − 5, s=10 − 5, b* =10 + 3 and Langmuir isotherm.
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1991
Fig. 6. Variation of fractional adsorption with time under microparticle transport control, showing the effect of b* with N = 30, m =7, g=10 − 5, s=10 − 5 and a Langmuir isotherm.
Fig. 7. Simulated micropore concentration profiles for a bidisperse sorbent particle under microparticle transport control, with N =30, m=7, g =10 − 5, s = 10 − 5, b* = 10 + 4 and a Langmuir isotherm.
6.2. Case 2: Simulation of the adsorbate transport with a Langmuir isotherm. We assume that 1. Langmuir isotherm, i.e. I*(C *M )= m
C*M (1+ K*) m ; 1+ K*C m * M
D*(C *) m m =
1 ; 1 + K*(1 − C*) m m
a*(C*)= 1+ K*(1− C*); m m m
2. external mass transfer offers negligible resistance, i.e. Bi .
1992
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
6.2.1. Micropore diffusion control We consider a single particle sorption model based on micropore diffusion control g = 10 − 5, s =10 − 5 with the Langmuir isotherm (K *m =1). Under this condition, the macropore dynamics is not of significance, and convergence with respect to the microparticle discretization parameter is of interest. Fig. 5 shows the effect of N, with m=7 and b*=103 on the USFMCMD fraction adsorption curve, demonstrating convergence for N] 20. It may be noted that in the microparticle control regime (g 1) the relative significance of the relaxation time is indicated by the parameter b*g, as the time can now be rescaled to gt. The convergence of the hyperbolic model with respect to N, even for the low value of b*g of 10 − 2, indicates the robustness of the present scheme. Fig. 6 shows the effect of b*. As indicated above the parameter b*g reflects the significance of the local non-equilibrium in this regime and, as seen in Fig. 6,
Fig. 8. Variation of fractional adsorption with time under macroparticle transparent control, showing the effect of b* with N =30, m =7, g= 10 + 4, s = 10 − 4 and a Langmuir isotherm.
Fig. 9. Variation of fractional adsorption with time, showing the effect of b* with N =30, m= 7, g= 1, s= 10 − 4 and a Langmuir isotherm.
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1993
Fig. 10. Variation of fractional adsorption with time, showing the effect of b* with N =30, m = 7, g=1, s= 10 − 4 and a Langmuir isotherm, during desorption.
with decrease in the value of this quantity below about 10 − 3 convergence to the classical parabolic model (with b*g= 0) occurs. The latter solution was obtained using a collocation scheme (Liu & Bhatia, 1999) in place of the upwind method described here. Even at b*g =10 − 2 the effect of local non-equilibrium is weak, but is very significant for b*g =0.1. Interestingly, with increase in b*g the system is initially sluggish, but has more rapid uptake at larger times compared with the case for b*g =0. Fig. 7 shows the simulated micropore concentration profiles for a bidisperse sorbent particle at various times, obtained using USFMCMD N= 30, m=7, b*= 104, indicating the moving steep profiles under the chosen condition. These steep profiles are akin to a moving step function, consistent with the wave nature of hyperbolic systems.
6.2.2. Macropore diffusion control We consider a single particle sorption model based on macropore diffusion control g = 10 + 4, s = 10 − 4 with the Langmuir isotherm (K *m =1). Fig. 8 shows the effect of b*, with N= 30 and m= 7, on the USFMCMD fraction adsorption curve. In this case since macropore transport controls, the micropore dynamics is not expected to influence the uptake significantly. This is evident from the figure, demonstrating insensitivity to the value of b*, even up to b*= 1. Convergence of the results for the chosen value of m was confirmed by varying this parameter. 6.2.3. Bimodal diffusion control When g= O(0), both macropore diffusion and micropore diffusion are important. The system dynamics therefore depends on the macrosphere as well as microsphere sizes. Figs. 9 and 10 show the effect of b* with N= 30, m=7, and g= 1 and s= 10 − 4 using USFMCMD for adsorption and desorption, respectively. In this case the effect or b* is not as strong as that under microparticle diffusion control; however, as before with increase in b* the initial response is more sluggish, but later on it overtakes the result for small values of b*. This occurs for both adsorption and desorption, and is due to the wave like nature of the transport in the presence of finite local relaxation time.
7. Conclusions Numerical techniques for adsorption problems involving non-linear hyperbolic equations and steep gradients in bidisperse solids have been described in the present article. Upwind schemes were applied to hyperbolic equations with source terms in the microparticle, and fitted mesh collocation technique was applied to the parabolic equation in the macroparticle. A combination of the upwind schemes, fitted mesh collocation method and the differential-al-
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1994
gebraic equation solver DASSL is applied to solve the problem. This method is numerically evaluated for comparison with the analytical solution (b* =0) when the adsorption isotherm is linear, and for convergence of the upwind scheme in the microparticle as well as the fitted mesh collocation in the macroparticle. The method is rapid and highly efficient. This method is also used to examine adsorption and desorption problems when the adsoption isotherm has the Langmuir form. Studies on the effect of the model parameter representing the local relaxation time show this to have significant influence on the dynamics particularly under conditions of microparticle transport control.
8. Notations Bi CM Cm Cmb Cb C0 Cm0 C *M C *m De Dp Dm Dm0 D *m FA Im I *m km Km K *m Nm N *m rM rm Rp Rm T Greek a(Cm) a* ab s1 b* oM g s j h t ts ts 1
Biot number (=kmRp/De) macropore concentration microsphere concentration (=Im(Cb)) bulk concentration during adsorption, and C0 during desorption initial (t= 0) macropore concentration microsphere concentration in equilibrium with Cb (=CM/Cb) for adsorption, (= CM/C0) for desorption (=Cm/Cmb) effective diffusivity particle macropore diffusivity (=De/oM) micropore diffusivity [= Dm(Cb)] dimensionless surface diffusion coefficient fractional adsorption a suitable micropore adsorption isotherm dimensionless isotherm external mass transfer coefficient Langmuirian constant (=Km×Cb) for adsorption, (= Km×C0) for desorption micropore flux (NmRm/[CmbDm(Cb)]) particle radical coordinate microsphere radial coordinate radius of particle radius of microsphere time coordinate letters local relaxation time (= a/ab) [= a(Cmb)] [= 3(1−s)g)] [= Dpsab/R 2p] macroporosity [= Dm0R 2p/(DpR 2ms)] {=oMCb/[oMCb+(1−oM)Cmb]} (= rm/Rm) (=rM/Rp) (= sDpt/R 2p) (=t/s) (=s1ts )
F. Liu et al. / Computers and Chemical Engineering 24 (2000) 1981–1995
1995
Acknowledgements This research has been supported by the Australian Research Council through the Large Grant Scheme. References Abarzhi, I. I. (2000). Wave mechanism of mass transfer for kinetics of adsorption in biporous media. Colloids and Surfaces A, 164, 105–113. Bermudez, A., & Vazquez, M. E. (1994). Upwind methods for hyperbolic conservation laws with source terms. Computational Fluid, 23, 1049 – 1071. Bhatia, S. K. (1987). Modelling the pore structure of coal. American Institute of Chemical Engineering Journal, 33, 1707 – 1718. Bhatia, S. K. (1997). Transport in bidisperse adsorbents: significance of the macroscopic adsorbate flux. Chemical Engineering Science, 52, 1377 – 1386. Bhatia, S. K., Gray, P. G., & Do, D. D. (1991). Modelling of sorption of gaseous sorbates in bidispersed structured solids. Gas Separation and Purification, 5, 49 – 55. Finlayson, B. A. (1972). The method of weighted residuals and 6ariational principles. Academic Press. Gray, P. G., & Do, D. D. (1990). Adsorption and desorption of gaseous sorbates on a bidispersed particle with Freundlich isotherm-III contribution of surface diffusion to the sorption dynamics of sulphur dioxide on activated carbon. Gas Separation and Purification, 4, 149–157. Godlewski, E., & Raviart, P. A. (1991). Hyperbolic systems of conservation laws, Mathematiques & Applications, N. 3/4, Ellioses-Edition Marketing. LeVeque, R. (1990). Numerical methods for conser6ation laws. Basel: Birkha¨user. Liu, F., & Bhatia, S. K. (1999). Computationally efficient solution techniques for adsorption problems involving steep gradients in bidisperse particles. Computers and Chemical Engineering, 23, 933–943. Liu, F., & McElwain, D. L. S. (1997). A computationally efficient solution technique for moving-boundary problems in finite media. IMA Journal of Applied Mathematics, 59, 71–84. Liu, F., & Zheng, X. (1986). Uniformly convergent difference scheme for the singular perturbation of a self adjoint elliptic partial differential equation. Applied Mathematical Mechanics, 7, 495–504. Liu, F., McElwain, D. L. S., & Donskoi, E. (1998). The use of a modified Petrov – Galerkin method for gas – solid reaction modelling. IMA Journal of Applied Mathematics, 60, 1–14. Luikov, A. V. (1980). Heat and mass transfer. Moscow: Mir Publishers. Petzold, L. R. (1982). A Description of DASSL: A Differential/Algebraic Equation System Solver, Sandia Technical Report SAND82-8637 Livermore, CA USA. Ruckenstein, E., Vaidyanathan, A. S., & Youngquist, G. R. (1971). Sorption by solids with bidisperse pore structures. Chemical Engineering Science, 26, 1305 – 1318. Ruthven, D. M. (1984). Principles of adsorption and adsorption processes. New York, USA: Wiley. Sohn, H. Y., & Szekely, J. (1972). A structural model for gas– solid reactions with a moving boundary-III. A general dimensionless representation of the irreversible reaction between a porous solid and a reactant gas. Chemical Engineering Science, 27, 763 – 778. Villadsen, J. V., & Stewart, W. E. (1967). Solution of boundary-value problems by orthogonal collocation. Chemical Engineering Science, 22, 1483 – 1501.
.
.