Equilibrium Constant Estimation and Model Distinguishability

Equilibrium Constant Estimation and Model Distinguishability Gary E. Blau,* Richard R. Klimpel, and Edwin C. Steiner The Dow Chemical Company, Computation Research Laboratory, and E. C. Britton Research Laboratory, Midland, Mich. &?S,$O

Both in estimating parameters and in discriminating among rival mathematical models of a chemical reaction system, the presence of experimental errors makes the generation of physically meaningful results difficult. It i s particularly important to know if the complexity of a model i s warranted by the data. Techniques are developed for estimating equilibrium constants and evaluating different models for the general chemical reacting system model BiAj-i A1 M; j = 1, ., N). Then these techniques are BiAj ( i = 1, used to determine the preferred state of aggregation of pure metal alkoxides and the effect of added alcohol on these aggregates, using experimental data obtained by ebullioscopic techniques.

Downloaded by UNIV OF GEORGIA on September 5, 2015 | http://pubs.acs.org Publication Date: August 1, 1972 | doi: 10.1021/i160043a007

+

Considerable effort has recently been devoted to developing experimental design techniques for estimating parameters (Draper and Hunter, 1966; Kittrell, et al., 1966), distinguishing rival models (Box, 1969;Box and Hill, 1967), or some combination thereof (Hill and Hunter, 1966), in nonlinear systems. These methods sdggest that a sequential experimentation-analysis scheme be-followed such-that the iocation of each new data point be governed by all previous results and the specific design criterion under consideration. Unfortunately, these methods cannot be used on problems i n which a large number of parameters and potential models are involved because of the excessively high demands on computer time and memory. This paper deals with such a system.

...,

..

where B1 and -4,are arbitrary chemical species and BiAj is a species comprised of i molecules of B and j molecules of A. The chemical reaction for the formation of BiAl is BiAj-I

+ -41If Bpi,

(1)

with associated equilibrium constant K i j given by

The development is sufficiently general that it could be extended to other systems, for example a three-component system in which

General Chemical Equilibrium Model

I n a recent paper (Blau, et al., 1972), the problem of parameter estimation in physicochemical models of equilibrium systems was formulated on a statistically sound basis using the maximum likelihood principle (Hogg and Craig, 1965). Statistical criteria were developed to distinguish modelling error from experimental error so that the best model consistent with the data could be selected. The method is here applied to estimate equilibrium constants and to distinguish models for the general chemical equilibrium system

324 Ind. Eng. Chem. Fundam., Vol. 1 1 , No. 3, 1972

It will be useful to introduce some notation to describe a specific model within the general one; that is, to define its size and shape. Unless stated otherwise, only “rectangular” models are examined so that by giving values to JI and N the model is uniquely described. Although the path by which any species is formed will not affect its equilibrium concentration, it is convenient for computational purposes to define such

a path. I t is arbitrarily assumed that, e.g., species BzA2 is formed by the sequence of reactions K 20

+ Bi IfBz Bz + -41 IfBzAi Bi

solved by the approach of Brinkley (1947) which is to exploit the structural features of eq 2 and express the equilibrium concentrations in terms of the monomer concentrations (B)k, (A)k and the K,, i

K?l

(B2)k

=

J-J

(i = 2, . . ., M)

K,o(B),,"

n=2

(6)

Kzz

BzAi

+ Ai _J B2Az

and

The notation for such a path is P(B2Az) : (B1, Bz, BzA, BzAz). A specific model can then be denoted simply (M, N ) , where the paths

, Bi, BiA, Bt-42,


P(Bi, A N ) :(Bij Bz,

BZAN) (i = 1,

.

,M)

are understood. For any general equilibrium system the problem is to select the best or "optimal" model ( M * , N * ) from experimental measurements. The usual situation is where the data have already been collected and are the only information available to identify this model. Criteria for determining the most suitable location for additional experiments, if necessary, exist but have not yet been applied to models with the degree of complexity of system I. Blau, et al. (1972), have shown that experimental data which can be described by the general chemical equilibrium model consists of (1) formal concentrations A0 and Bo representing the amount of material added to the system initially and (2) some physicochemical property measurement yo which can be related to the overall concentration of the various species present a t equilibrium. In actual practice, yo has taken the form of such measurements as boiling point elevation, total conductivity, and radiation levels. If the kth data point is represented by B BO)^, ( A o ) ~YO)^) , and the equilibrium concentration of the species BiAj is given by (Bt.A,)k, then for some model ( X , N ) material balance equations relating the formal concentrations to the equilibrium species concentrations are given by

n i

=

j

Kno

n=1

J-J Kirn(B)k'(A)n'

m=l

(i = I, . . ,, M ; j

(&)k

=

( N k

+c

N

c j(B*Aj)k

a=1 j = l

(4)

(Bo)k=giv,M(Kij, (B)k, ( A ) J

(A0)k = hN,M(KI,,

= fN,M(K*j,

(BiAJk)

(5)

where the form of f N , M is determined by the specific system and measurement involved. Mathematical Analysis

Analysis of the data in terms of model I involves three different but related problems: (1) given a set of equilibrium constants K i J , to compute the equilibrium concentrations of the various species; (2) given a specific model of the form I, to determine a set of equilibrium constants Ki,*which best correspond to the experimental data; and (3) to select the model most suitable for representing the data. (1) Chemical Equilibrium Problem. The determination of the equilibrium concentrations (BtAJk for the kth data point given specific values of the equilibrium constants K,, is called the chemical equilibrium problem (Brinkley 1947; White, et al., 1958). The problem may be

(7)

(A),

(B)x

=

+

M

n i

i i=2

Kno(B)ki

n=2

+

(B)k, @)!J= M

N

+ c c j ri Kno fI n=l

i=lj=1

m=l

Kirn(B)ki(4kj

(9)

Solving this system of two equations for the two unknowns, (B)k and ( A ) k , allows the equilibrium concentrations (BiAj), to be obtained from eq 6 and 7. Shapiro and Shapley (1965) have shown that for any positive set of equilibrium constants, K i j > 0, and positive formal concentrations (BO), > 0 and (Ao)k> 0, the solution to this nonlinear equation system is unique. It may be obtained by an extension of the self-generation procedure described by Blau, et a2. (1970), or by some other nonlinear equation solving algorithm of the Newton Raphson type. (2) Parameter Estimation. The second problem is t o determine the best or optimal estimates of K , describing T sets of data ((Bo)kl( A o ) ~ (, ~ 0 ) ~ ) ( k = I, . . ., T ) . For any data point k, let the measured value of aggregate concentration (yo)k, be related to the true value uk by ulc

+(dk

(10)

where (ey)* is the error in measurement of yo)^. Similarly, let the formal concentration measurements be given by

(BO)*= (PI,

The physicochemical property measurement (yo)k is related to the equilibrium concentrations by the function (YOL

1, . . ., N )

Substituting eq 6 and 7 into the material balance eq 3 and 4 gives

(YO)k =

M

=

(Adx

=

(a)x

+ (€elk

+

(11)

(12)

(EA)k

where (p)k, are true values and ( E A ) I , are experimental errors. For any specific model ( X , Y) the observed data are related to their calculated values and experimental error by (yo), = f N , M ( K * j ,

(BO)k

W k ,

= gN,M(KijJ

(A),)

+

(€M)K

(B)k, (A),)

( ~ 4 0 )= ~ h N , M ( K < j J( B ) k ,

+

+ (ev)a

('B)k

(A),) f (€A),

(13)

(14) (15)

where f N , M and gn,m,and hn,rnwere defined by eq 5, 8, and 9, respectively, and ( E . ~ ) represents ~ the modelling error. Equations 14 and 15 describe a nondeterministic version of the chemical equilibrium problem. For any set K i j the best estimates of (B)k and (A)k are readily obtained by solving these equations simultaneously with ( E B ) ~= ( E A ) , = 0. In mathematical programming terms, these optimal estimates ((B)k,(A)k)are said to be feasible points. Ind. Eng. Chem. Fundam., Vol. 1 1 , No. 3, 1972

325

The estimation of appropriate values for the equilibrium constants may be accomplished by the use of the maximum likelihood principle (Bard and Lapidus, 1968; Hogg and Craig, 1965; Reilly, 1970). The approach is restrictive in that it requires explicit assumptions concerning the form of the probability distributions of the errors. However, once the form has been assumed, any unknown parameters appearing in the probability distribution can be generated along with the unknown parameters in the model. This implies that two mathematical models are being examined: one describing the chemical reaction system, and the other, a probability distribution for the errors. If we restrict our attention to feasible N ) with parameter estimates points only, for any model (M, K i j ,residuals corresponding to these points become

(YO), - f N , M ( K i l * , (B)k, (A),)

=

(eu)k

(16)


Upon substitution of these residuals for the experimental errors in the joint probability density function of all the errors in the observed values (yo)&,the maximum likelihood approach yields the following weighted least-squares minimization problem when the experimental errors are known

subject to the nonnegativity conditions (B), 2 0

(4, 2 0 Ki, 2 0

(18)

Equation 18 restricts the concentrations and equilibrium constants to those values which are physically meaningful. The following assumptions are necessary to arrive a t this mathematical program. (i) Both (B)&and ( A ) k used in calculating the residuals are solutions to the chemical equilibrium problem for the specific Ki,under consideration. That is, all points examined by eq 15 are feasible. (ii) The experimental errors from point to point are uncorrelated. (iii) The experimental errors for each data set are independent, are normally distributed with zero means and identical covariance matrices, and are known a priori. That is E ( e J k 2 = ru2 and E[(€,),. (eu)j] = 0 ( k = 1, . . ., T ; j = 1, . . .,T ; k # j). The development of eq 17 by the maximum likelihood approach is available in the literature (Bard and Lapidus, 1968; Blau, et al., 1972). The consequences of relaxing assumptions (ii) and (iii) are discussed by Bard and Lapidus (1968). An extensive designed experimental program is necessary to determine lack of independence. In many cases this is either impossible (large between sample variation) or too expensive (excessive experimental time). Estimation techniques are available but it will suffice to invoke the independence assumption suggested by exploratory experiments on the system. The validity of the normality assumption will be discussed in a later section. The optimal parameter estimates Ki,*are those values which minimize the “objective function” S M , Nsubject to the “univariate constraint” eq 18 and the feasibility requirement. Several methods are available, a t least theoretically, for solving minimization problems of this type (Wilde and Beightler, 1967; Zoutendijk, 1966). The nonnegativity constraint coupled with the feasibility requirement render most of these methods impractical. Blau, et‘al. (1970), have shown how this constrained minimization problem can be transformed into an unconstrained problem, albeit of larger dimension, which lends itself readily to solution by unconstrained techniques such as the conjugate gradient method (Davidon, 1959; Fletcher and Powell, 1963). Starting from some initial estimates Kilo, successive param326

Ind. Eng. Chem. Fundam., Vol. 1 1, No. 3, 1972

eter estimates K i t , K J , . . . , K i t are generated by the minimization algorithm until a minimum value for S M , Nis obtained corresponding to the optimal parameter estimates K i j * = K i t . That is

SM,N(Kjj(’))-

5 60

SM,N(Kij(’-’))

(19)

and

] K i p - K i p ’ )I

< 6a3 (i

= 1,

. . ., M;j

=

1,

where the 6ii are some arbitrarily small positive real numbers. In general there is no guarantee that the minimum defined by eq 19 and 20 is unique. It is frequently necessary, therefore, either to examine the results with regard to the physical system involved or to perform additional experiments to verify the parameter values obtained (3) Model Selection. Having developed a method for generating optimal parameter estimates Kij* for a given model ( M , N ) , it remains to find a procedure for selecting the best model (M*, N * ) from the general system I. An inherent assumption is that such a model exists. That is, the general model contains elements defining a specific model which adequately describes the experimental data. The problem, therefore, is to discriminate rival models from the general chemical equilibrium system. KO attempt will be made to distinguish different general model systems. It is often relatively easy from physical considerations to determine the minimum model size which might explain the dat’a. However, determining how large the model needs to be is a nontrivial problem. Let the minimal and the best model sizes to be denoted ( X O , NO) and (,If*, N * ) , respectively. By definition the best model will be the smallest model such that the data is explained within experimental error a t a given confidence level. Given T sets of experimental data (Bo),, (A,Jk) (k = 1, . . ., T), the problem is to develop a procedure for generating ( M * , N * ) from NO) and obtain the equilibrium concentration distribution by determining Kij* for all reactions along the paths P(BtAN*): (Bi, BP, . . ., Bi, BiA, BiAz,

(i = 1, . . ., M * ) A procedure similar to stepwise multilinear regression (Draper and Smith, 1966) may be used to distinguish the best model. The stepwise add analog would be executed by starting with the minimal model ( X O , NO) and adding species to form a sequence of progressively larger models. For any model (M, N ) , optimal parameter estimates Kij* are found by solving the minimization problem described in the previous section. Blau, et al. (1972), have shown that a t the minimum S M , N ( K j j *it) is possible to test for the adequacy of the particular model under consideration by examining the mean square ratio T

RM,N*=

C k=l

(e,*)k2/(T

-P

-

1)

(21) SY

+

where P = X ( N 1) - 1 is the number of equilibrium constants, Zk_lT(e,*)r2 is the sum of squared residuals evaluated a t Kij* and sy2 is an estimate of the experimental error in yo generally available from independent measurements on the system. Two conditions can occur. The first condition is

RM,N*2

Fn,,n, (a) =

5

(22)

-

where F is the probability distribution with nl = T P - 1 the number of degrees of freedom for the numerator and n~ the degrees of freedom for the denominator, and a is the probability that the ratio is greater than the tabulated F value. Under this condition there is a significant difference between modelling error and experimental error. This implies that the current model is inadequate or, in statistical terms, the null hypothesis of no significant model error is rejected. The second condition is


RM,N*

< Fn,,,,

(a) = 5

(23)

I n this case it is impossible to distinguish modelling error from experimental error or alternatively, the null hypothesis of no significant model error cannot be rejected. Therefore, the present model is large enough. An examination of still larger models is physically meaningless with the present level of experimental error a t the specified confidence level. By definition, the smallest model for which R M , N * is approximately equal to or slightly less than 5 is the best model (M*, N * ) . The parameter estimates Kij* and corresponding feasible points ((B),*, (A),*) can then be used to calculate the equilibrium concentrations (BiAj)k*. The stepwise delete analog could be executed by starting with a model obviously more sophisticated than warranted, for example, taking P > T and dropping terms until the smallest adequate model is found. This method is computationally inefficient since larger models are being used throughout the procedure. The model adequacy test condition, of course, depends on the confidence level assigned in the F test. Thus if a is taken as 0.05 and R M , N * > 5, the probability is 95% that the model error is greater than experimental error and a larger model is justified. On the other hand, if R M , N * < 5, the physical system may contain species of a larger model but the probability of their existence is less than 95% so that invoking a larger model is doubtful. The condition also depends on the validity of the assumption of normal error distributions. As was discussed in a previous paper (Blau, et al., 1972), this assumption generally becomes valid as the modelling error approaches experimental error, Le., at the time the F test is actually applied. This will be demonstrated in the discussion of computational results.

combining intuitive physical chemistry argumentis with an analysis of the experimental error, no attempt being made to generalize or quantify this concept. It was discovered that the tetramer is the energetically preferred state of aggregation but that aggregate species up to the octamer must be assumed to describe the data properly. An experimental investigation has now been made to determine the effect of added alcohol on the state of aggregation of the alkoxide. Ebullioscopic techniques were used to determine the aggregate concentration at equilibrium. This physicochemical property is related to the equilibrium composition and to the average aggregate size N Oin the following way

(Ad,

M

N

i=l

j= 1

+ c ( ( B A + c (BiAj),)

(24)

where BI represents an alkoxide moiety, AI an alcohol moiety and BIAj an aggregate comprised of i molecules of B and j molecules of A. Solving for each of the aggregates in terms of the equilibrium constants and the concentrations of the alkoxide and alcohol moieties by eq 6 and 7 gives

Data were collected over the concentration range 0.0026 5 (Bo) 0

5

5 0.3

m

(26)

(Ao) 5 0 . 5 m

The upper and lower bounds result from experimental limitations. A sequential experimental procedure was followed in a manner designed to concentrate data collection where the state of aggregation changed the most rapidly. Figure 1 depicts the location of the data points and various experimental “runs” on a graph whose coordinates are the alkoxide and alcohol concentrations. Some of the 165 data points are

Application

Chemical Problem Description. Alkali metal alkoxides are used extensively in synthesis and in theoretical studies of chemical reactions, yet relatively little is known about them as a class of compounds. An investigation was undertaken t o learn about the nature of alkoxides in media of relatively low dielectric constant and particularly to learn the effect of added alcohol on the solutions. Polyether alcohols and their alkoxides were chosen for the study because of their high solubility in organic solvents, the low volatility of these alcohols, and the commercial interest in them as the reactive species in alkylene oxide polymerizations. One important phase of this investigation was the ebullioscopic determination of the average state of aggregation of pure alkoxides as a function of concentration in diglyme solvent. Blau, et al. (1970), were able to analyze the data in terms of successive association equilibria and to obtain estimates of the equilibrium constants. Implicit in the process was the need to identify the minimum number of aggregate sizes which were necessary to explain the data satisfactorily. The distinction between different models was achieved by

0.001

0.0 1

0.10

Bo ( ALKOXIOE CONCENTRATION I MOLAL

Figure 1 .

Distribution of alkoxide-alcohol data Ind. Eng. Chem. Fundom., Vol. 1 1 , No. 3, 1972

327

Table 1. Experimental Ebulliorcopic Data

Run no.

1


2

3

4

5

Initial alkoxide concn, (Bo)kt m

Initial alcohol concn, (AD)& m

0.0026 0.0057 0.013 0.020 0.038 0.062 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.30 0.0272 0.0272 0.0272 0.0272 0,0272 0.0272 0.0627 0.0627 0.0627 0.0627 0.0627 0.0627 0.0990 0.0990 0.0990 0.0990 0.0990 0,0990 0.0023 0.0243 0.0540 0.1010

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0068 0.0197 0.0510 0,1419 0.3245 0,4805 0.0036 0.0249 0.0573 0.1075 0.2745 0.4930 0.0107 0.0457 0.0861 0.1395 0,2894 0.4383 0.3085 0.3069 0.3042 0.3003

Molecules/ aggregate at equilibrium, (Ndk

2.29 2.89 3.33 3.60 3.97 4.17 4.32 4.48 4.50 4.64 4.72 4.79 4.84 4.89 3.35 2.76 1.94 1.44 1.27 1.23 4.09 3.73 3.13 2.40 1.68 1.45 4.28 3.90 3.30 2.72 2.01 1.75 1.09 1.28 1.59 2.11

Aggregate concn, (Ydk

((Bdk

=

f

(Ao)k)/ (Ndk

0.00114 0.00197 0.00390 0.00556 0.00957 0.01487 0.02083 0.02679 0.03333 0.03879 0.04991 0.05010 0.05579 0.06185 0.0102 0.0148 0,0403 0.1174 0.2769 0.4128 0.0162 0.0235 0.0383 0.0708 0.2006 0.3832 0.0256 0.0371 0.0561 0.0879 0.1935 0.3074 0.2851 0.2588 0.2213 0.1902

Run no.

5 6

7

8

9

10

Initial alkoxide concn, (Bdkt

m

0.1542 0.2063 0.2544 0.0032 0.0313 0.0619 0.0947 0.1489 0,2001 0.0049 0.0168 0.0292 0.0703 0.1084 0.1784 0,2252 0.0032 0.0080 0.0150 0.0524 0.1257 0.2321 0.0027 0.0082 0.0357 0.0764 0.1297 0.1726 0 0009 0,0029 0.0106 0.0269 0.0434 0.0602 0.0873

Initial alcohol concn, (A&, m

0,2959 0.2915 0.2875 0.1657 0.1634 0.1614 0.1593 0.1558 0.1525 0,0892 0.0889 0.0887 0.0877 0.0868 0.0853 0.0842 0.0011 0.0027 0.0050 0.0174 0.0418 0.0772 0.0027 0.0082 0.0358 0.0765 0.1294 0.1730 0.0028 0.0087 0.0317 0.0807 0.1301 0,1809 0.2618

Molecules/ aggregate at equilibrium, (Na)k

2.71 3.30 3.75 1.12 1.54 2.08 2.68 3.59 4.16 1.19 1.45 1.78 2.93 3.68 4.29 4.51 1.95 2.43 2.80 3.66 4.19 4.55 1.62 2.04 2.76 3.11 3.39 3.55 1.29 1.38 1.51 1.67 1.78 1.84 1.93

Aggregate concn, (Ydk = ((6O)k f ( A o ) k ) / (No)k

0,1661 0.1508 0,1445 0.1508 0.1264 0,1073 0.0948 0.0849 0.0607 0.0816 0.0729 0.0662 0.0539 0,0530 0.0615 0.0686 0,0022 0.0044 0.0071 0.0191 0.0400 0.0680 0.0033 0,0080 0.0259 0,0492 0.0765 0.0974 0,00287 0.00841 0,0280 0.0644 0.0975 0.1210 0.1809

tabulated in Table I. Runs 2-4 were carried out a t constant (Bo) and are shown in Figure 2; runs 5-7 were carried out a t constant (A,) and are shown in Figure 3; and runs 8-10 were carried out a t a constant ratio of (A,) to (Bo). The general shape of the three-dimensional surface, taking the state of aggregation as the dependent variable and the concentrations as the “independent” variables, may be visualized from Figures 2 and 3. 5

1 Bo (ALKOXIDE CONCENTRATION1 MOLAL

-I\

No (AVERAGE MOLECULES

T 11

\\ I

I

.01

Figure 3. Aggregate changes with alkoxide concentration at constant alcohol level

0.10 Ao(ALC0HOL CONCENTRATION I MOLAL

I .O

Figure 2. Aggregate changes with alcohol concentration at constant alkoxide level 328

Ind. Eng. Chem. Fundam., Vol. 11, No. 3, 1972

Prior to analyzing the data in terms of the chemical equilibrium model, the experimental error was determined by standard statist,ical means. There is a restriction on the randomness of the data since they were generated by gradually increasing the initial concentrations of hlcohol and/or alkoxide and observing yo a t equilibrium. The sources of error within a run are small relative to the between run varia-

Table II. Analysis of Pure Alkoxide System (Run 1) Model

Mo = 5

c (e,)k2/(24 - p - 1)

M = 6

M = 7

M* =

a

M = 9

24

RluM,N*=

k=l

SY2

28.4

4.42

3.78

3.91 (1)-(2) Indisting. 0.098 0.35

3.92 (2)

3.93 (2)

0.087 0.35

0.093 0.35

168 2655 5681 318 178 494 829

192 2352 5424 397 156 329 515 694

Fz4--p-1,8(0.05) = 5 Termination condition

3.89 (1)

3.90 (1)

DT

0.23 0.35

0.131 0.35

2478 339 275 76196

275 6185 1584 196 3268

Acceptable DT a t 95% significance level Kzo Parameter values a t K30 termination K40 Kso K60 K70 KW

KgO Downloaded by UNIV OF GEORGIA on September 5, 2015 | http://pubs.acs.org Publication Date: August 1, 1972 | doi: 10.1021/i160043a007

6.82

tions introduced by different operators, stock solutions, or changing environmental conditions. By considering the intersection points of the various runs shown in Figure 1, the between run variations are determined.. Duplication of runs by both the same and different operators permits estimation of reproducibility and within run errors. Considering these factors, the following error estimates were obtained from seven reproducibility experiments for the pure alkoxide system and 16 reproducibility experiments from the alcoholalkoxide [A-B] system sY2 = 2.92 X where sy2 is the estimate of the variance of y, r y 2 .The method described in this paper will be applied first to find the best model and parameter estimates for the pure alkoxide system (run 1) and then for the -4-B system. Pure Alkoxide System. As was mentioned above, this system was analyzed previously without the aid of the statistical methods described herein (Blau, et al., 1970). It will now be shown that the intuitive procedure for model selection used before is consistent with the statistical criteria. The system is a special case of the total system I in which A0 = A = 0. All equations in the general development are valid on making the substitution. Using 24 data points, the optimization was carried out for models such that .M = 5, 6, 7, 8, 9, starting from initial approximations K,o(O) = 1000, i = 2, . . , Ai. The results are presented in Table 11. It is seen that the optimization algorithms for = 5, 6, 7 are terminated with RM,o*> 5 indicating models that the models are too small, while those for models ,M = 8, 9 terminate with R,,o* < 5 . The best model by definition is thus Ji = 8. This is also the model in which the abnormal magnitude of the highest association constant has almost vanished (Blau, et al., 1970). The choice of the best model does not preclude the possibility that the physical system consists in reality of larger aggregates. It says only that larger aggregates are not required to explain the data within experimental error a t the 95% confidence level. The normality assumption was examined for the residual error distribution of each model a t the minimum by both the x2 and Kolmogorov-Smirnoff (DT statistic) goodness of fit tests. As seen by the value of DT (Lindgren and McElrath,

189 3290 4429 242 311 1244

3.72

1966), no departure from normality was detected a t the 95% confidence level for any of the models. Thus the use of the F statistic is valid. Alcohol-Alkoxide System. Inspection of the data in Table I and Figures 2 and 3 provides some useful information relative to the type of model which needs to be examined. The alcohol by itself is essentially monomeric (average measured value is 1.089 over the whole concentration range), which makes it unnecessary to consider equilibria of the type 2A $ Az. Addition of small amounts of alkoxide to solutions of high alcohol concentrations yields data which suggest that a maximum of two alcohol molecules associates with each alkoxide moiety except a t very low alkoxide concentration. Thus species B,A, in which j > 2i are probably unnecessary. Addition of small amounts of alcohol to high concentrations of alkoxide causes gradual reduction in average state of aggregation. This suggests that alcohol molecules associate with the alkoxide aggregates without causing massive destruction of the aggregates, and that species such as BtA,, where i > 1, will be required in the model. The number of equilibrium constants implied by the above reasoning is large and some simplification would be desirable. For instance, it would be helpful if certain concentration regions could be examined in which only a porbion of the equilibria are important. Values obtained for the corresponding equilibrium constants could then be used as known parameters in analyzing the rest of the data. This has been done for the pure alkoxide system, Le., the region where (Ao) = 0. Unfortunately, there is no other experimentally accessible region where it can be assumed that a well-defined subset of equilibria are the only important ones. Consequently, the whole A-B system must be examined a t once. The process of choosing the “best” model for the A-B system is not straightforward since both size and shape are involved. Analysis of the pure alkoxide system indicated that M* = 8, but that the concentration of Bs was small over the feasible concentration region. iiddition of alcohol tends to reduce the number of alkoxide moieties in the aggregates so that i1P = 7 is probably appropriate for the mixed system. There is no compelling evidence which can guide the choice of the maximum value of j for each i in BiAf and the computaInd. Eng. Chem. Fundam., Vol. 1 1 , No. 3, 1 9 7 2

329

Table 111. Analysis of Alcohol-Alkoxide System (Runs 2-10) Model (7,N)

N = l 18

RM9N*

=

N = 2

N = 3

465

133

N = 4

N = 5

N = 6

N = 7

(e,)nz/(165 - NM - M)

k-1

1200 SY2

F M - N M - M , ~ ~(0.05)

= 5 Termination condition Number of equilibrium constants

DT

Acceptable DTa t 9570 significance level

2.10 Not app. 13 0.31 0.17

2.10 Not app. 20 0.28 0.17

21.6

2.10 Not app. 27 0.22 0.17

2.10 1 34 0.12 0.17

Table IV. Optimal Equilibrium Constants for Various N (M

4.18

1.91

1.35

2.11 1 41 0.10 0.17

2.11 2 48 0.18 0.17

2.10 2 56 0.12 0.17

= 7)

Equilibrium constants

Klj*

K*j *

Kaj*

Kdj*

Kej*

Klj*

j = O

...

189

3290

4429

242

311

1244

N = l

j = l

18

14

57

1654

1154

115

1140

N = 3

j = 1 j = 2 j = 3

11 17 15

12 13 12

32 32 20

77 155 109

109 71 32

34 24 16

80 41 24

11 10 10 9 8 9

30 27 19 12 8 8

396 65 23 9 4 6

82 46 34 15 6 7

29 22 1 13 10 10

39 27 21 14 10 10

Model


Pure alkoxidea

N = 6

(1

j = 1 4 j = 2 5 j = 3 4 j = 4 2 j = 5 3 j = 6 6 Values for Kio* are the same for all the models.

tional burden of examining all possibilities is not warranted. A heuristic approach has therefore been taken wherein the model has been enlarged stepwise by increasing N in the rectangular model (7, N). The results of the optimization are recorded in Table 111. For the smaller models N _< 3, the model error is considerably larger than experimental error and both the x 2 and Kolmogorov-Smirnov test show a significant

-0

IO

-.

20 AX) 40 60 NUMBER OF EOUlLlBRlUM CONSTANTS

70

Figure 4. Validity of statistical criteria for alcohol-alkoxide system 330

Ind. Eng. Chem. Fundam., Vol. 11, No. 3, 1972

K5j*

departure from normality a t the minimum. For N 2 4, however, these tests fail to reject the normality assumption so that the ratio RM,.,,*follows the F distribution and tests for model adequacy are applicable. Figure 4 amplifies Table 111. As the number of equilibrium constants (model size) increases three distinct regions are encountered. When N 6 3, that is 27 equilibrium constants or less, the modelling error is so large that the statistical criteria are invalid. In the second region 4 _< N 5 5, where the number of equilibrium constants is between 28 and 47, a significant difference between the modelling and experimental error exists as specified by the now valid condition R M a , , * > 5 . Finally for N 2 6, the number of equilibrium constants greater than or equal to 48, it is impossible to distinguish the modelling error from the experimental error since R M , N * < 5. Thus i t is established that the model need be no larger than (7,6). The equilibrium constant estimates for three of the models are given in Table IV. Some further information may be obtained by considering the starting vectors for each model size. For the optimizations of models (7,1), (7,2), (7,3) and (7,4) all of the equilibrium constants except K i Owere started a t 10. For (7,5) and (7,6) the optimum values for the next smaller model were used while the new constants were placed a t 10. A general pattern appears in which K , , is adjusted the most when 1 is small and when i corresponds to the largest Kto. Some of the constants were not changed appreciably, which implies that the corresponding species are unimportant in the calculated equilibrium system. This indicates that the model should be reduced in size by eliminating those constants. Another test of this came by using nonrectangular models and a starting vector in which all values except Kio were placed a t unity.

Table V. Nonrectangular Models Starting with K,j

=1

1,

(i

....Mij > 0 )

Equilibrium conrtantr

i

Klj*

0

Pure alkoxidea

1 2

Model P

Kaj*

168

2655

kj* 568 1

...

3.9 1.5 L-----1.1 1.0 1.0

1

296.7 45.9 36.6 9.0 2.6

...

1.0

I

1.3

717 2670

9398 408000 1032

2.0 1.3 1.5 r--1-1---1 ------a 1.1 1.0 ... 1.0

3 4 5 6

... ... 359000 44.1

1

Model Q

Kpj*

2


Values for

kio*

I

L------l

0.0 ... are the same for both models.

3 5

{

Ksj*

Kej*

K7j*

Ksi*

318

178

494

829

182.9 31.3 14.4 3.0 :--1-1-

’

1.0

- -9.6 - - - - -5.8 - - - - - 12.8 - -I

-’I

1.6 1.1 1.0 1.o 1.0

1.4 1.1 1.0 1.o 1.0

1.5 1.1 1.0 1.o 1.0

... ... ...

...

...

...

...

... ...

... ...

...

...

... ...

This corresponds to assuming that there is no free energy change for the solvation reactions. The results are shown in Table V, Model P, where the unimportant constants are clearly apparent and are separated by the dotted boundary. Chemical reasoning still casts some doubt as to the physical reality of the model. The manner in which alkali metal alkoxides tend t o aggregate is shown in the following diagrams indicating the probable structures of alkoxide aggregates.

magnitude of the equilibrium constants will not change radically as long as the model is big enough, and further work does not seem warranted. The chemical implications and a more detailed discussion of the nature of the physical system described by the model will be deferred to a subsequent paper.

M-0-R

A methodology has been developed for estimating parameters and distinguishing rival models for general chemical equilibrium systems. The method has been shown to be applicable for the analysis of ebullioscopic data which gives information regarding the nature of alkoxides in media of low dielectric constant and the effect of added alcohol on the alkoxides. The tetramer of the alkoxide is the preferred, but not exclusive, state of aggregation of the alkoxides over a wide range of concentrations. The alkoxide aggregates may be solvated by a small number of alcohol molecules, regardless of the size of the aggregate. Limits have been placed on the number of equilibria which are necessary to describe the data. Equilibrium constants have been estimated which allow the model to fit the data within +6y0 over the entire concentration range of alcohol and alkoxide which was studied. Many models exhibit mathematical properties similar to those of the models examined here (Berthon and Luca, 1970; Tucker, et al., 1969). The techniques developed in this paper can be used both to formulate these problems as unconstrained nonlinear programming problems and to achieve statistically meaningful solutions.

R - 0 - M/

/

Bz R R M-0’

R-O< M-0 ‘R B3

B4

In Bz and B3 the oxygens each have one pair of electrons available for hydrogen bonding with an alcohol. In Bb they do not. Higher aggregates or a linear form of B4 might have a limited number of bonding sites. h minimal model from a chemical point of view probably should include the species B A , BAz, B ~ A IB , z A ~B&, , B A , and BgAa. Furthermore, a starting vector for the optimization could assume a free energy change of 4-6 kcal/mol for each hydrogen bond formation. The results of optimization with this model are shown in Table V, Model Q. The characteristic results of too small a model are apparent in the way certain constants are made very large. Furthermore, the general shape of the calculated surface does not conform to the data, particularly in the low concentration region, and the residuals are not normally distributed. Based on this analysis, therefore, the higher aggregates are indeed solvated by alcohol molecules. There is no obvious reason why, if the B4 aggregate can be solvated, the higher aggregates cannot also be solvated. Therefore, i t is concluded that Model P in Table V after excluding the equilibria below the dotted line is the “best” of the models examined. That is, it is the smallest model which adequately explains the data. Further examination of models can be envisioned. It is apparent, however, that the

Conclusion

literature Cited

Bard, Y . , Lapidus, L., Catal. Rev. 2 ( l ) ,67 (1968). Berthon, G., Luca, C., Chim. Anal. (Paris)52 (4), 392 (1970). Blau, G. E.. KlimDel, R. R., Steiner, E. C., IND.ENQ.CHEM., F U N D A M . 9, 334-(1970). . Blau, G. E., Klimpel, R. R., Steiner, E. C., “Parameter Estimation and Model Distinguishability of Physicochemical Models at_hemical Equilibrium,’’ to appear in Can. J . Chem. Eng., 1Y7Y.

Box, G. E. P., Hill, W. J., Technometrics 9, 57 (1967). Box, W. J., A p p l . Statist. 18, No. 3 (1969). Brinkley, S. R . , J . Chem. Phys. 15, 107 (1947). Davidon. W. C.. U. S. Atomic Energv Commission Research and

Develo ment’Report ANL-5990 o“959). Draper, R., Hunter, W. G., F m e t r i k a 53, 525 (1966). Dra er, N . R., Smth, H., Applied Regression Analysis,” d l e y , New York, N. Y., 1966.

8.

Ind. Eng. Chem. Fundom., Vol. 1 1 , No. 3, 1972

331

Fletcher, R., Powell, M. J. D., Comput. J . 6 (2)’ 163 (1963). Hill, W. J., Hunter, W. G., University of Wisconsin, Department of Statistics Technical Report 69 (1966). StaHogg, R. v.7 Craig, A. T.,L‘lntroduction to tistics,” 2nd ed, Macmillan, New York, N. Y., 1965. Kittrell, J. R., Hunter, W. G., Watson, C. C., A.I.Ch.E. J . 12, 5 (1966). Lindgren, B., McElrath, G., ‘‘Introduction to Probability and Statistics,” 1st ed, Macmillan, New York, N. Y., 1967. Reilly, P. M., Can. J. Chem. Eng. 48 (2), 168 (1970).

Shapiro, N. Z., Shapley, L. S., J . SOC.Indust. Appl. Math. 13 (2), 353 (1965). Tucker, E. C., Farnham, S. B., Christian, S. D., J . Phys. Chem. 73, 3820 (1969). White, W. B., Johnson, S. &, Dantzig, I., G , B., J , Chem. Phys. 2 8 , 751 (1958). Wilde, D. J., Beightler, C. S., “Foundations of Optimization,” Prentice-Hall, Inc., Englewood Cliffs, N. J., 1967. Zoutendijk, G., J . SOC.Indust. Appl. Math. 4 (I), 194 (1966). RECEIVED for review April 21, 1971 ACCEPTED April 27, 1972

A Model for Predicting the Permeation of


Hydrogen-Deuterium-Inert through Palladium Tubes

Gas Mixtures

Frank J. Ackerman” and George J. Koskinas Lawrence Livermore Laboratory, University of California, Livermore, Calif. 94550

This paper describes a mathematical model of the permeation of hydrogen-deuterium-inert gas mixtures flowing through long palladium-silver alloy tubes. Permeation experiments using Dz-He, Hz-Dz, and Hz-DzHe mixtures agree well with the model. The permeation experiments were performed with a Pd-25 wt % Ag tube at 400°C under a total driving pressure of 1000 psia and a total back pressure of 20 psia. The model can readily be extended to other conditions.

Palladium has been used to produce high-purity hydrogen and its isotopes in the laboratory since the early 1920’s, but the hydrogen-palladium permeation process has been used commercially only in the last 10 years. A major cause of this delay has been the structural weakening of the palladium by a phase change in the hydride. Now, however, the development of palladium-silver alloys has overcome this problem. Even with such a long history of laboratory application, few investigators have reported hydrogen or deuterium permeation studies at common industrial operating pressures (above 50 psig) (deRosset, 1960; Hunter, 1965; Rubin 1966). However, many commercial uses for this process have been suggested. They fall into three categories (Darling, 1958; Serfass and Silman, 1965) : (1) purification of the hydrogen supplied to a process; (2) recovery of hydrogen from process streams for recycling or as a by-product; and (3) isotope separation. I n order to be practical, processes for purifying or recovering hydrogen must be designed to recover most of the hydrogen. No studies involving the substantial recovery of hydrogen from gas mixtures have been found. deRosset (1960) has studied the permeation of H2-Hz and H2-CH4 mixtures] but under conditions where most of the hydrogen remains in the mixture. The use of palladium for separating hydrogen isotopes has been studied by several investigators (Lewis, 1967). In these studies, the separation is accomplished by gas chromatographic techniques that use the difference in the solubilities of the isotopes to accomplish the separation. This is essentially 332

Ind. Eng. Chem. Fundam., Vol. 1 1 , No. 3, 1972

a batch separation method. The use of palladium in a continuous separation process has not been reported. This study is concerned with the efficient recovery of Hz or Dz from gas mixtures containing one or more nondiffusing contaminants. The geometry of the system is similar to a single-tube heat exchanger. The entering mixture becomes depleted in hydrogen as it flows through the center of a heated palladium-silver tube, and pure hydrogen is withdrawn from the outside of the tube. A mathematical model of this process has been formulated to examine the influence of the process variables (composition, temperature] driving pressure, back pressure, bleed rate, and tube geometry) on the recovery of hydrogen or deuterium. The isotope separation that is obtained during the simultaneous permeation of two hydrogen isotopes can also be calculated with this model. The model has been tested with mixtures of DZ-He, Hz-Dz, and HrDz-He. Helium is a convenient nondiffusing contaminant since it is readily available in pure form and has properties similar to hydrogen and deuterium, thereby simplifying the PVT calculations. The permeation experiments were performed with a Pd-25 wt % Ag permeation tube a t nominal conditions of 4OO0C, 1000 psia of total driving pressure, and 20 psia of total back pressure. Theory

Pure-Gas Permeation. Our purpose was to develop a mathematical model for predicting the permeation of gas mix-

Equilibrium Constant Estimation and Model Distinguishability

Recommend Documents