I n d . Eng. Chem. Res. 1989, 28, 108-119
108
Facilitated Transport with Unequal Carrier and Complex Diffusivities O s m a n A. Basaran,*,+Peter M. B u r b a n , a n d Steven R. Auvil Air Products & Chemicals, Inc., P.O. Box 538, Allentown, Pennsylvania 18105
Operation of facilitated transport membranes is analyzed to answer the long-standing question of the effect of differing carrier and complex diffusivities, DB and D, on system performance. For the reaction scheme A B e AB, where A is the permeant, B is the carrier, and AB the complex, the governing equations are either two coupled differential equations augmented by an integral constraint (method 1) or two coupled integro-differential equations (method 2). The governing nonlinear equations are solved by asymptotic analysis and also, numerically, by the Galerkinlfinite element method. Steep boundary layers that arise for operation near equilibrium are resolved numerically by means of an adaptive mesh. Method 1leads to an arrow-shaped matrix, which can be LU decomposed more efficiently than the full matrix resulting from method 2. The efficient algorithm arising from method 1makes it possible to analyze situations in which DB # DM a t roughly the same computational cost as previous algorithms that are restricted to situations in which DB = Dm-a boon when performing repeated engineering evaluations. The results show that relative t o the case in which the ratio of the diffusivities, R = DB/DAB,is equal t o one, the facilitation factor-the flux of A divided by the Fickian flux of A-is increased when R < 1 and is decreased when R > 1. Furthermore, an optimization is carried out to determine the maximum facilitation factor and the corresponding dimensionless equilibrium constant, thereby extending the results of Kemena e t al. to the situation in which DB # DAB.
+
1. Introduction A facilitated transport membrane is a substrate that is filled with a solution of a reactive molecule called a carrier. The carrier reacts reversibly with a permeant, for example, a gas in a mixture of gases to be separated. First, all gases pass through such a membrane, albeit a t different rates, by the mechanism of solution/diffusion. However, the reactive gas, A, associates with the carrier, B, on the high-pressure side, and AB, or complex, also diffuses through the membrane. The reactive gas then disassociates from the complex on the low-pressure side. This chemical reaction enhanced diffusion can lead to a high permeability to thickness ratio and a high selectivity. Facilitated transport has received a great deal of attention in the literature because of its importance in life sciences and its potential in industrial separations (see, e.g., Ward (1970), Goddard et al. (1970, 1974), Schultz et al. (1974), and Kemena et al. (1983)). Almost all of the attention has been focused on the reaction scheme kf
A + B Z A B kr
(1)
where kf and k, are the forward and reverse rate constants. Furthermore, all except three of these previous analyses (see below) have been restricted to the situation in which the diffusivity of the carrier, Dg, and that of the complex, DAB, are equal. Unfortunately, analyses of facilitated transport for which DB # D m are insufficient to draw any general conclusions, as shown below. A goal of this paper is to remedy this situation. Interestingly, it is not obvious at the outset what the effect of the ratio DBIDm would be on transport through a membrane because the molar fluxes of B and AB are equal in magnitude but opposite in direction. Aris (1975), in his treatise on catalysis, has analyzed by means of regular perturbation analysis facilitated transport
* Author to whom correspondence
should be addressed. Present address: Oak Ridge National Laboratory, Oak Ridge, T N 37831. Oak Ridge National Laboratory is operated by Martin Marietta Energy Systems, Inc., under Contract DE-ACOS840R21400 with the US Department of Energy. 0888-5885/89/2628-0108$01.50/0
when DB # DAB in the limit in which the importance of diffusive flux is large compared to the rate of reaction, Le., in the limit as the Damkohler number P 0 (see section 2 ) . His analysis is complete to first-order in P. Earlier, Smith et al. (1973) had analyzed the situation in which Dg = Dm in the limit as P 0; their perturbation expansions included terms that are of second-order in P. In section 3, it is shown that the perturbation expansion for the facilitation factor-a measure of the enhancement of the diffusion of A due to chemical reaction-does not involve the ratio DB/Dm until it is second-order in P. That is, the results of Aris (1975) and Smith et al. (1973) are identical up to and including terms first-order in P, and Aris’s (1975) analysis has to be extended to second-order in P in order to understand the effect of the ratio DB/Dm on the flux. Kaper et al. (1980) have analyzed by means of singular perturbation analysis facilitated transport when DB # D m in the limit in which the importance of the rate of reaction m. Though is large compared to diffusive flux, Le., P their elegant analysis is, for the most part, correct, their conclusions are in error. This too is remedied in this paper. The general case, for which P is neither small or large, has been considered briefly by Jain and Schultz (1982), who analyzed the problem numerically. Unfortunately, Jain and Schultz (1982) have drawn no general conclusions as to the effect of the ratio DBIDAB on the solutions they have generated. Moreover, the algorithm they have used is moderately efficient when DB = DAB but is inefficient when DB # DAB-see section 6. Section 2 presents the equations and boundary conditions that govern the steady-state operation of facilitated transport membranes. It is shown that there are two alternative forms of the governing equations. The first one is a standard boundary-value problem comprised of two simultaneous nonlinear differential equations, but is augmented by an integral constraint. The other one is two simultaneous nonlinear integro-differentialequations. The second formulation has not previously appeared in the literature. In section 3, the governing equations are solved by means of regular perturbation analysis. In section 4,the analysis of Kaper et al. (1980) is corrected.
-
-
-
0 1989 American Chemical Society
Ind. Eng. Chem. Res., Vol. 28, No. 1, 1989 109 et al., 1970). Therefore, an additional condition is needed. If CT is the concentration of carrier dissolved or initially present in the membrane before any gas is passed through it, conservation of mass requires that
LOW PRESSURE OR PERMEATE SIDE
4
'r CA
= c;
t I ! J
A.AB
t
CA
= c:
i
B
(7)
I
HIGH PRESSURE OR FEED SIDE
Figure 1. Facilitated transport membrane.
Section 5 presents numerical methods for solving the governing equations. The Galerkinlfinite element method is used to discretize the equations. The resulting set of nonlinear algebraic equations is solved by Newton's method. Newton's method leads at each iteration to a set of linear algebraic equations, the structure of which depends on the ratio DBIDU. When DB/DAB = 1, the resulting matrix is sparse and banded about the main diagonal. We describe elsewhere an efficient moving or adaptive finite element method to solve this problem ((Basaran et al., 1987) hereafter referred to as I). When DB # DAB, the structure of the matrix that results at each Newton iteration depends on which of the two formulations is used. The first formulation leads to a matrix which is sparse but banded and is augmented by a dense row and a dense column; i.e., it is arrow-shaped. The second formulation leads to a full matrix. In each case, the matrix is LU-decomposed by a different technique. Also, a moving finite element method is implemented for accurate resolution of concentration boundary layers that develop for operation near equilibrium. Section 6 presents results of numerical analysis. Computational efficiency of the two formulations for DB # D m are compared against one another and the benchmark results reported in I, where DB = D A B . Furthermore, an optimization is carried out to determine the maximum facilitation factor and the corresponding equilibrium constant for a given set of operating conditions. Concluding remarks are the subject of section 7. 2. Governing Equations and Boundary Conditions For the reaction scheme a t eq 1, conservation of mass inside the membrane requires that the concentrations of A, B, and AB throughout the film shown in Figure 1 be governed by the following three nonlinear differential equations (cf., Smith et al. (1973)): D A d2CA/dX2 = %? E kfCACB - krCAB (2)
DB d2CB/dX2 = R DAB
d2Cm/dX2 = -2
=
CA'
dcB/dx = dCAB/dX = 0
+
c* is a constant of integration. The integration constant,
c*, can be solved for solely in terms of cB by substituting eq 8 into eq 7:
Of course, eq 9 can now be substituted into eq 8 to yield an expression for CAB as a function of cB: CAB
=
CT
DB DAB
- -CB
+ (DB- DAB)/DAB S,"CB dx L
(10)
Alternatively, it can be shown that cB can be expressed as a function of cm Thus, either CB or cm is an independent unknown, but not both. Because once C A and CB are known everywhere CAB can be solved for from eq 10, only differential eq 2 and 3 need to be solved. At this point, there are two alternatives: 1. Equation 8 can be solved for CAB in terms of c* and cB and the result substituted into eq 2 and 3. One would then have to solve simultaneously eq 2 and 3 augmented by the constraint at eq 9 for C A ( X ) , c&), and c*. 2. Equation 10 can be substituted into eq 2 and 3 to yield two coupled integro-differential equations for cA(x) and cB(x). The equations for the two formulations are summarized below. But first, the governing equations and boundary conditions are cast into dimensionless form using film thickness, L, as the characteristic length and concentration of A on the high-pressure side, cAO, as characteristic concentration. When the diffusivities DB and DM are equal, the equation set is much simplified and is given by (cf., I) d2CA/df2 = P[(QFA + ~ ) E B- T ] When DB #
DAB,
(11)
formulation 1 leads to
(3) (4)
Here D stands for diffusivity, a constant, and c for concentration. x is distance measured from the high-pressure side of the membrane. Equations 2-4 are solved subject to the following boundary conditions: x = 0: CA
The global mass constraint, eq 7, is the additional condition that is needed. If eq 3 and 4 are added, their sum is integrated twice, and if use is made of the no-flux conditions at eq 5 and 6 to simplify the resulting expression, it follows that DBcB DABCAB = C* (8)
and formulation 2 leads to
(5)
x = L:
dcB/dX = dCAB/dX = 0 (6) However, if eq 2-4 are integrated across the film thickness, it can be shown that only three out of the four no-flux conditions of eq 5 and 6 are independent (Goddard CA
=
C A ~
In eq 11-17, variables with bars over them are dimensionless counterparts of those without bars introduced earlier, i.e., cAE cA/cAo, F~ E c B / c A o , F* E c*/DUcA0, and f x / L . P E k A 2 / D Aand s 1 k C 2 / D Bare Damkohler numbers that measure the relative importance of reaction
lncl Eng. Chem. Res., Vol. 28, No. 1. 1989
1J 0
rate to diffusive flux. Q = k s A o / k r is a dimensionless equilihrium constant. T = cT/cA0 is the dimensionless concentration of carrier initially dissolved in the membrane. i'i 2 DB/DAB is the ratio of the diffusivity of the carrirr tu that of the complex The choice of dimensionless g w ~ p isw d here is, of course, not unique. Instead of S aiid R. one may choose equally well the diffusivity ratios DR/.1)4= P I S and DAB/DA= P / ( S R ) . The dimensionless hniiwJnn* conditions follow from eq 5 and 6:
Substitution of the expansions of eq 23 into the nonlinear system eq 16 and 17 and boundary conditions eq 18 and 19 leads to a sequence of linear boundary value problems: d2CA(nJ/df2=
RA(n)
i2
dzCg(n'
-= RB(d = - RA(fl)
4
in (f: 0 < x < I ) , n = O.l,etc. (251
df2
-
y
5'4
=1
dCB/d.f = 0
(18)
f -. FA
=
FAL
dC,/df
=0
where 6,, is the Kronecker delta and the right-hand sides, RAtn),up to second-order, are
(19)
mhere PA', E cAL/cAo. Though formulation 1 involves an additional equation (hermw of the presence of the unknown E * ) , it is shown in section 6 that eq 13-15 can be solved numerically more efficiently than eq 16 and 17 of formulation 2. However, as far as the asymptotic analyses of sections 3 and 4 are crrncerned, using either formulation requires about the saim amourit of labor. Of course, when R = 1, ?* = T. and both ftormulations 1 and 2 reduce to eq 11 and 12 of I. Puhtracting eq 3 from eq 2 or adding eq 2 and 4 and mtegrating the resultant expression once leads t o (note: z), dl R / d"L - D A B dCAB/dX) I ) , d ( ' ~ / d X DB dCB/dX = DAdcA/dX DABdCAR/dX (20)
RA(0) =
0
(28)
S,
RA") = (QcA'O) -t- R)CR(O' -- T - ( R - 1)
1
cB'O'
df
(29)
RA(2)= (QcAi0'
s,'
+ R)c,'" + QcB(O)CA(~)( R - 1) -
cR")
df (30)
It is also assumed that S I P = DA/DBremains an O(1) quantity as P -* 0. The perturbation problem a t zerothorder is, of course, homogeneous. The perturbation problems at first-order and beyond are inhomogeneous. The condition for a solution of the inhomogeneous problem at each order n 3 1 is (see, e.g., Friedman (1956))
(11
E "o the no-flux conditions of eq 5 and 6, it follows that the facilitated flux of A is equal to the negative of the coiutant of integration in eq 20, -al. A second integration (lip
3
gi
LIRA(")da = 0
n 2 1
(31)
Calculation of the facilitation factor correct to O(P2) requires only the zeroth- and first-order corrections to the concentration of the carrier T cR(o) = ____-___ (32) 1
wlieic! a2 IS another constant of integration.
+ Q(l+ ?*L)/2
%w a useful quantity, called the facilitation factor, can i e f i n d as the ratio of the facilitated flux of A to the ils11,e cbr Fickiaa flux of A, DA(cA0 - cAL)/L:
but requires the zeroth-, first-, and second-order corrections to the concentration of the reactive gas: CA(O1
Some authors have chosen to define the facilitation factor AS 5 = F - 1 (see. e.g., Smith et al. (1973) or Kaper et al. (
= 1f
(CAI' -
1)f
(34)
IOPIi) i
3. Regular Perturbation Analysis of the
cA(2'
Equations of Facilitated Transport T wo distinct regimes of reaction/ diffusion exist in a facilitated transport membrane depending on whether the Uamkohler number P is small or large. In the limit as P 0, p q set 16 and 17 or 13-15 is amenabie to solution by regular perturbation analysis. The results reported here are an extension of the excellent analysis of Smith et al. (1973) t o the situation in which R z 1. In this section, the analysis is carried out with eq 16 and 1 e . formulation 2. The analysis with formulation 1 P\ identical results with those reported below. Yhe unknown concentrations are represented as power < P I i p s in the small parameter P
= A 1f f A f
+ A3X3 + A$ + A5f5-i-A6f6
The constant coefficients A1-A, in eq 36 are given by
+
m
m
,A
= 5, ,,-n
pllc4(fl)
pH
ZI
t:p n C A ( n ) -0
(23)
(36)
A -1 -(cAL s - 240 P
-
1)QcB'O' 5(R + Q ) +
I
Ind. Eng. Chem. Res., Vol. 28, No. 1, 1989 111 et al. (1980) have concluded that F > F1 if R > 1 and F < F1 if R < 1. It is demonstrated below that these conclusions are in error. Kaper et al. (1980) have shown that the facilitation factor F 3 1 + 3, where 3 is given by
3 = -CY us 1 A4 = --QCB(O)(CA~ 48
- 1)
1
+ Q) + QcB(0)
(49)
1011
(40)
In eq 49,
where po,pl, and p2are PO
= -dd
u2 + YoYid
+ ' / z ( ~ o+ ~ i ) d+l d) 7011
p1 = u6d - [l
From eq 22 and the foregoing results, the facilitation factor can be shown to equal
+ (d - l)c~w]7o11
P2 = d1011
(51) (52) (53)
and v0, vl, and w are defined as 1 F = 1 + -PQc,(O) 12
1 - -(PQcB(o))2 720
v0 = u
+ dro
w = (70 -
Equation 43 shows that the dimensionless group R does not enter the regular perturbation expression for the facilitation factor until the second-order problem. When R = 1, eq 43 and 36 reduce to the expressions given by Smith et al. (1973). Now the effect of the parameter R on the facilitation factor can be calculated-correct to O(P)-for the first time: (44) Letting Fl = FIR=1,
F = F1 + (aF/aR)IRS1(R - 1) + etc.
(45)
Therefore, F < F1 when R > 1and F > F1 when R < 1; i.e., the facilitation factor decreases with increasing R and increases with decreasing R. It should be emphasized that all dimensionless groups other than R, i.e., P, Q, S, T, and EAL, are held fixed in this comparison. The predictions reported in this section are confirmed in section 6 by means of finite element calculations. The integration constant F* too can be expressed in a power series in P m
= Cpnc*(n) n=O
(46)
q1 = u
d-l(ln
10
+ drl
- In 71)
(54) (55)
The dimensionless groups of Kaper et al. (1980) in eq 49-55 are related to the ones used here as d = 1/R S=P/SR u = l/QT 71 = FAL/T (56) 7 0 = 1/T Also, F* of this paper is related to CY as C*/RT = a. To calculate the effect of the ratio of the diffusivity of the carrier to that of the complex on the facilitation factor, they linearize eq 49 about d = 1:
3 = 31[1- cP(d - 1) + O(d - l)']
(57)
where 31 3 ( d = 1). They then go on to show that 9 is always positiue. However, in constructing the expansion a t eq 57, they treat 6 as a constant, which is incorrect. When a proper linearization of eq 49 is done, it is found, letting p ( ~ l + ~ -Q)/2, ~
+ Pa6 Yo Y1 +-+--2 7011 + 10 11
1011w
@,=a
(58)
where all the variables are to be evaluated at R = d = 1. By contrast, Kaper et al.'s (1980) linearization incorrectly predicts -2 in eq 58 to be -1 (see eq 43 in their paper). It can now be shown that 9 is always negative. On physical grounds, yo 2 yl. Then qo 2 q1 and 7o = v1 (1 + z ) for some z 2 0. Therefore, w = (v1z)-l In (1 + z ) and p = q1-l(1 + 1/2z)(l z)-l. From the inequality (see Kaper et al. (1980)),
+
Substituting eq 32 and 33 into eq 15 yields C*(O) = T
+ ( R - 1)cB")
(47)
it follows that w d
p.
Then, from eq 58,
9 I up 4. Singular Perturbation Analysis of the Equations of Active Transport In the limit as P m, eq set 13-15 or 16 and 17 is amenable to solution by singular perturbation analysis. Working with eq set 13-15, Kaper et al. (1980) have developed asymptotic expansions that are complete to the lowest order in P. Their expressions reduce to those of Goddard et al. (1970) when R = 1. Furthermore, Kaper
-
Yo Y1 ++10 11
2
(60)
Using the definition of p, it is easy to show from eq 60 that 9 I-up < 0 (61) Kaper et al. (1980) have also shown that the concentration profiles of A, B, and AB are given by T FA = - ( 6 [ ( ( t i Z + i))2 + F2)ll2 - (Ex + 6)] - u) (62) d
112
Ind. Eng. Chem. Res., Vol. 28, No. 1, 1989 CAB
= '{a - [Ef d
+ 6 + ((6%+ 6)' + F2)'/']1
(63) FORMULATION 1
Cg =
OIT- dC,
FORMULATION 2
(64)
where P 2 = acr/d and
5. Finite Element Analysis The computational domain 0 G f 6 1 is partitioned here into N finite elements. The elements may or may not be all of equal size. A method of adaptively sizing the elements is described below. The unknown concentrations EA and CB are expanded in this paper in a set of linear basis functions {a')(see, e.g., Strang and Fix (1973)): N+l
N+l
C cBi$'(f) (66) i=l In eq 66, cAiand cBLare the values of the concentrations a t the nodes of the finite element mesh. Equations 11 and 12, 13-15, and 16 and 17 are solved in this paper by the Galerkin/finite element method. In all three cases, the residual equations arising from the mass balance on species A-eq 11, 13, and 16-are C*(f) =
C cA'@(f)
CB(%)
=
i=l
i
= 1, ..., N+1 (67)
where C2 = (QCA f R)CB - a - b(R - l)Zc
(68)
In eq 68, (69) and a =
(
T, F*,
T,
ifR = 1 if R # 1, formulation 1 if R # 1, formulation 2
(70)
and b = 1for formulation 2 and b = 0 otherwise. Similarly, the mass balance on species B-eq 12,14, and 17-leads to the following residual equations:
i = 1, ..., N+1 (71) In formulation 2, an additional residual equation arises from eq 15: R* = F* - T - (R - 1)Z, = 0 (72) Letting c be the vector of nodal values of the concentrations of A and B, i.e., c = (cA', CB', ..., cAN+', cBN+'), and R , be the vector of residuals of the mass balances, i.e., R , (RA', R B I , ..., RAN+',RBN+'), the residual equations can be written in abbreviated form as R ( u ; p )= 0 (73) where R 3 R , when R = 1 and for formulation 2 and R (R,,R*)for formulation 1. u is the vector of unknowns; i.e., u c when R = 1 and for formulation 2 and u (c,c*)
Figure 2. Types of Jacobian matrices that arise in numerical analysis of facilitated transport.
for formulation 1. p is the vector of parameters; Le., p = (P,Q,S,T,CAL,R);R = 1 for DB = D A B and hence Is not a parameter in that case. Equation 73 is a set of N = 2(N + 1)equations when R = 1 and for formulation 2, and it is a set of N 2N + 3 equations for formulation 1. Equation set 73 is solved here by Newton's method (e.g., Isaacson and Keller (1966)), where the formula for the (k + 1)th iterate is J(u@+l)- d k ) = ) - R ( d k ) ; p ) k = 0, 1, etc. 0 4 )
J is the Jacobian matrix of partial derivatives. Its ijth entry, Jll dR'/aul, is calculated analytically. The Newton iteration was continued until the square root of the sum of the squares of both the correction vector and residuals was less than lo+. This is usually satisfied in 2-4 iterations. In all three cases considered, the Jacobian matrix is asymmetric. When R = 1, J is a sparse matrix that is banded about the main diagonal, with a bandwidth b, = 7-Figure 2a. The integrals occurring in the residuals and the Jacobian in this case are standard and are calculated by means of Gaussian quadrature (Strang and Fix, 1973). The linearized set of 2(N + 1) simultaneous equations at eq 74 is solved by means of a banded matrix solver, BANDED, using Gaussian elimination. For additional details, see I. When R # 1 with formulation 1, J is a sparse matrix that is banded about the main diagonal, again with a bandwidth b, = 7, but is augmented by a dense row and a dense column, as shown in Figure 2b. The dense column, dR,/dC*, arises because E* enters all the residual equations. The dense row, dR*/dc, is due to the integral constraint nature of residual eq 72. The integrals occurring in R and J are evaluated in the same way as when R = 1. Evaluation of R , and dR,/dc requires the same amount of work as when R = 1. In this case, however, additional work is necessary to evaluate R* and the dense row and the dense column. Evaluation of R* requires at each iteration the calculation of the integral at eq 69. Evaluation of the dense row and the dense column of J necessitates at each parameter step-see below-the calculation of the integral over the computational domain of each basis function, @Ld f . The linearized set of 2N + 3 simultaneous equations at eq 74 is solved in this case by means of a matrix solver, ARROW, that fully exploits the "arrow" structure of J (Thomas and Brown, 1987). ARROW, too, performs an L u decomposition using Gaussian elimination. When R # 1 with formulation 2, J is a full matrix, as shown in Figure 2c. This is because the integral, I,, inside the outer integral in the residuals and the Jacobian leads to a global coupling of the unknowns. The entries of R and J are calculated numerically as described previously. However, these evaluations require more computational work in this case than with formulation 1. The two computationally taxing integrals are
soi
Ind. Eng. Chem. Res., Vol. 28, No. 1, 1989 113
I , = J1@i(?)Ic dx
(75)
As before, I , is evaluated once a t the beginning of each Newton iteration. The basis function integrals in eq 76 are evaluated once at each parameter step and stored. The linearized set of equations at eq 74 is solved by means of a full matrix solver, FULL, employing Gaussian elimination. Critical to Newton's method in solving the nonlinear eq set 73 is a good initial estimate u(O)which must be accurate enough to fall within the domain of convergence of the method. First initial estimates are provided by the asymptotic or perturbation solutions described in sections 3 and 4. Thereafter, a converged solution, u ( p , ) , is available. Steady solutions to the governing equations are obtained by changing sequentially the parameters. Initial estimates for solutions a t parameter set p = p , + Ap are obtained by first-order continuation over the increment Ap (Kubicek, 1976; Brown et al., 1980). The initial estimate is approximated as
The sensitivity derivative or the continuation vector du lap,,,, where p mis the mth parameter, is calculated by solving the linear system au dR J-=-(78) aPm aPm where J is the converged Jacobian of the previous solution. Computationally, eq 78 is solved simultaneously with eq 74 in the last Newton iteration. The updated or new finite element mesh at the parameter set p , A p , is calculated by equidistributing over the elements in the new mesh the linear interpolant of the square root of a constrained integral estimate of the local spatial discretization error in the explicitly predicted profile u(O)(p0 A p ) . The procedure that is used is a variant of a moving finite element method developed by Benner ((1983) see, also, Benner et al. (1987)) and Heiba (1985). This is described in detail in I. Once the updated finite element mesh for the new parameter step is computed, the continuation-predicted solution is simply linearly interpolated onto the new mesh and becomes the initial estimate to the first Newton iteration. The nodes remain fixed during the Newton iterations. Thus, the basis function integrals in eq 76 need only be evaluated once, before the first Newton iteration at that parameter step. Kemena et al. (1983) have shown that when P, T, S, and EAL are held fixed and R = 1, the facilitation factor attains a maximum value, Fopt, when the dimensionless equilibrium constant reaches a critical value, Q, In this paper, the optimum facilitation factor, Fopt, and &e corresponding dimensionless equilibrium constant, Qopt,are determined when R # 1. To do this, the condition for the occurrence of the optimum is expressed as an additional residual equation, R, (79) Ropt = dF/dQ = 0
73 is found a t two nearby points Q1, Qz = Q1 + t, by Newton's method, as shown in eq 74, wi = (ui,QJ w2 ( ~ 2 9 Q 2 ) (80) At each point, dFldQ is calculated exactly by means of eq 22, via
The derivatives dEB(l)/dQ and dEB(0)/dQ are available from the Fist-order continuation calculation at eq 78. The location of the optimum is then estimated as
where (dFldQ), and (dF/dQ)2are eq 81 evaluated a t Q1 and Qz. Steps 80-82 are repeated until the change in Qopt between two successive iterations is less than a prescribed tolerence. By way of contrast, Kemena et al. (1983) did not use first-order continuation and consequently had to compute three (instead of two) solutions, wl, w2, and w3, to estimate dFldQ at two points, Q1 and Q2. Therefore, the algorithm used in this paper is more efficient that the one that they used. Following Kemena et al. (1983), optima in Q were tracked as S was varied, holding T , P I S , EAL, and R fixed. Once the optimum equilibrium constant is known for two values of the Damkohler number S, S1 and S2, first approximations ( u 3(0),Qopt,3(0)) for calculating the optimum value of Q, Qopt,S,a t S = S3 are constructed by the finite difference analogue of the continuation method (Ungar and Brown, 1982):
+
+
Equation 79 is solved simultaneously with eq 73 for u and €J at the optimum. Because the sensitivity of Roptwith respect to u and Q cannot be calculated analytically (as the Jacobian matrix at eq 74), the augmented equation set is solved by a chord technique. To do this, solution of eq
(84)
The algorithms were programmed in FORTRAN, and calculations were carried out on the Vax 111785 and the IBM 3090 at Air Products and Chemicals, Inc. Execution times on the IBM 3090 were roughly 22 times faster than those on the Vax 111785. Computation times reported in the next section are measured in Vax 111785 central processor seconds. The correctness of the residuals was attested by comparison of the predictions of finite element analysis with those of the asymptotic solutions where the latter are valid, i.e., small and large P. The correctness of the Jacobian was determined by confirming that the convergence of Newton's method attained its asymptotic quadratic rate as the iteration proceeded (Isacson and Keller, '1966). 6. Results The execution times for single Newton iterations using BANDED, ARROW, and FULL are shown in Table I. The execution time is the sum of (1)the time spent assembling the finite element equations, i.e., the time required to form J and R ,and (2) the time spent solving the resulting system of linear equations, i.e., the time required to solve eq 74. Table I makes it plain that as the number of unknowns, N , becomes large, the computational cost of formulation 2 becomes prohibitive. This rapid increase in computation time with N for formulation 2 is associated with the cost of solving a full system of algebraic equations,
114 Ind. Eng. Chem. Res., Vol. 28, No. 1, 1989 Table 1. Execution Times for Single Newton Iterations time/iteration (Vax 11/785 cp s) ( % total time) 20 elements 40 elements
form J, R solve Jb = -R total time
BANDED
ARROW
FULL
BANDED
ARROW
FULL
0.05 (71.4%) 0.02 (28.670) 0.07
0.06 (66.7%) 0.03 (33.3%) 0.09
0.15 (35.7%) 0.27 (64.3%) 0.42
0.09 (69.2%) 0.04 (31.8%)
0.11 (68.8Yo) 0.05 (31.2%) 0.16
0.45 (18.8%)
0.13
1.94 (81.2%) 2.39 I
”.. I
r
a
to
20
so NUMBER OF ELEMENTS
40
I
sa
Figure 3. Execution time to reach a converged solution.
which requires O(N 3, operations (Strang, 1976). The increase in computation time with N for formulation 1and for the situation when R = 1 is modest because solving the system of algebraic equations in these two cases requires O(NNbW2)operations (Strang, 1976; Thomas and Brown, 1987). The total execution time to reach a converged solution as a function of the number of elements is shown in Figure 3. The execution times shown are average times for several calculations, each solution requiring four Newton iterations for convergence. The execution times shown in Figure 3 include CPU time spent during program initialization, generation of initial estimates to Newton’s iteration, etc., in addition to the times shown in Table I. With formulation 1 (ARROW), it takes 11-17% more computer time to analyze the performance of a membrane in which R # 1 than a membrane in which R = 1. In contrast, with formulation 2 (FULL), it takes 71-146570 more computer time to analyze the performance of a membrane in which R # 1 than a membrane in which R = 1. Of course, calculations with both formulations 1 and 2 give the same numerical results. Usually, 30-40 moving elements are required for accurately analyzing the performance of membranes operating near equilibrium. This fact renders formulation 2 unattractive to use in engineering evaluations. Jain and Schultz (1982) solved the equations of formulation 1 by the method of successive substitutions or Picard iteration. That is, their algorithm was (i) solve eq 13 and 14 holding e* fixed, (ii) update c* from eq 15, and (iii) if the change in F* between two successive iterations is less than a prescribed tolerance, stop; otherwise, go back to (i). Because several iterations (i)-(iii) are necessary before e* converges, the computation time with their algorithm would be several times that of BANDED even if their method of solution at step (i) was as efficient as BANDED. However, we show in I that BANDED is at least an order of magnitude more efficient than their algorithm for step (i). Moreover, it is now well-known that Picard-type iterative schemes can lead to divergent iteration in certain ranges of the parameters (see, e.g., Saito and Scriven (1981) or Ettouney and Brown (1983));the full Newton method used in this paper where all the equations are solved simultaneously in the same iteration is free of convergence problems.
h
P=o 1
0
PZ05
94
0 92
os
0.7
09
11
I S
I5
17
19
CARRIER/COMPLEX DIFNSIVITY, R
Figure 4. Facilitation factor versus carrier/complex diffusivity (near diffusion limit).
CARRIER/COMPLEX
DIFFUSIVITY, R
Figure 5. Facilitation factor versus carrier/complex diffusivity (near equilibrium limit).
The variation of the facilitation factor with the ratio of carrier to complex diffusivity for situations in which the importance of diffusive flux is large compared to reaction rate is shown in Figure 4. The results show that F falls as R rises, in accordance with the prediction of eq 44. When P = 0.1, the predictions of finite element calculations and perturbation results of section 3 are virtually identical. However, when P = 0.5 and 1,the perturbation analysis predicts that the facilitation factor is less than one for some values of R. This is physically not possible and signals the breakdown of the asymptotic expansion at eq 43 for these values of P. Nevertheless, the asymptotic analysis is valuable because its predictions are qualitatively correct and it offers a means for checking the accuracy of the numerical calculations when P is suitably restricted; i.e., P 1 or R < 1. This intuitive observation is consistent with eq 43, which is valid for P 0. However, whether it is consistent with eq 49, which is valid for P ---* m , is not apparent. After some algebraic manipulations, it can be shown that as Q 0. CY 1, and eq 49 reduces to --1
+
I
04
1
01
I
I
1
10
100
1000
+
lOOGO
DAMKOHLER NUMBER. S
Figure 13. Optimum facilitation factor versus Damkohler number.
-
Thus, as P m and for R = 1, Qoptapproaches m, 4.47, and 2.24 when EAL = 0, 0.05, and 0.2. The finite element computations plotted in Figures 9-11 are in agreement with the asymptotic results at eq 85 and 86. The analysis of this paper shows for the first time that Qoptis a strong function of R for operation near equilibrium and that Qopt increases with R. How the optimum facilitation factor varies with the Damkohler number S is shown in Figures 12-14. From eq 43 and 85, lim Fop,= 1 +
P-Q
PT
6 (1 + EAL)
I- O ( P ?
(87)
When EAL = 0,0.05, and 0.2, eq 87 predicts that Fop, = 1.83, 1.79, and 1.69. The values of Fop,predicted by finite element analysis for these three values of EAL at S = 0.1 are Fop,= 1.62, 1.59, and 1.52. Thus, the value of the optimal
-
Thus, the diffusivity ratio ( R ) also has a big effect on facilitation when P m. The effect of R on the concentration profiles of A, B, and AB and the values of CY and 3 are shown in Figure 15 for situations in which Q = and P m. The concentration profiles and the values of a and 3 shown in Figure 15 and also in Figure 16 were obtained from eq 62-65 and 49 and 50. Figure 15 makes it plain that the profile for CAB is essentially unaffected by R , and F,therefore, responds in accordance with intuition. It is surprising that R seems to have no effect on facilitation at large Q. This observation is consistent with eq 43, which is valid for P 0. Once more, whether it is is not consistent with eq 49, which is valid for P m, immediately apparent. It can be shown that as Q CY d , and eq 49 reduces to
--
-
- 05,
--+
(90)
Ind. Eng. Chem. Res., Vol. 28, NO. 1, 1989 117 (1)IT: 0 99980 cg
+
L
c U
c U
z
t
0
c
w
a :
io00
i-0497
? J
I
1c-2
E,( r
t
chg ( 1 ) I T * 1 9 9 7 9 x 1 0- d
0
02
04
06
08
10
DISTANCE, A
rB(0)/T 9 9897x10
....
1
0.2
0
0,4
0.6
08
C,,(OllT = 1.0003
I
I CAEIT
~
~:IOOO.
1
i-00500J
F
0
I
-1
10-2
c,(0)/T = 1 0 0 0 3 10 - 3
.
0.9924 1
I
i t
1: 10-4
10
x
-
1
i
10-3
DISTANCE 0)iT:
1) I T
4 9 7 5 1x 10-3
10-3 0
02
0 4
0 6
DISTANCE,
08
10
k
t z
t
i
lo-'
U 0
.
1-
t r
E
i
C r
c4, 9(1)IT 8 1 6= ~ 1 0 - ~
10
1 0
--_L0 2
-
10 4
-
0 6
08
10
DISTANCE x
Figure 15. Effect of R on concentration profiles and facilitation factor in the limit in which the reverse reaction dominates the forP m , PIS = lo-*, T = lo3, cAL = 0.2. ward reaction. Q =
-
a result which is independent of R. Nevertheless, a high Q implies a low cB/cAB-Figure 16-and it is surprising that the facilitation is not improved if DB/DAB> 1 or R > 1. Intuition and observation can be reconciled by examining the effect of R on the concentration profiles of the diffusing species for Q a. Figure 16 depicts situations in which Q = lo3 and P a. At high Q, the values of FB are very low and Figure 16 shows that the profile for cB is little affected by R. Thus, DB dcB/dx is little affected
-
-+
c g (O)/T = 9 9 8 3 8 ~ 1 0 . ~ 10
I 0
04
02
06
08
10
DISTANCE, x
Figure 16. Effect of R on concentration profiles and facilitation factor in the limit in which the forward reaction dominates the rem, PIS = T = io3, = 0.2. verse reaction. Q = io3, P
-
by R. Hence, by eq 8, DAB dcm/dx is also little affected by R. Therefore, F is nearly constant. These arguments can be quantified by careful scrutiny of Figure 16. Figure 16 shows that cB(l)is lowered and cB(0) is increased when DB/DAB= 2 relative to when D B / D m = 1, which is to be expected. Because Q is high, this causes the concentration
118 Ind. Eng. Chem. Res., Vol. 28, No. 1, 1989
of the complex on the high-pressure side of the membrane, F m ( 0 ) , to be larger when R = 2 than when R = 1. Hence, a larger fraction of A is being transported across the membrane in the bound state, as AB, when R = 2 than when R = 1. This hurts the facilitation in the former case compared to that in the latter case because Dm is lowered when R = 2 relative to when R = 1. The extent of the lowering can be quantified by careful examination of the concentration profiles for A and AB in Figure 16. If flux is measured in units of D A C A o / L , the dimensionless fluxes of A and AB-at any point in the membrane are N A EdCA/dZ and N A= ~ -(l/R)(P/s)dc,/df. When R = l, N ~ ( 0 . 2= ) 1.011 and "(0.2) = 0.190; when R = 2, N ~ ( 0 . 2 ) = 1,010 and NAB(0.2) = 0.185. Here the location f = 0.2 in the membrane was chosen solely for the purpose of illustration. Also, using eq 22 it can be shown that the facilitated contribution to the flux decreases from 0.497 when R = 1 to 0.495 when R = 2. In other words, for both Q 0 and Q m, the constituent present at low concentration levels will be virtually unaffected by R. At high Q, this nonetheless can change the gradient of c, in order to satisfy the condition that the fluxes of B and AB be equal and opposite.
-
-
7. Discussion and Concluding Remarks According to the foregoing results, the effect of the ratio of the diffusivity of the carrier to that of the complex, R, on transport through a membrane becomes increasingly pronounced as the ratio of the reaction rate to diffusive flux increases provided that the rate of the forward reaction is not too large compared to that of the reverse reaction. R can differ greatly from one, especially if the carrier and the reactive gas molecules are of similar size. Under most circumstances, R > 1. When R > 1,the actual facilitation factor would be less than that predicted by a simplified model that assumes that R = 1. The analysis and results reported in this paper properly address this situation. What is depicted in Figure 1 is just a small section or a differential area of an actual membrane module, called an FTM (Facilitated Transport Membrane) based module in I. To determine the performance of an FTM-based module, eq 13-15 subject to the boundary conditions 18 and 19 have to be solved repeatedly for each differential area in turn. This was done in I for situations in which R = 1. The differential area program ARROW developed here has been integrated into the MODULE program of I to analyze the performance of commercial-scale FTM-based separations in which R may not equal one. Moreover, the analyses reported in this paper for planar membranes (a flat-sheet or a spiral-wound membrane configuration) have been extended to cylindrical (a hollow-fiber membrane configuration) and spherical membranes. The governing eq 2-4 are readily generalized to any one of these three geometries by replacing x by r and d2/dr2by ( l l r j ) dldr ( r j dldr). Here, r is the radial distance measured from the geometric center of the cylindrical or the spherical membrane and j = 0 for a plane, j = 1 for a cylinder, and j = 2 for a sphere. Boundary conditions 5 and 6 are applied at r = ri and r = ro,where ri and ro stand for the inner and outer radii of the membrane. An interesting outcome is that, for operation in the near diffusion limit, asymptotic analysis shows that the dimensionless group R does not enter the perturbation expansions for the facilitation factor until O ( P 2 )in any geometry. Furthermore, the Galerkinlfinite element method that is used in this paper is ideally suited to solve the governing equations in any of the three geometries. The effort involved to do so amounts to no more than
having to vary the index j between 0 and 2.
Acknowledgment We thank (1)Dr. P. L. T. Brian of Air Products and Chemicals, Inc., for his critical review of the paper, (2) Prof. R. A. Brown of MIT for providing us with a preprint of his paper and the computer code of the ARROW solver on diskette, and (3) one of the reviewers (i) for encouraging us to examine the effect of R on facilitation in the limits as Q 0 and m and (ii) for offering a number of valuable comments that helped to improve the presentation of the materials of section 6.
-
Nomenclature A, B, AB = reactive gas, carrier, complex Al-A6 = constants of integration defined at eq 36, dimensionless a, b, = coefficients in eq 68, dimensionless 6, b, E2 = coefficients defined in eq 62-65, dimensionless -al = facilitated flux of A (see eq 20), mol m-2 s-l u2 = constant of integration defined in eq 21, mol m-l s-l b, = bandwidth, dimensionless cA,cg, CAB = concentrations of A, B, AB, mol m-3 cA, cg = cA/cAo, cg/cAo, concentrations of A, B, dimensionless cT = concentration of carrier initially dissolved in membrane, mol m-3 cA0= concentration of A on high-pressure side of membrane, mol m-3 cAL= concentration of A on low-pressure side of membrane, mol m-3 cAL= cAL/cAo = concentration of A on low-pressure side of membrane, dimensionless c* = constant of integration defined in eq 8, mol m-l s-l e* = C*/DABCAo, constant of integration, dimensionless cA("), ~ g ( " )c*(") , = nth correction to concentrations of A, B, and c* in regular perturbation expansions in eq 23, dimensionless , . cA(,cgt = nodal values of concentrations, dimensionless DA, Dg, DAB = diffusivities of A, B, AB, m2 s-l d = 1/R, dimensionless F = facilitation factor defined in eq 22, dimensionless 3 = F - 1, alternative definition of facilitation factor, dimensionless F, = facilitation factor F evaluated at R = 1 Y1 = facilitation factor F evaluated at d = 1 Fop,= optimal facilitation factor, dimensionless I , = integral of concentration of B across the membrane, dimensionless J = Jacobian matrix, dimensionless kf = forward rate constant, m3 mol-l s-l k , = reverse rate constant, L = membrane thickness, m N = number of finite elements N = number of unknowns N A ,NAB= fluxes of A and AB, dimensionless P = kJ2/DA, Damkohler number, dimensionless p = (P,Q,S,T,cAL,R), vector of parameters, dimensionless po, p l , p 2 = terms in quadratic equation for a (see eq 51-53), dimensionless Q = k$Ao/k,, equilibrium constant, dimensionless Qopt = optimal equilibrium constant, dimensionless R = Dg/DAg, diffusivity ratio, dimensionless 2 = generation term in eq 2, mol m-3 s-l RAi,Rgi,R* = residuals, dimensionless R A (n), RB(") = right-hand sides in eq 24 and 25, dimensionless R , , R = residual vectors, dimensionless S = k&'/DB, Damkohler number, dimensionless T = cT/cAo, concentration of carrier initially dissolved in membrane, dimensionless
Znd. Eng. Chem. Res. 1989,28, 119-121 u , c = vectors of unknowns, dimensionless u (O) = initial estimate t o Newton's method, dimensionless w = ( u ,Q), solution vector augmented by Q, dimensionless x = distance, m f
= x/L, distance, dimensionless
Greek L e t t e r s a = t * / R T , dimensionless yo, y1 = 1/T, tAL/T, dimensionless
S = P/SR, dimensionless t = small number, i.e., c
