Stochastic modelling of reactors for process analysis and design — I isothermal steady state well mixed and plug flow reactors

ClumicolEngineering Science, 1973, Vol. 28, pp. 1257-1272.

Pergamon Press.

Printed in Great Britain

Stochastic modelling of reactors for process analysis and design-I Isothermal steady state well mixed and plug flow reactors J. E. BERRYMAN

and D. M. HIMMELBLAU

Department of Chemical Engineering, The University ofTexas at Austin, Austin, Texas 78712, U.S.A. (Received8

June 1972; accepted6

October

1972)

Abstract- In most conventional methods of reactor analysis and design the stochastic nature of the process inputs and parameters is ignored. Safety factors are added to the deterministic results to allow for uncertainty. Stochastic treatment of a reactor model is more ditlicult than a deterministic treatment but can provide more significant information concerning the design parameters. Both linear and nonlinear models of isothermal well-mixed and plug flow reactors have been investigated in which random variables were introduced to represent random input concentrations and flow rates, as well as uncertainty in the model coefficients. This study indicates under what circumstances the sample estimate of the expected value of the output of a stochastic model of the well-mixed tank differs from the output of the related deterministic model based on replacing the random variables by their expected values. It also shows how to compute safety factors for design and in what ways the response of the stochastic well-mixed tank differs from the response of the stochastic plug flow model. INTRODUCTION MATHEMATICAL

models developed to describe simulated or actual reactors formally consist of various types of algebraic and differential equations together with the appropriate boundary and initial conditions. Although the equations employed in design evolve from the fundamental continuum equations describing the conservation of mass and energy, so many simplifying assumptions must be made in the evolution that the resulting process design equations must incorporate one or more empirical coefficients if they are to have any predictive ability. The values of the coefficient are obtained from actual operation or from laboratory and pilot plant data. Such models provide a convenient and economical basis for studying the effect of changes in variables and operating conditions upon the output of a real process in so far as the model can represent the real process. In formulating a mathematical model of a process, two different approaches can be discerned; one deterministic and the other stochastic. Stochastic models evolve through two main routes. Some stochastic models can be derived by directly employing probabilistic concepts to small physical subdivisions of the process. An alter-

nate method of obtaining stochastic models is to take a deterministic model and transform it into a stochastic model by introducing random variables either as inputs or as coefficients, or both. It is this alternate approach that will be used here. Neither route of approaches is more “correct” than the other, because both types of models are approximations in any case. What is important in each approach is to make the introduction of randomness in the model reflect what is happening physically in the process in so far as possible. The engineer using a model would like to be able to relate the statistics of the response of the model to the statistics (such as the mean and standard deviation) of the random parameters and random inputs of the model. It would be nice if he could solve the stochastic model directly to obtain either the probability density functions of the dependent variables or the statistics of the variables. Bharucha-Reid [5] presents a brief survey of the theory of stochastic equations for algebraic equations, difference equations, differential equations, and integral equations having random variables for inputs and for coefficients. Unfortunately, but not unexpectedly, the formulation of the more representative stochastic

1257

J. E. BERRYMAN

and D. M. HIMMELBLAU

models results in equations with random parameters whose dependent variables are nonMarkovian, non-stationary or both. Consequently, while it is usually not too difficult to write down the basic (differential) equations describing a reactor, or other process, how to obtain a complete, partial or approximate, analytical solution for a stochastic model process seems to be a stumbling block. For these reasons, only a few attempts have been made to apply an analytical stochastic analysis to reactors and other models of interest to the chemical engineer (Acrivos [ 11, Aris and Amundson [2], Homan and Tiemey [6], Katz [7], King [81). It is possible, though, to obtain some insight into the stochastic character of models via estimates of the expected values and standard deviations of the random responses through Monte Carlo simulation. This investigation used Monte Carlo simulation to obtain sample estimates of the expected values and variances of the output variables of reactor models incorporating random inputs, random coefficients, and random boundary conditions. The relation between the mean and standard deviation of the output variables was computed and interpreted in terms of reactor performance and computation of overdesign factors. Both linear and nonlinear models of continuous flow well-mixed tank reactors and plug flow reactors are used to illustrate the physical and computational features encountered in stochastic analysis that differ from the usual deterministic treatment. In some instances the results obtained differ from those reported in the literature, and the reasons for these differences are explained. Studies using stochastic models can be helpful in deciding if the model coefficients have been determined with sufficient accuracy for the use in design or simulation. They can also materially assist in providing a more quantitative value for the magnitude of the safety or overdesign factor to use in the design of process equipment. REACTOR

MODELS

STUDIED

Table 1 lists the types of reactors investigated.

Table 1. Models studied Single, first order differential equation: A. Continuous flow, well-mixed, isothermal tank reactor 1. First order, irreversible reaction- Model 1A 2. Second order, reversible reaction- Model 1B B. Steady state, plug flow, isothermal reactor 1. First order, irreversible reaction- Model 2A 2. Second order, reversible reaction- Model 2B

(a) Model 1

Well-Mixed,

Isothermal

Reactor

CAf qf “f +-z (b)

Model 2

Steady-State, Reactor

Plug Flow, Isothermal

Fig. 1. Two reactors.

The continuous flow, well-mixed, isothermal tank reactor represents a good starting point for stochastic analysis since the deterministic material balance for one component A is a simple ordinary differential equation in which the independent variable is time: ydC/, -= dt

4iC,4*-4fC.4 -4&A + vR(kv C*) (1)

with initial condition c, (t = 0) = c,,. The volume of the reactor remains constant; however, the moles of product may be different from the moles of reactant and for a gas phase reaction this results in the flow rate at the outlet being different from the rate at the inlet. To

1258

Reactors for process analysis and design- I

account for any volumetric flow rate change it was assumed that the change in rate is proportional to conversion. The outlet flow rate is related to the inlet rate by G= qr(l +&CA)

(2)

The unit time will be the residence time as determined by the reactor volume and mean inlet flow rate. By following a similar analysis for the isothermal plug flow reactor, the differential equation that represents a material balance on reactant A is

which becomes after substitution

NAi- NA

GA=

=

A1

After substituting Eq. (1) one obtains

C.4

(3)

(8)

with boundary condition CA(Z = 0) = CAi.

(4) Introduction of Eqs. (2)-(4) into Eq. (8) plus the following expression that permits a change in flow rate

qi* = qJV and Eq. (2) into

q” CAi(CAi-CA) t c,,+E,c,

dt=

-

CAi+ EACA

qiCAi( 1 +EA) CAi+EACA .

“’

dCA

CAi

N,

d+=NkvC)

gives

(9

+R(k,CA)

G(l ‘*

CA (t = 0) = CA,,,

+EA)

dCA

(CAi+EACA)2z=

R(k’

‘)

(9)

with Introducing the following rate expressions

CA(Z = 0) = CAi.

Model A: R (k, CA) = -kCA Model B: R (k, CA) = -klCACB i- k&p for a first order irreversible reaction and a second order reversible reaction, respectively, leads to

For a constant crossectional area reactor let x = z/L. Then the models for the isothermal plug flow reactor are Model 2A

dCA -_

Model 1A

dx -z

dCA

cAi(cAi-CA)

_kC

SF

CAi+EACA

cAi(t = 0) = CA,

dC.4 -=

(6)

qz”C;i(l

Model 1B

dCA -=4: dt

CAi(CAI-C~)

c

+E

Ai

- k,CA (CBi CA(t = 0) = CAQ

c

A

+kz(CPi+CA*-CA) A

CAi

+

CA)

(7)

+EA)

(10)

Model 2B

(CA,+EACA)2[k~(C,,+CAi-cA)-kk,CA(C,i-cAi+C~)I

dx

Cii(l

CA(X = 0) = CAie

A

dt

-=

-k c,(c,i+&cA)2

(11)

+EA)

In solving the deterministic equations for a plug flow reactor it is easier to work with fractional conversion than concentration. However, conversion was not used as the dependent variable in this study because of the need to maintain the identity and statistics of the random variables which can actually be measured. If the statistics

1259

J. E. BERRYMAN

and D. M. HIMMELBLAU

for the variable “fraction conversion” had been obtained, they could not have been easily reduced to the statistics of the output concentration except for the case in which CAiwas a constant. MONTE

CARLO

SIMULATION

The Monte Carlo method was used to obtain the statistics of the model response in terms of the statistics of the random inputs and coefficients. The Monte Carlo method essentially consists of the generation of a large number of repetitive solutions of the process model from which sample statistics of the response can be calculated. Random inputs and random coefficients were introduced into the model equations in a manner which conceptually simulated the way that uncertainty in the respective input or coefficient enters the real process. The question thus arises: how should the aspect of randomness be introduced into a deterministic model? When the objective of the analysis is to calculate the effect of stochastic inputs and coefficients on the response of the process model, then it is necessary to be quite careful to relate the model to physical reality lest the results of the analysis be erroneous. Processes that are represented by the same deterministic differential equation and boundary or initial conditions can in fact be quite different in a physical sense, and hence must be treated differently if they are to be represented by stochastic models. A fluid element entering the well-mixed vessel becomes mixed with the fluid already in the vessel, and the input variable CAiassociated with the fluid element entering the vessel is not isolated from other values of CAi in the time sequence as would be true in the plug flow reactor. Each value of C,i in the time sequence of discrete values, as it mixes, causes the random error associated with the variable to mix with the accumulated mix of random error in the vessel from the start up. Random values of qi do not mix since the volume of the vessel was assumed to be constant and the instantaneous flow rate out was equal to the flow rate entering. Thus, in the well-mixed tank model qi and CAi were the random variables that were functions of

time. On any one run of a series of repeated runs, the differential equation(s) in the model were solved numerically. Values of the time varying random variables were represented by a sequence of discrete (piecewise constant) random numbers. The net effect of this treatment of the well-mixed tank model was to discretize the time varying random variables and solve the equations in the deterministic model for each time interval At. During each time interval At (selected to be of the order of 0.005-O-04 residence times in order to achieve the desired numerical accuracy), fourth order Runge-Kutta integration was used. During one time interval, At, the random parameters were regarded as having fixed values that might or might not change in the next time increment. As long as At was very much less than the time constant of the tank (in the same units), this treatment corresponds to what presumably is the physical mixing existing in the tank. In the plug flow model qi and CAiwere random variables that were functions of time. In the more complicated cases, u was also permitted to vary imitating a gas phase reactor with a change in number of total moles. To correspond to the physical nature of the mixing in the plug flow vessel, each input for the period At was applied to the fluid element entering during the interval At. However, once the fluid element started to travel down the vessel, it became isolated from all other input variables. Consequently, a deterministic solution could be used for linear differential equations, or otherwise a numerical integration was used, to obtain a deterministic solution for the output concentrations. By this treatment, in the plug flow model, it made no difference whether the input variables were time dependent or not. Each of the values of the random variables employed were computed by adding the generated random variable to a deterministic quantity, as for example CAi = ( CAi) + random error.

(12)

Table 2 lists the random variables used and their characteristic distributions. The method of

12!60

Reactors for process analysis and design-

I

Table 2. Classification of random variables

Function of time Models 1A and 1B: k k, 4T c*i CSi CP, c A0

Independent normal distribution

no no yes yes yes yes no

Models 2A and 2B: k k, 41’ C.4, Cl3i CPi

Independent uniform distribution

yes yes yes yes yes no yes

Yes

yes

Yes

Normal amplitude but auto-correlated no no yes yes yes no no

yes yes yes 5-s

no yes

Yes

yes

yes yes yes

Yes

yes yes

no

no

no no no no no no

*Always zero in this study.

generation of these variables and the validation of their characteristics is summarized in Berryman and Himmelblau [4] and described in detail in Berryman [31. The auto-correlated variable was not used in the plug flow model because its effect on any trial was exactly equivalent to that of the normal variable as should be clear from the discussion above on mixing. The coefficients in the models were assumed to be normal random variables only because they were presumed to be estimated from experimental data. It was necessary for the statistical analysis to obtain samples of the process output of adequate size in order to obtain a good estimate of the ensemble mean and ensemble standard deviation of the model response. The required sample size is related to the variance of the sample and, thus, to the variance of the ensemble. It was demonstrated by the use of the Chi-squared statistic that a sample size of 300 was more than adequate for the desired precision.

corresponding deterministic response, if at all? Because such large samples would be required to investigate this point by simulation, an analytical treatment of some of the simpler models is the most clear cut away to obtain an answer. If we consider the simplest plug flow model, Eq. (10) with EA = 0 and a first order reaction

dC,b) -= dx

-K,(x)

C,(x=O)

= CAi

4i* (13)

where q: , CAi and k are all random but not functions of x, we can demonstrate the difficulty of obtaining an analytical expression for the ensemble mean of CA. Multiply both sides of Eq. (13) by l/ ( CA1), and, after rearranging, take the expectation of each side

E

HOW THE EXPECTED VALUE OF THE RESPONSE IS RELATED TO THE DETERMINISTIC VALUE

(14)

The first question of importance to examine is: under what circumstances does the expected value of the output of a reactor deviate from the 1261

J. E. BERRYMAN

Since k and q$ are independent,

E[d kZ)l

and D. M. HIMMELBLAU

Eq. ( 14) becomes

= -E[k]E[l/q,*]

Note that R,ln(E[&])=

dx

R,lnw

(15) E[cA;;A;jo)]

=

= 1.

Equation (15) has the solution (after integrating and exchanging expected value and integral operators)

1.06

3

----

(CAi)

R’lnm=O

Y, =0,05

1.05-

c i ‘= N

(16)

1.04-

:

In order to obtain the desired E[ CA (x) /( CAi)] in terms of the expected values of CAI, k and @, it is necessary to express E[ln(CA(x)/(CA1))] in

W

te~sOf~n(E[C~(~)/(~Ai)l),E[ln(~Ai/(~Ai))l

g

in terms of ln(EICA1/(CAI)]), and E[l/qt] in terms of l/E[ qt]. These three relationships enable one to determine the difference between the deterministic solution for CA(x)/(CAi) and the ensemble mean solution for CA(x) /( CAO. Figure 2 shows the ratio of E[ln (2) ] to In (E[Zj) where 2 is any normal random variable as a function of the magnitude of 2 between 0 and 1 for coefficients of variation for Z of O-05 and 0.10. Note that if Z is less than 0.6 with a coefficient of variation of less than O-10, any error introduced by the exchange of In (E[ CA(x)/(CAi)]) for E[ln (CA(x)/(CAi))] k kSS than one perCent. Figure 3 shows how the ratio of E[ 1/Z] / l/E[ Z] varies with the coefficient of variation of Z. Note that the difference between using E[ l/Z] is less than one percent up to a coefficient of variation of Z equal to O-10. Replacing E[ln (CA(x)/ ( CAi ))I by RI In (E[CA(X)/ (cAi)l),

E[ln(CAi/(CAi))I

by

ij

1.03-

2 1.02-

I.01 -

l.OO!LZS 0

0.2

0.4

7

L

OS

0.6

Fig. 2. Ratio of E[Ln(Z)J to Ln(E[Zj) for a normal

ID random

variable. l-06-

&ln(E[CAi/(CAi)l)

and E[ I/q?] by Z?,/E[q?] in Eq. (16) gives

0 E[kl =XR”E[q,*l.

alo

Coefficient

(17)

of Variation

Q20 of Z,

yz

Fig. 3. Ratio of E[ l/Zl to l/E[Z] for a normal random variable.

1262

Reactors for process analysis and design- I

Both (qf ) and ( CAi) are constants because qi* and CAi were ergodic random variables. Also, k is a constant because it is not a function of time. Integration of Eq. (6) after the substitution of the relations for qf ( f) and CAi(t) gives

and Eq. ( 17) reduces to

The corresponding deterministic solution of Eq. (10) with the random variables replaced by their respective deterministic quantities is (19) The conclusion can be drawn that although the expected value of the output of the simple plug flow model is not the same as that obtained from the corresponding deterministic solution, the ratio RJR, will be almost unity and can be neglected in most practical cases. When the conversion in the plug flow reactor is small, the ratio will not be one but the difference from one should still not be too large. If one attempts to solve Eqs. (10) and (11) with EA not equal to zero by taking the expectation of each side as was done with Eq. (13), one finds that the equation cannot be rearranged to separate CA from the random variables CA*,k and 4:. As a result, it is impossible to show analytically whether the expected value of CA(x) approaches the value obtained using a corresponding deterministic solution in which the expected values of CAi, k and 1/qF, respectively, are introduced as deterministic values. Thus, Monte Carlo simulation is required. If we consider next the simplest isothermal well-mixed model, Eq. (6), with EA = 0, for a first order reaction, because of the character of the mixing that takes place, we must integrate as follows. If k, q:, and CAiare all random variables, split the random variables into two parts, the ensemble mean and a random deviation Zi (t)

+_(* Z&f)&,,(f)

z,,(f)C~(t)

dt. cw

If y designates the deterministic concentration of A, the deterministic integral equation obtained from Eq. (6) is I U~t*)dy(t)

+[k+(ql*)l

/d’~(f,

dt

JOtzdf=O.

-(qf)(cAi)

(21)

Subtraction of Eq. (21) from Eq. (20) gives

+

=

c,(t) df-j-Ot* y(t) (,*,)(,,ot2

W+

(cAi>

+

d(t)

df-I,”

sd’

Jot’

z,:

z,,(t)

(t)ZcAi(d

dt+

(q?)

dr-

I,”

dt

Jot2 ZCA~(~)

-&f(t)

CA(t)

> dt

dt.

= (sT)+Z,,(O (22)

cAi(t) =

tcAi>

E[Z,tl = 0 E&,1

= 0.

+zCA,(f) Let ZcA (t) = CA(t) -y(t), where ZcA(t) is a measure of the divergence of CA(r) from the deterministic solution, so that Eq. (22) becomes 1263

J. E. BERRYMAN

0 +I,”

Z,? (t)Zc,,(t)

and D. M. HIMMELBLAU

response of the corresponding deterministic solution. However, because of the small -yk values used in this study, the divergence of C, (t) from y(t) as a result of a random k could not be detected by significance tests. To sum up, for either reactor, although in principle the ensemble mean response might not be the same as the corresponding deterministic solution, in so far as design is concerned, there was no practical difference in the two responses.

@ dl-[;*

Z&A(t)

dt.

RESULTS ISOTHERMAL

(23)

Evaluation of terms 3, 4 and 5 on the righthand side of Eq. (23) is possible under some conditions. If the time period of integration is long enough, the time averages of Z,: , ZcAi and Z,:ZcAi (since Z,: and Zc,, are independent) approach zero by definition. (The use of ergodic random inputs and parameters in this study permits the interchange of ensemble averages and time averages.) Note that the remaining terms, 1, 2 and 6 do not contain (C,i) or ZcAi, hence, the divergence of CA(t) from the deterministic solution is independent of the noise in CM. In order to obtain an estimate of the ensemble average values of the integrals in Eq. (23), 400 simulations were run using normal random variables for CAi and 4; in Eq. (6). During each run the values of each of the three integrals 1,2 and 6 of Eq. (23) were recorded at eight time periods (fz = 1,2, 3,4,5, lo,15 and 20 residence times), and the sample average of the three integrals determined for each r,. Term 1 was essentially the same for each t2 at - 0.0025, Term 2 approached zero as t2 approached 20, and Term 3 was essentially the average of Term 1 multiplied by tP. Thus, for the case in which @ and CAi were the random variables it was determined that the ensemble mean response of model 1A should essentially be the same as the deterministic solution. Extension of the above analysis to include instances in which k was a random variable indicated that E[ C,J( CA,)] should diverge from the

OF SIMULATION WELL MIXED

FOR THE TANK MODEL

The influence of the type of distribution of the random inputs and the coefficients on the distribution, sample mean response, sample standard deviation, and the range of CA/( CAi) was investigated for a number of cases of models 1A and B. Although all the simulation runs were carried out by calculating from an initial condition CA,,= 1.0 through the start-up period to the pseudo-steady state (5-20 residence times), of most interest were the distributions and statistics of the output variable in the steady state. Results for the unsteady state are reported in Berryman and Himmelblau [4]. 1. E#ect of the type of distribution of the random parameters on the distribution of the response, C.4lfCA)

The parameters k, CAI, and qT were individually and jointly permitted to the random. Chisquare tests were used to ascertain the results listed in Table 3. 2. Effect of the type of distribution of the random parameters on the sample mean and standard deviation of the response, CA/( CAi) Significance tests for comparing sample means were used to demonstrate that there was no significant difference in the means computed from the trials using the uniform generator vs the normal generator. F-tests for comparing variances demonstrated that the standard deviations produced in the response by the random parameters generated by either generator were not significantly different.

1264

Reactors for process analysis and design-

I

Table 3. Effect of the distribution of the random parameters on the steady state distribution of C.,/(C,,) Random variables Variable(s) Distribution k k 4i Yi 4i

Cdi

c*i

CAi k qr. C/ii k qi, CAi k,, k,, Gi, CB~.qi k,, k,, CM, C~ir qi

Normal

Response Uniform

Other

X

normal uniform normal uniform correlated normal uniform correlated normal uniform normal uniform

xt

X X X X X X X X

X

XS

X X

tFor (k) V/( qi) greater than 6; *For (k) V/( qi) greater than 1.

1A with EA = 0. If the sequence of discrete values of CAigenerated to represent the input concentration changed rapidly, the net effect on the output concentration was small, so that the output concentration approached a constant value. In effect, as a result of the mixing of the input fluid with the much larger volume of fluid in the tank, the well-mixed tank acted as a damper in the sense that the dispersion of the concentration of the outlet stream is less than the dispersion of the concentration of the inlet stream. The values of both s,/(C,J and yCA(CAL)are less than the corresponding values of sc,,/( CA0 and yc,,( CAi) of the input at any frequency of change. It was found that the values of Ykr ‘y,:, and -ycki had a linear or almost linear relation to the value of ycA so that plots such as Fig. 4 could be used to predict the dispersion when more than one parameter was simultaneously varied through use of the “propagation of error” formula. This finding makes it unnecessary to run large numbers of simulations-all the engineer need do to compute the confidence limits for CA in general is to use a few basic plots such as Fig. 4 and prepare a figure such as Fig. 5. Figure 5 shows (1) the pseudo steady state values of C,/(CAi) plotted as a function of the rate factor and (2) the upper confidence limit for 3. Results for Model

Frequency -.O-’


of Change 2

Changes/Residence

Time

-0-

IO

Changes/Residence

Time

-e-

50

Changes/Residence

Time

Changes/Residence

Time

-‘o-‘200

A2 (x08

/

0

0.10

Coeffment

of Voriotlon

0.20

of CAi

, 7,-.,

Fig. 4. Relation of ycAi to yc,/(C,J at (k)V/(qJ = 1 for Model 1A, E, = 0, pseudo steady state.

a confidence coefficient of 0.99 for the case in which Yk = 3/47= -rCAi= O-10. Upper confidence limits are shown for frequency of change of qp and CAi of 10 and 50 changes per residence time. The upper confidence limit for the case in which k was the only random variable (Yk = 0- 10) coincides with the line shown for a frequency of change of 50. Plots such as Fig. 5 can be used in the design of reactors by knowing or assuming the mean and standard deviation of the random

1265 CES Vol. 28 No. 6-C

h ._

0” ” ‘4

J. E. BERRYMAN

?_ 1.00 .

3

z5 E

-

and D. M. HIMMELBLAU

frequency of change of 50 we find that the rate factor is 3.98, equivalent to a volume of 39.8 ft3, or an overdesign factor of 33 per cent for a confidence coefficient of 099.

Upper Confidence Limit for Confidence Coefficient = 0.99 -

Frequency

of

Change

= IO

-

Frequency

of

Change

= 50

0

-

Sample Values of C*/

._ 0

--.005

0

Deterministic

2

Solution

I 4

Rate Factor,

I 6

8

1 IO

V/

Fig. 5. Effect of random k, qT and CAi on the pseudo steady state response of Model lA, EA = 0.

variables k, qi and C’Ai,and using the standard deviation of C,/( CAi) to set up confidence limits about (CA) / ( CAi) for a given confidence coefficient. The difference between the volume of the vessel at the upper confidence limit for (C,)/ (C,i) and the volume of the vessel at C,/(CAi) gives a reasonable estimate of the required overdesign factor at the given confidence level. As a specific example of the type calculation for overdesign, consider the following data: yh_ =

0.10

(k) = lO/min

-yqi= 0.10

( qi ) = 100 cubic ft/min

ycAi = 0.10

( CAi) = 3 moles Of A/ft”.

EA = 0 The frequency of change of qi and CAi is 50. Suppose that the desired conversion is 0.75. The deterministic solution at a concentration of (C,) / ( CA1) = 0.25 corresponds to a rate factor (k) V/ (qi) of 3.0; therefore, the deterministic volume of the reactor required is V = (3) (lOO/lO) ft3 = 30 ft3. Extending the line at C,/( CAI) = 0.25 to intersect with the upper confidence limit at a

4. Results for Model 1A with EA = 1 If the density of the fluid changes as a result of the reaction, the model equation becomes nonlinear. Two cases of Model 1A were run with EA = 1 to determine if the nonlinearity of the model yielded results different in principle from the results obtained for Model 1A with E.4 = 0. The cases were run with k, qT and CAi being simultaneously random and qi* and CAi having a frequency of change of 10 per residence time for each case. The nonlinearity of the model had no unusual effect on the dispersion of C,/(C,i) during the unsteady state and the values of qc, changed smoothly from 0 at t * = 0 to the pseudo state value. The values of qc,/( CAi) were slightly less than the values of qc,/( CAI ) for the corresponding cases in which EA = 0 because the effective frequency of change of qT and CAi was larger than 10. (The frequency of change was based on residence time as calculated from the entering flow rate and when E.., was greater than 0, the actual residence time was increased.) It was concluded that the nonlinearity of the model did not cause the character of the sample statistics of the response to differ in principle from the character of the sample statistics of the response of the linear model. 5. Results for Model 1B Studies of Model 1B can be separated into two categories, first, the cases in which k,, k, , q?, CAi and CBj were simultaneously random variables and second, the cases in which k, and k, were simultaneously random variables. The unsteady state sample mean response for each case investigated followed the deterministic solution and the calculated values of Tc,/( CAi) changed smoothly from zero to the pseudo-steady state values. No unusual results were observed. Overdesign factors based on the upper confidence limit can be computed as described above but it should be

1266

Reactorsfor

process analysis and design-

az 12,cW--- --g H - - 3 0 5 10”&ooo:

noted that they would be large, primarily, because of the nature of the reversible reaction; that is, a relatively small concentration decrease occurs as a result of a unit increase in the rate factor. It was concluded that the results for Model 1B were no different in principle from the results obtained for Model 1A. Specific results for Model 1B and appropriate figures for design can be found in Berryman [3]. RESULTS OF SIMULATION FOR ISOTHERMAL PLUG FLOW REACTOR MODEL

I

THE

Distribution

r ii

?, b B

/ 9:

r

E

oii

C‘4l(C‘4,)

3. Results for Model 2A with EA = 0 Figure 7 summarizes a typical set of results for

Normal

d L?j\ 7 87 ‘i?,\

# - r7 $ ._a

1

2. Effect of the type of distribution of the random parameters on the sample mean and standard deviation of the response, CA/( CAi) A considerable amount of data was collected for models 2A and 2B with EA = 0 and EA = 1, and the significance tests employed indicated that for the distributions used the sample means and sample standard deviations of the model outputs were not influenced by the type of distribution used for the parameters k, CAi and qr.

Variables Variables

r4,cQo-

b

In an effort to get some idea of the nature of the probability distribution function for CA(x), 10 sets of 10,000 values (a total of 100,000 values) of CA were computed at the reactor outlet for two cases: (a) CAi, 9: and k were uniform random variables and (b) CAi, qT and k were normal random variables. Each of the ten sets of values for both cases were tested by the Chi-Squared Goodness of Fit Test using the hypothesis that the 10,000 values represented a normal distribution, and none of the sets passed the test. Figure 6 illustrates the resulting frequency distribution for Model 2A. Details of the results of the simulations can be found in Berryman [3] but the type of distribution of the input random variable was found to have a significant effect on the distribution of the model response.

Random

rj 11

ri

1. Effect of the type of distribution of the random parameters on the distribution of the responses,

Random

Normal

Reference

i 21 Z ’

Uniform

I

-i 1

rl --1

I, I I I I I I I 9 I 9 8 11 1 I 8 0 18

V/

5

= 2

\

15

IO

20

Subinterval

Fig. 6. Distribution

of

100,000 values of (G-cc,) Model 2A.

-- -

Deterministic

for

cA/cAi

Random Variables

-0.005 0

I

2

Rate Factor,

3

V/
4

I 5

>

Fig. 7. Effect of random k on the response of Model 2A, E,, = 0.

Model 2A with EA = 0. The upper confidence level was calculated from the sample mean and standard deviation assuming that CA/(C,i) was

1267

J. E. BERRYMAN

represented

by a normal distribution because:

and

D. M. HIMMELBLAU

4. Results for Model 2A with EA = 1. The results for this model can be summarized as follows. The sample mean response of Model 2A with E.., = 1 for the cases in which random coefficients k and q: were used was statistically different from the deterministic solution; however, the difference was small. When CAiwas the only random variable, the sample mean response was not significantly different from the deterministic solution. The dispersion of C, was less for the nonlinear cases (EA = 1) than for the linear cases (EA = 0) for the same values of the rate factor ( k)V/(qi). When the dispersion of C.,, for the cases in which E.,, = 1 was compared at the same conversion of CA for the cases in which EA = 0, the values of Tc /( CAi) for the nonlinear cases were slightly smaller.

(1) The actual count of the values of C,/( CA1) which exceeded the critical value of C,/(C,,) associated with the confidence level of 0.99 agreed closely with the 3 values (out of 300) that would be expected to be greater than the critical value. (2) When the null hypothesis was used to test the hypothesis that the frequency distribution of C,/( CAi) obtained using values of (k)V/(q) less than 2.5 was the same as a normal distribution by the Chi-Squared Goodness of Fit Test, the hypothesis was accepted for almost every case. Based on the actual sample distribution of C,/( CAi) the assumption of a normal distribution provided a conservative value of C,J(C,i) at a confidence level of 0.99 when the rate factor had a value of less than 3 and a slightly optimistic value of C,/ 5. Results for Model 2B. (C,$> when the rate factor had a value greater We do not have the space to provide a detailed than 3. analysis of Model 2B, but wish to merely point Figure 7 and similar figures in Berryman[3] out two significant results that appeared to be can be used to determine the confidence limits different than for Model 2A. First, the effect of on the concentration of the outlet reactant, and random inputs and random coefficients on the hence to calculate the overdesign factor for plug mean response and standard deviation of the flow reactors. In Fig. 7, for example, if the value response of a reactor with a second order reof (li)V/(q) is known to be 3.0, then CA/( CA{) action is substantially greater than their effect on should have an expected value of 0.05. If the the mean and standard deviation of the response assumption about C,/( CAi>being normally distriof a reactor with a first order reaction at comparbuted holds, C,/( CAi) = 0.069 at the upper confiable levels of conversion. Second, over-design dence limit (for a probability of O-99). The differfactors required for safe design must be larger for ence between the volume of the reactor at the Model 2B than for Model 2A. Berryman [3] proupper confidence level for (CA)/ ( CAi) and the vides a number of charts to facilitate design and volume of the reactor at the ensemble mean of comparison among the plug flow reactors. (C,>/(C,i) gives a reasonable estimate of the required overdesign factor at the selected probability level. As a specific numerical illustration, COMPARISON OF RESULTS FOR THE ISOTHERMAL WELL MIXED TANK WITH Figure 7 can be interpreted as follows. If k is the THOSE OF THE PLUG FLOW REACTOR only random variable, Yk = 0.05 and the desired We now compare certain of the results obconversion is 95 percent so that (C,)/(CAi) = 0.05,thevalueof (k)V/(qi) at (C,>/(C,i) = 0.05 tained from the two models. If the random cois 3.0. At the confidence level of 0.99 the value of efficients and random inputs had the same coefficient of variation, the sample mean response (k)V/(qi) = 3.33 for (C,)/(C,,) = 0.05. Hence, the required overdesign factor is 3*33/3*0 = l-1 1 of the plug flow reactor was more likely to diverge from the deterministic solution than was the or an increase of 11 percent in the length of the sample mean response of the well mixed reactor. reactor. If yk = 0.10, then the overdesign factor However, divergence of the sample mean reswould be 3.6613.0 = l-22. 1268

Reactors for process analysis and design-

ponse of the plug flow reactor was quite small for the analogous reactions of Models A and B. Fluctuations of CAi entering the well-mixed reactor were damped with the amount of damping depending on the frequency of change of the random variable CAL.For the well-mixed reactor sc,I(CAi) was less than ~~,,/(CAL) and ?c,/(CAi) was less than qc,,/( CAi). For the plug flow reactor the term damped needs to be carefully defined because sc,/(C,,) was less than s,~/(C.~~) but p~.,/(C.~i) was equal to Tc,4i/(C,i). In both types of mixing the range of C.4 in the reactor was always less than the range of CAi (except at the inlet where they were equal) for the cases in which CAiwas the only random variable. The effect of the dispersion of qi* on the response of the plug flow reactor was not adequately treated in this study because of the extensive computational time required, hence a good comparison of the effect of a random & on the responses of the plug flow reactor and the wellmixed reactor cannot be made. Nevertheless, the effect of dispersion and the frequency of change of 4: on the response should be about the same for either type of mixing because a random & affects the residence time of a plug flow reactor or a well-mixed reactor in about the same manner. High frequencies of change of 4: almost negate any effect that the dispersion of 4: might have on the response of either type of reactor. A random k had a much greater effect on the response of the plug flow reactor than a random k had on the response of the well-mixed reactor. Figure 8 shows how for the same dimensionless output concentration the coefficient of variation of the plug flow model is considerably higher than for the well mixed tank. COMPARISON

OF RESULTS OF OTHERS

WITH

THOSE

No work has been reported on the plug flow reactor so that the comments here refer only to Model 1. Some investigators have reported that the ensemble mean response of the stochastic differential equation for the well mixed tank using random coefficients is not equal, in general,

Q6

I

r

\ \

Yk = 0.10

\ \

---

Well - Mixed

-

Plug Flow

Reactor Reactor

\\

\

.o E E z s 0

\ \ 0.2-

\

2 A! 5 'W

\ \

QI-

\

E

6

\

\

a3-

\

I

01

\

0

I

Coefficient

0.1 of Voriotion

I

0.2

of CAN&>,

I

ycA,scAi,

Fig. 8. Comparison of dispersion of response for well-mixed and plug flow reactors when k was the only random variable.

to the deterministic response, as was approximately found to be true in this work. To understand this paradox it is necessary to keep in mind that the dispersion of the random parameters used in this study was limited to a range such that (a) the fluctuations in the entering flow stream were physically realizable, and (b) the uncertainty associated with the coefficients corresponded to that reasonably experienced with empirical data. In contrast, other investigators have, in general, used white noise to generate the random coefficients. If particular care is not exercised, the use of white noise can lead to unrealistic values of the dispersion in the random coefficients and inputs of a model of a real process, and thus lead to erroneous results. For example, King[81 used a stochastic Green’s function to analytically calculate the ensemble mean response and the ensemble autocorrelation function of CA(t) for the well mixed tank. The random flow rate 4: was separated into two parts, qi* = (4T > + 2 , and two types of noise were introduced for the Z,: component: (a) white noise, Z,(I) and (b) the auto-correlated random noise ZA(I). Because white noise was added to (qt ) to generate the random flow rate,

1269

J. E. BERRYMAN

and D. M. HIMMELBLAU

sT, the values of & could and did assume negative values in King’s work. Since the numerical value of 4T at t = Of determined the numerical value of the tracer concentration CA at t = 0+ in King’s work, it is apparent that the realization of C, (O+) included some negative values, and thus the procedure of King did not represent a model of a real stirred tank when 2 was white noise. Even if the white noise term had been replaced by a normal random variable with a finite variance (Y,the values of CX/(C,) = 1 and 1.5 used by King were so large that negative flow rates and negative values of CA(0+) would still have occured. Figure 9 shows the analytical results of King for different noise levels (Us) , and the simulated response obtained in this work (the open symbols) at the same noise levels. Note the good agreement between the results of King’s analytical study and the sample points obtained by simulation in this study for these two noise levels, but keep in mind that the dashed curves

/
do not represent a physically realizable wellmixed tank; only the solid curve does. Pell and Aris [9] presented a study of a wellmixed reactor with the reaction sequence A*B+C that included a random absolute temperature term in the expression for the reaction rate coefficients. Although Pell and At-is presented their results for the variance and mean response of CA and CB for (a) both the sample statistics of the simulation runs and (b) the ensemble statistics from their analytical calculations, as a function of the power spectrum of Z,(t), it was decided that a more meaningful presentation here would be to show the variance and mean response as a function of the coefficient of variation of T, yT. Figure 10 shows the good agreement between the values of the sample mean response obtained in this work and the ensemble values reported by Pell and Aris, although the tendency of C, (/3 = 10) to decrease as YTincreases is much more pronounced in this study. It is interesting to note that the values of both (C,) and C, for j3 = 40 are, in effect, an extension of the p = 10 data when plotted as a

> 0.34-

Results

_

3,75

0.50

0 ‘g\

-

Deterministic

Solution

CA/CO

,o”

‘\

om-

0 \ ‘\

IS 90

0.26 -

\

“b,

8 i 4 ‘= E .o 5

This \

Pal and Ads’s Work 0,24_

E $ s o 6

‘\

0.26-

--

p= IO

--

P=10 -.-

0,22--

#g=40 j3=40

0

F,

r”

I

Dimensionless

Time,

, 0,02

0.20 I 0

2

Coefficient

y’

Fig. 9. Comparison of results for the unsteady state response of a stirred tank to a dirac pulse input of tracer.

of Variotlon

I 0.04

1

of T, YT

Fig. 10. Mean concentrations of A and B as functions of the coefficient of variation of T.

1270

Reactors for process analysis and design-1

function of YT. When (C,) were plotted as representing (C,) for line representing (C,)

4:

Pell and Aris’s values of a function of YT, the line p = 40 coincided with the for /3 = 10.

-* Y R

RI, Rz, Rs MAIN

QA

CONCLUSIONS

In general the relative frequency distributions of the steady state values of C,/(C,i) appeared to be normal or approximately normal regardless of distribution used in a simulation trial for the random inputs and coefficients. The well-mixed tank acts as a “damper” relative to the plug flow reactor (which acts as an “amplifier”) to reduce the effect of dispersion of random coefficients and random inputs on the dispersion of the output stream regardless of whether a reaction is taking place in the vessel or not. The frequency of change of the time dependent random parameters affects the extent of the damping; a high frequency of change gives substantial damping whereas a low frequency of change results in less damping. The sample mean response or Models A and B agrees for all practical purposes with the deterministic solution in all the cases investigated, but the dispersion is in general greater for the nonlinear reactions.

T t

12 t* V V x

x(t) i Z CA

Greek symbols

variance a constant, a parameter in the autocorrelated random number generator that adjusts the frequency of change of the random number ensemble coefficient of variation of YZ any random variable Z sample coefficient of variation of any % random variable Z A difference o..4 ensemble standard deviation of A

NOTATION

a deterministic

constant in general of reactant A of reactant A at

molar concentration C::molar concentration

C*i

inlet expected value of X E* fractional change of system volume between no conversion and complete conversion GA conversion of reactant A k, kl, k, reaction rate coefficient L reactor length NA the number of moles of reactant A NAi the number of moles of reactant A at inlet flow rate 4 exit flow rate 9s flow rate at inlet qi

EL-U

flow rate at inlet in reduced time units time average flow rate (reduced time) reaction rate ratios given by graphs sample standard deviation of CA temperature time integration limit dimensionless time (time divided by the residence time) reactor volume velocity fraction of reactor volume a deterministic input variable a deterministic dependent variable random variable (C, -Y>

Subscripts A B ; N P 0

reactant A reactant B inlet condition outlet condition a normally distributed able product P initial condition

random vari-

Other symbols

1271

(2

X

ensemble average (overlay bar) sample mean

J. E. BERRYMAN

[l] [2] [3] [4] [5]

[6] [7] [8]

[9]

and D. M. HIMMELBLAU

REFERENCES ACRIVOS A., Chem. Engng Sci. 1960 12 279. ARIS R. and AMUNDSON N. R., Chem. Engng Sci. 1958 9 250. BERRYMAN J. E.,Ph.D. Dissertation, University of Texas, Austin, Texas 197 1. BERRYMAN, J. E. and HIMMELBLAU D. M., Ind. Engng Chem. Proc. Des. Deu. 1971 10441. BHARUCHA-REID A. T., Elements ofthe Theory ofMarkou Processes and Their Applications. McGraw-Hill, York 1960. HOMAN C. J. and TIERNEY J. W., Chem. Enpng Sci. 1960 12 153. KATZ S., Chem. EngngSci. 1958 961. KING R. P., Chem. Engng Sci. 1968 9 1035. PELL T. M. JR. and ARIS R., Ind. Engng Chem. Fund/s 1969 8 339.

1272

New