Finite element modelling of multisolute activated carbon adsorption

Finite element modelling of multisolute activated carbon adsorption M. Akram Hossain Dept. of Civil Engineering

and Construction,

Bradley University, Peoria, IL, USA

David R. Yonge Department of Civil and Environmental Pullman, WA, USA

Engineering,

Washington

State University,

Multisolute model@ of activated carbon adsorption in afixed bed reactor

can be complex and expensive computationally. The computational techniques that have been utilized over the years are jmite differences and a global method of orthogonal collocation. The global method of orthogonal collocation, which is currently gaining popularity over the more traditional finite difference techniques, has some inherent disadvantages. In an attempt to eliminate these disadvantages, finite element models are developed herein by using the characteristic-Galerkin method and the Pertov-Galerkin method. Both the models were found to be sensitive to changes in film transfer coefficients but were relatively insensitive to changes in dtffusion coefficients. The models were alsofound to be stable and computationally efficient. The model predictions were compared to bisolute adsorption data, and it was found that the PertovGalerkin model gives a better description of experimental data and, in addition, gives convergent results using a significantly larger time step than the characteristic-Galerkin method.

Keywords:

activated carbon, modelling, multisolute adsorption, finite element

Introduction Activated carbon adsorption in a fixed bed reactor (FBR) has been found to be an effective technology for removing a wide variety of hazardous organic contaminants from water and wastewater. l-4 The process leading to the design of such an adsorption system is frequently time consuming and, as a result, expensive.2.s*6 Mathematical process modelling is considered to be an effective tool that can facilitate the design of FBR by reducing the number of pilot scale tests required to evaluate various operating conditions and develop design parameters. lm4 Obviously, a model that yields accurate predictions under a wide variety of system conditions and is computationally efficient is highly desirable. Current modelling techniques, however, fail to satisfy these criteria. The global method of orthogonal collocation (GMOC) is currently favored over finite difference methods in modelling carbon adsorption in FBRs.‘~‘-‘~ Address reprint requests to Professor Yonge at the Department of Civil and Environmental Engineering, Washington State University, Pullman, WA, 99164-2910, USA. Received 23 October 1990; revised 23 April 1992; accepted 5 May 1992

630

Appl. Math. Modelling,

1992, Vol. 16, December

In the GMOC a trial solution is constructed by using orthogonal polynomials as the basis function. The solution is then substituted into the governing differential equations to obtain the residual, which is made to vanish at the collocation points. These collocation points are fixed and are determined by the roots of the highest-degree polynomial that is a member of the space spanned by the trial solution. The fixity of the collocation points may lead to an inaccurate determination of the field variable in regions where it is changing rapidly. ‘I It has been reported that if one attempts to take arbitrarily spaced collocation points in regions of rapid change, the method loses its efficiency.i2 The GMOC is usually employed in conjunction with integration routines such as GEAR7s8.13or semi-implicit Runge-Kutta methods.‘4.‘5 These integration routines are relatively costly and require small time steps to achieve desired levels of accuracy.‘5-‘7 Moreover, the discretization matrices arising from the use of GMOC do not have any computationally advantageous properties. The matrices are neither banded nor symmetric. This adds to computational inefficiency and makes stability analysis difficult. ‘* The GMOC has also been shown to result in oscillatory results.” Furthermore, researchers in modelling the carbon adsorption in an FBR are confronted with a

0 1992 Butterworth-Heinemann

Multisolute activated carbon adsorption: M. A. Hossain and D. R. Yonge numerical problem of getting negative solid phase concentrations.‘4~‘s~20 The negative solid phase concentration causes the algorithm to break down because real powers of a negative quantity are undefined. The breakdown of the algorithm is commonly avoided by arbitrarily setting the negative concentrations to zero.*OThis is not desirable because it violates the principle of local equilibrium and is physically incorrect. In single-solute modelling, an asymmetric finite element method (FEM) has been found to be effective in eliminating the negative solid phase concentrations, a result of oscillation in the solution of the liquid phase equation.*’ The liquid phase equation is an advective transport equation, and certain FEMs have been shown to be efficient in suppressing oscillations inherent to the solution of such an equation.‘1,22-28 Moreover, the use of an FEM usually results in a straightforward stability analysis and faster convergence.29 Additionally, the use of piecewise linear interpolating functions (PLIFs) results in tridiagonal discretization matrices that are often symmetric and positive definite, making the algorithm efficient and stable. This paper presents the development of two finite element models employing the characteristic-Galerkin (CG) and the Pertov-Galerkin (PG) finite element methods for multisolute adsorption in a FBR by using PLIF and compares their performances.

Equation (1) describes the concentration change in the liquid phase, and equation (2) describes that in the solid phase. Equation (2) is commonly referred to as the homogeneous surface diffusion model (HSDM). The HSDM has been reported to be successful in predicting the breakthrough curves for many single-solute and multisolute adsorption systems.5v6*8 The liquid phase mass balance equation (equation 1) is subject to the following initial and boundary conditions: t = 0,

O
tro,

z = 0,

The equations that describe adsorption dynamics in a FBR are developed from a mass balance. If one assumes that plug flow occurs through the FBR, the rate of adsorption is diffusion limited, film transfer describes the liquid diffusion resistance that occurs at the external adsorbent surface, the carbon particles are homogeneous spherical solids, and surface diffusion is the major mode of intraparticle transport, then one obtains the following liquid and solid phase equations for the ith component of a solution containing N solutes*‘: dCi -= at

- “gQJ1-E -(C; az R E

1

asi___ Dsi a ,*!!s? at

r*

ar [

dr

- C,i)

(1) (2)

where C = liquid phase concentration of the solute [M/L31 r= time [T] ll= velocity of flow [L/T] distance measured along the axis of the Z= FBR [L] kf = film transfer coefficient [L/T] R= radius of a carbon particle [L] E= porosity of the fixed bed c, = concentration of solute in the boundary layer [M/L’] 4= solid phase solute concentration [M/M] D, = surface diffusivity of the solute [L*/T] r= radial distance from the center of a carbon particle [L]

c; = 0 Ci= Cp

(3) (4)

where L = length of the FBR [L] CO = concentration of the solute at the entrance of FBR The HSDM (equation 2) is subject to the following initial and boundary conditions: t = 0,

Osr
tro,

r = 0,

9_ -0 3r

tzo,

r = R,

DJJ~

r = R,

The mathematical model

L,

qj = 0

C,i = f(q,,

3 qs2,

= k,;(Ci - C,J . * . 3 c?sN)

(5)

(7) (8)

where p = apparent density of the carbon particles [M/L31 qS = concentration of solute on the surface of the carbon particles [M/M] Equation (8) couples the intraparticle transport equation, equation (2), with equation (1). In adsorption terminology, equation (8) is called the isotherm equation or equilibrium model. There are numerous single-solute and multisolute equilibrium formulations, both theoretical and empirical in nature, that vary considerably in applicability and utility. l4 Among the single-solute models the Freundlich isotherm equation often yields the best description of experimental data.4 Consequently, this semiempirical model is frequently used in modeling FBR adsorption dynamics and is shown below for reference: qsi = t&k

(9)

where K and n are characteristic constants. Ideal adsorbed solution theory (IAST) is frequently used to describe multisolute equilibria.30 The model requires single-solute isotherm data for the prediction of multisolute equilibria. If the single-solute isotherm data are defined by the Freundlich equation, then the IAST can be represented by the following equation3’:

(10)

Appt. Math. Modelling,

1992, Vol. 16, December

631

Multisolute

activated

carbon

adsorption:

M. A. Hossain

Equation (10) was used to couple the liquid and solid phase equations.

and D. R. Yonge

defined in equations (1 l)-(22) results in the following nondimensional equations and initial and boundary conditions:

Nondimensional model

JG=

In numerical computation it is always advantageous to deal with normalized variables. Normalization of the variables in the mathematical model describing the adsorption dynamics in an FBR reduces the number of independent parameters to be dealt with. The following variables were introduced to nondimensionalize the model equations:

dT

a~;-

- Dp~

T= 0,

o
TzO,

2 = 0,

@i -=

E,ii

8T

72

3D,Sti(Ci

1,

(23)

- C,;)

C; = 0 c;=

(24)

1

(25)

aT; [ a7 1 ---3

(26)

T= 0,

osr<

T2 0,

F=O,

@i -_=O ar

Tz 0,

‘i;= 1,

z

1,

(11)

4; = 0

(27 ) (28)

(12)

z!_

qi =

= Sh;(c;

(13)

$,

9ei

2 D

1-

=

r=

4&

Ei=l

(103>

--

Is

E

c

1,

qeii?si

.=

”

(14)

5

- es;>

njqejGsj

n;K; 4ej4sj

j=l

[

$&co

(29)

ai

I

(30)

Selecting the interpolating function

(17)

22

(18)

L

An interpolating function defines the variation of the field variable over an element. To ensure interelement continuity, the interpolating functions selected must satisfy compatibility and completeness requirements.32 An element of at least Co continuity is necessary in this case to meet the above-mentioned criteria. The PLIF satisfies these requirements and results in tridiagonal coeffkient matrices. Consequently, the PLIF was used herein. The following is the equation of a PLIF: 0 X-Xi-1

(19) 4; = Rk,iCP Sh; = qeiD,ip(10-3) St,

=

(20) (21)

RE

qei =

K;Cp;

The application of the chain rule of differentiation to the model equations in conjunction with the variables 0

1

x - xi-1 _ 3a (x -

44=

x; - xi-, x, rtl

--x

xi + 1 - xi

xi-

3a

(X

-

Appl. Math. Modelling,

-

xi)

x5x;_,

if

X;_l(XlX;

XJCx

-

xi + 1) xi)*

if

~~_~Ix--(x~

x;

if

x~SXIX~+~

if

x5x;+,

(31)

X;

PLIF as the weight function in the liquid phase equation oscillatory results, which can as the basis of the trial soluthe asymmetric basis function

(32)

Xj_1,12

(Xi + 1 -

b

632

l)Cx

(X; +

xi+1 Xi+l -

The use of an asymmetric for the spatial derivative (equation 1) can prevent occur when using PLIFs tion.22*25The equation of is given below:

if

X5X;-I

Xi-1

! 0

k,;T(l - l1

I

Xi-

if

if

~~5x5~~~~

if

x2x;+,

1992, Vol. 16, December

Multisolute

activated

carbon adsorption:

Here, LYis a free parameter that defines the asymmetricity and &?-educes to the standard PLIF when (Yis set to zero. When an asymmetric PLIF is used as a weight function, the method is referred to as the PG method.

TrO,

M. A. Hossain and D. R. Yonge

v= 1,

5

= Sh;(l$ - lZsi)

s

njqejGsj

j=l

Characteristic-Galerkin formulation

QK;

i

The CG formulation for the multisolute adsorption of activated carbon adsorption in an FBR is given below: - 30, Sri (C; - C,J (33)

1,

pi = 0

T= 0,

o
TzO,

z = 0,

ci = 1

TrO,

I?= 1,

s= a2

(34) (35)

Edi -- a 723av2 a7 [

T= 0,

osr<

TrO,

‘7:= 0,

[Ah]

1,

(37) qi = 0

(38)

aa -0 z-

+ 3St, Dg[ AhI

1 ni

(41)

In equation (33), $is aconstant, A Tis the time step, and the second quantity on the right-hand side is an artificial diffusion term. If $ = 1, then the formulation is equivalent to the CG formulation, and Cc,= 0 reduces it to the original formulation. The CG formulation does not change the solid phase equation because a smooth solution of the solid phase equation can be obtained by using the PLIF.

Finite element formulation

O

31

aqi -= aT

(40)

+ Dg[ Bh]{c;k} - $y[

There are a number of finite element methods.32,33 The semi-discrete Galerkin finite element method (SDGFEM) was used to discretize the partial differential equations (PDEs) in the mathematical model because of its versatility. The application of the Galerkin principle to the model PDEs and appropriate mathematical manipulations resulted in the following ordinary differential equations (ODES) in time:

k= 1,2,...,N

Eh,{c’,} = (03

+ Ed, ShJ IF’]{&} = Ed,Jhk (c’k,NR - csk){DZ)

(42)

k= 1,2,...,N

(43)

where I ai =

$i (Zl4j (T)dZ

I

1,2 , . * f 3NA

i,j=

(44)

0

r2+i @)+j @)dF

a$ =

i,j=

I,2 ,...,

NR

(45)

i,j=

1,2 ,...,

NA

(46)

0

b;. =

’

j = 1,2, . . . , NA if if

i,j = 1,2,. . . ,NR (47 ) i,j = 1,2,. . . , NR (48)

i#NR i=NR

(49)

(50)

Here NA + 1 is the number of nodes along the axis of the FBR and NR is that along the radius of a carbon particle. The above formulation becomes the CG formulation when + = 1 and ‘pi = +j. For asymmetric finite element formulation, 4 = #and Cc,= 0. The problem now reduces to solving equations (42) and (43), coupled through equation (41), with the conditions defined in equations (34) and (38).

Appl. Math. Modelling,

1992, Vol. 16, December

633

Multisolute activated carbon adsorption: M. A. Hossain and D. R. Yonge

Solution technique The system of ODES obtained after the application of Galerkin principle is stiff in nature. It was determined that an implicit backward discretization of the temporal derivative would result in a stable algorithm for such a system. The following equations are the results of the implicit temporal discretization: [L;]{c+

‘> = [ k@{~} + {Di}

lIJ%1{q”,+‘~ = [MmzI

(51)

+ m1

(52)

where l& = (1 + 3St,D,AT)a$

D2A T + DBATbi - I)+ ez (53)

m$= ai

(54)

lfij = a$.+ Ed,A Tb$

(55)

m$= a$

(56)

Here, s is an index for the time step. The matrices [ Li] and [ L$] do not change with time and were decomposed into upper and lower triangular matrices by employing the principle of lower and upper (LU) decomposition.34 At each time step the systems of equations was then solved by forward and backward substitution.34

Results and discussion Computational effort The use of PLIF as the basis function provided two distinct advantages: the discretization matrices were tridiagonal, and the elements of the matrices were determined exactly by analytical integration. Tridiagonal discretization matrices reduce the computational effort necessary to solve the model equations, and the analytical integration saves computational time in obtaining the element matrices. The discretization matrices were also constant, and as a result, it was necessary to compute them only once. The CG formulation differs from the PG formulation in that the former requires the determination of an additional matrix, [ Eh]. The matrix [ Eh] is tridiagonal, symmetric, and positive definite, and only three floating point arithmetic operations are

necessary to compute all of its elements. The two formulations therefore require essentially the same amount of computational effort with regard to discretization. The vast majority of the computational effort is spent in solving the fully discretized equations. The solution technique, as discussed in the preceding section, consists of one-time LU decomposition of the matrices [ Li] and [ L$f] and forward and backward substitution at each time step for each organic contaminant that constitutes the multisolute mixture. Since the matrices are tridiagonal, the computational effort necessary to perform all these operations is of order N, N being the number of equations in the system. In addition, it is necessary to compute the equilibrium liquid phase boundary layer concentration from the solid phase surface concentration by employing equation (41). This operation requires approximately 3N + 2 floating point multiplications, 2N floating point additions, and one exponentiation. The CG formulation, however, utilizes a smaller time step than the asymmetric formulation. The CG formulation was found to provide convergent results for A T = 10A5, which is an order of magnitude smaller than AT = 10e4 required by PG formulation to converge. So the PG formulation is computationally less intensive than the CG formulation. Stability analysis The stability of the algorithms was checked in two different ways: by a matrix method and by the von Neumann procedure. The matrix method states that the algorithm is stable when the spectral radii of the matrices [Lt]-‘[Mi] and [Lg]-‘[Mi] are less than or equal to unity. The spectral radii of the matrices were determined for varying time steps and number of nodes. It was found that the spectral radii were always less than or equal to unity for the input parameters given in Tables 1 and 2. This method, however, is specific to the input parameters used. In other words, it is necessary to perform the analysis for each set of input parameters. In an attempt to derive a generalized stability criterion, the von Neumann procedure was used to get an expression for the amplification factor. The amplification factor should be less than or equal to unity for the algorithms to be stable. The amplification factors that were obtained by applying von Neumann procedure are given below:

*; - =

(2 + 6StD,AT)LR + lZ~D,Cos’

PAZ +6aD 2

R

(CC= 1 + 3StD,AT+

634

Appl. Math. Modelling,

3

R+3

1992, Vol. 16, December

(57 1 kTY2-6LTD “AT AT

’

Multisolute

activated

carbon adsorption:

fer coefficient, k,. Tables I and 2 contain the reference values of the input parameters used to study the impact of varying mass transfer parameters on model predictions. The impact of varying the intraparticle diffusion coefficient on model prediction for solute 1 is presented in Figure 2. The diffusion coefficient was varied by an order of magnitude from the reference values, and everything else was kept constant. For the two higher diffusion coefficient values, the breakthrough curves are essentially the same. For the smallest value, however, an earlier breakthrough is observed. This suggests that for a given film transfer rate there is a limiting value of the diffusion coefficient beyond which the model is relatively insensitive to additional increases in the diffusion coefficient. The curves in Figure 3, which show the impact of the diffusion coefficient on model predictions for solute 2, reveal the same trend as was observed in Figure 2 for solute 1. The curves in Figures 4 and 5 depict the impact of varying film transfer resistances on the model prediction for solutes 1 and 2, respectively, while maintaining the other parameters at their reference values. The film transfer coefficients were varied by ? 10% from the reference values. Although the initial portions of the breakthrough curves are not affected, it is observed that at later times a higher film transfer

Table 1. Reference values for isotherm and mass transfer parameters and feed concentrations for sensitivity analysis Solute 2

Solute 1

Parameter K

2. 0.12 7.2 x 10m3 cm/s 1.20 x 10e8 cm*/s 0.14 mmol/L

2.70 0.18 6.8 x 1O-3 cm/s 1.30 x 10-s cm*/s 0.16 mmol/L

n

4 DS c0

Table 2. Physical and discretization parameters sensitivitv analysis Value

R

0.017 cm 0.80 g/cm3 0.0001 0.60

P

:T ;vA

44.0 25 s

NR

10

where tpG = amplification formulation tc” = amplification formulations R=

for

Parameter

1 + 2 co??

factor for the PG factor for the CG PAZ

C

Y=

Cos (y)-

P=

a real number

M. A. Hossain and D. R. Yonge

It can be shown that the conditionally stable because unity regardless of the choice finite element formulation is St,Azl.

1

tn

1.6

(y)

0.6

CC formulation is ungCG is always less than of p and the asymmetric stable when LY~ % 11 -

0.6

0

On asymmetricity

0

Asymmetricity in the PG Iinite element formulation is introduced through the parameter CLThe intent is to reduce oscillations in the solution. However, numerical analysts are concerned about the numerical dispersion that might be introduced as a result of using asymmetric basis functions. The extent of numerical dispersion was studied by obtaining spatial solute concentrations at varying degrees of asymmetricity. The value of (Ywas varied from 0.5 to 1.Owith an intermediate value of0.75. The curves in Figure I provide the graphical representations of the results obtained for T = 0.3. It is obvious that the asymmetricity introduces very little numerical dispersion. So any value of u up to 1.OOcan be used to reduce the extent of oscillations in the solution of the liquid phase mass balance equation as a means of avoiding breakdown of the PG finite element algorithm.

0.75

0.50

1 .oo

: Figure 1.

Effect of varying a on spatial oscillation

,

2.0

1.6

1.2 G 0.6

----D,, ....._- 4, -D,,

0

Model sensitivity to mass transfer parameters

There are two mass transfer parameters: the intraparticle diffusion coefficient, D,, and the film trans-

0.25

Figure 2.

Appt.

0.4

0.6

=1.3x10”cm2/s =l .3x10~‘cm2 la =l .3x10’ecm2/s

1.2

1.6

2.0

Effect of varying Ds on model prediction for solute 1

Math. Modelling,

1992, Vol. 16, December

635

Multisolute

activated

carbon adsorption:

IV. A. Hossain and D. R. Yonge

2.0

1 ---.... ....... D, D,=l.3x10”cm2/s =, ,3x, 9scm*/s

1.6 -

0, =1.3x10~0cm2/.s

1.2 5, 0.6 -

0

0.4

0.6

, 1.6

1.2

2.0

T Figure 3.

Effect of varying 0, on model prediction for solute 2

2s I

Comparison of model predictions with experimental data The experimental results presented here were obtained from published literature.i4 Table 3 contains the mass transfer and isotherm parameters, and Table 4 contains the physical and discretization parameters that were used as model inputs. Figure 6 provides a graphical comparison between the model predictions and the experimental data points for a bisolute mixture of pchloro phenol (PCP) and dodecylbenzenesulfonate (DBS). The PG finite element formulation is found to furnish a better description of the experimental data. The significant differences in the prediction of the two models may be explained by the fact that the CG formulation adds an additional term to the liquid phase PDE that is inconsistent with the mathematical model. The failure of the models to predict the experimental results accurately can be attributed to the fact that the IAST cannot account for solute-solute interaction, un-

Table 3. Isotherm and mass transfer parameters, concentrations for PCP and DBS4

1.5 E*

Parameter

1.0 -

K n kf

0.3

0.6

0.9

Solute 2 (DEW

2.62 0.15 4.35 x 10m3 cm/s 1.42 x 10mg cm2/s 0.1642 mmol/L

D* C0 “0

Solute 1 (PCP)

and feed

1.61 0.31 2.03 x 10e3 cm/s 0.36 x 10-gcm2is 0.1722 mmol/L

1.2

T Figure 4.

Table 4. Physical and discretization parameters PCP and DBSq4

Effect of varying k, on model prediction for solute 1

jl

jj

1.0 -

0.5 0

2.0

0

Parameter

Value

R P

0.018 cm 0.99 g/cm3

;vA

0.47 0.0001 (PG) 0.00001 (CG) 30.1 25 s

NR

IO

;\T

Ee

0.3

0.6

0.9

1.2

for

,

I

1.6 -

T Figure 5.

Effect of varying k, on model prediction for solute 2

coefficient results in smaller effluent concentrations. This is explained by the fact that a higher film transfer results in a higher concentration gradient at the surface of a carbon particle, which in turn results in a greater rate of adsorption and a lower liquid phase concentration. Although the above discussion makes particular reference to the results obtained from the CG formulation, the same was found to apply to the PG formulation.

636

Appl.

Math. Modelling,

1992, Vol. 16, December

Figure 6. Comparison mental data

between

model predictions

and experi-

Multisolute

activated

equal competition, and irreversibility. These factors were accounted for by Liangr4 by using an empirical Freundlich-Langmuir type isotherm model and determining the isotherm parameters by fitting the model with actual bisolute experimental data. Obviously, the isotherm parameters have a significant impact on the model predictions. So a better tit could be obtained if IAST were modified to account for irreversibility, unequal competition, and solute-solute interaction.

carbon adsorption:

11 12 13 14 15

Conclusions 16

In light of the discussion presented in the preceding section it can be concluded that the CG and the PG models are very sensitive to changes in isotherm parameters and coefficient of film transfer and are relatively insensitive to changes in the diffusion coefftcient. The models are also stable and convergent and are not computationally intensive. The computational effort necessary to solve the model equations is of order N. The PG formulation is computationally less intensive than the CG formulation because the former was found to provide convergent results for A T = 10m4whereas the CG formulation required AT to be 10p5. We finally note that the PG finite element formulation yields a better description of the experimental breakthrough curves of a bisolute mixture of PCP and DBS.

17

18

19 20

21

22 23

References 1

2

Crittenden, J. C. and Hand, D. W. Modeling of adsorption, desorption and displacement in fixed bed adsorbers. Fundamental of Adsdrption: Proceedings of the Engineering Foundations Conference, Schloss. Elmau, Bavaria, West Germany, May 6-ll”, 1983, pp. 185-194 Crittenden, J. C., Luft, P., Hand, D. W. and Friedman, G. Predicting the removal of known organic compounds in unknown mixtures and total organic halogens in fixed-bed adsorbers. Proceedings of the 1984 Specialty Conference, ASCE Env. Enara. Div.. Univ. of Southern Calif.. June 25-27. 1984. pp. 370--373 Thacker, W. E., Crittenden, J. C. and Snoeyink, V. L. Modeling of adsorber performance: Variable influent concentration and comparison of adsorbents. J. WPCF 1984,56(3), 243-250 Weber, W. J., Jr. and Smith, E. H. Criticalreview-Simulation and design methods for adsorption systems. Environ. Sci. Technol. 1987, 21(1 l), 1040-1049 Lee, M. C., Crittenden, J. C., Snoeyink, V. L. and Ari, M. Design of carbon beds to remove humic substances. J. Environ. Engrg. Div. ASCE 1983, 109(3), 631-645 Crittenden, J. C. and Weber, W. J., Jr. Model for design of multicomponent adsorption systems. .I. Environ. Engrg. Div. ASCE 1978, 104(EE6), 1175-l 195 Crittenden, J. C., Luft, P. and Hand, D. W. Prediction of fixedbed adsorber removal of organics in unknown mixtures. J. Environ. Engrg. ASCE 1987, 113(3), 486-498 Hand, D. W., Crittenden, J. C. and Thacker. W. E. User-oriented batch reactor solutions to homogeneous surface diffusion model. J. Environ. Engrg. Div. AS& 1983, 109(l), 82-101 Speitel, G. E., Jr.. Lu, C.-J.. Zhu. X. J. and Turakhia. M. H. Biodegradation and adsorption of’s bisolute mixture in GAC columns. J. WPCF 1989, 61(2), 221-229 Weber, W. J. and Liu, K. T. Determination of mass transport

24

25

26

I

3 4 5 6 7

8 9 10

27 28

29 30 31

32 33 34

M. A. Hossain and 0. R. Yonge

parameters for tixed-bed adsorbers. Chem. Engrg. Sci. 1980,6, 49-60 Carey, G. F. and Finlayson, B. A. Orthogonal collocation on finite elements. Chem. Engrg. Sci. 1975, 30, 587-596 McGowin, C. R. and Permuller, D. D. Regions of asymptotic stability for distributed parameter systems. Chem. Engrg. Sci. 1971, 26, 275-286 Crittenden, J. C., Wong, B. W. C., Thacker, W. E. and Hinrichs, R. L. Mathematical model of sequential loading in fixed-bed adsorbers. J. WPCF 1980, 52(1 l), 2780-2795 Liang, S. Mathematical modeling of sorption in heterogeneous system. Ph.D. dissertation, Univ. ofMichigan, Ann Arbor, 1984 Liu, K. T. Determination of mass transport parameters for fixed-bed adsorbers by a micro-column technique. Ph.D. dissertation, Univ. of Michigan, Ann Arbor, 1980 Bui, T. D. and Bui, T. R. Numerical methods for extremely stiff differential equations. Appl. Math. Modelling 1979,3,335-358 Hindmarsh, A. C. Linear multistep methods for ordinary differential equations: Method formulation, stability, and the methods of Nerderieck and Gear. Rept. UCRL-51186 Rev. 1, Lawrence Livermore Laboratory, Univ. of Calif., Livermore, 1972 Villadsen, J. and Sorensen, J. P. Solution of parabolic partial differential equations by double collocation method. Chem. Engrg. Sci. 1969, 24, 1337-1349 Neretnieks, I. Adsorption in finite bath and counter current flow with systems having a concentration dependent coefficient of diffusion. Chem. Engrg. Sci. 1976, 31, 465-471 Thacker, W. E. Modeling of activated carbon and coal gasitication char adsorbents in single-solute and bisolute solute systems. Ph.D. thesis, Univ. of Illinois at Urbana-Champaign, 1981, Hossain, M. A. and Yonge, D. R. Finite element model of single-solute activated carbon adsorption. J. Environ. Engrg. Div. ASCE 1992, 118(2), 238-252 Huyakorn, P. S. An upwind finite element scheme for improved solution of the convective-diffusion equation. Res. Rept. 76WR. 2, Dept. of Civil Engrg., Princeton Univ., 1976 Huyakorn, P. S. Solution of steady state convective transport equation using an upwind finite element scheme. Appl. Math. Model@ 1977, 1, 187-195 Huyakorn, P. S. and Nilkuha, K. Solution of transient transport equation using an upwind finite element scheme. Appl. Math. Modelling 1979, 3, 7-17 Lapidus, L. and Pinder, G. F. Numerical Solution of Partial Differential Equations in Science and Engineering. John Wiley, New York, 1982 Lee, J. H. W., Peraire, J. and Zienkiewicz, 0. C. The characteristic-Galerkin method for advection-dominated problemsAn assessment. Comput. Methods in Appl. Mechanics and Engrg. 1987, 61, 359-369 Lucas, T. R. and Reiddien, G. W., Jr. Some collocation methods for nonlinear boundary value problems. SZAM J. Numer. Anal. 1971, 9(2), 341-356 van Genuchten, M. T. On the accuracy and efficiency of several numerical schemes for solving the convective-dispersive equation. International Conference on Finite Element in Water Resources, Princeton Univ., 1976 Russel, R. D. and Shampine, L. F. A coliocation method for boundary value problems. Numer. Math. 1972, 19, l-28 Radke, C. J. and Prausnitz, J. M. Thermodynamics of multisolute adsorption from dilute liquid solutions. AIChE J. 1972, 18, 761-768 Crittenden, J. C., Luft, P., Hand, D. W., Oravitz, J. L., Loper, S. W. and Ari, M. Prediction of multicomponent adsorption equilibria using ideal adsorbed solution theory. Environ. Sci. Technoi. 1985, 19(11), 1037-1043 Huebner, K. H. and Thornton, E. A. Finite Element Methods for Engineers. 2nd ed., John Wiley, New York, 1982 Zienkiewicz, 0. C. The Finite Element Method. 3rd ed., McGraw-Hill, New York, 1982 Burden, R. L. and Faires, J. D. Numerical Analysis. 3rd ed., PWS, Boston, 1985

Appl. Math. Modelling,

1992, Vol. 16, December

637