0097-4943/91 $3.00 + O.00 Pergamon Press plc

SOME APPLICATIONS OF COMPUTERIZED SYMBOLIC MANIPULATION IN THE ANALYSIS OF CHEMICAL ENGINEERING SYSTEMS WILLIAM R. THISSELL Department of Chemical Engineering, Washington University 1 Brookings Dr., Campus Box 1198, St. Louis, MO 63130-4899 PATRICK L. MILLS Du Pont Company, Central Research and Development Experimental Station, WiLmington, DE 19880-0262

(Received Februa~'y 1990)

Abstract--The potential advantages of using symbolic manipulation software in the analysis of Mathematical ModelS that describe typical chemical engineering systems is described. It is suggested that the systematic application of computerized symbolic manipulation would permit application problems in chemical engineering to be solved with greater efficiency and accuracy than using tedious manual operations which are susceptible to error. In addition, a greater degree of mathematical sophistication is possible which can provide greater insight into the physical situation when compared to purely numerical techniques. To introduce the application of symbolic manipulation to chemical engineering systems, several example problems are briefly described and solved using the software package MACSYMA. The particular ones used here describe the flow of thin liquid films, the exit age probability density function for nonideal mixing in a continuous-flow stirred tank, the steady-state reactant concentrations that exist in a train of backmixed reactors with first-order reaction, and the steady-state concentrations of an organic species in a countercurrent flow multistage extraction column. The use of these for pedagogical purposes and as an entry point for analysis of more sophisticated applications is suggested and discussed.

NOMENCLATURE

B CAi C1 c~ ~O~t 7Your(°) -j E(O 9~ k~ L~ M P

= matrix, first appearing in eqn. (25). ---- column vector, first appearing in eqn. (25). ----- matrix, first appearing in eqn. (19). = column vector, first appearing in eqn. (19). = concentration of species A in reactor i, mol cm -3. -- concentration of tracer in the first stirred tank, first appearing in eqn. (ga), reel cm -3. = concentration of tracer in the second stirred tank, first appearing in eqn. (ga), reel cm-3. = concentration of tracer in the exit stream, defined by eqn. (11), reel cm -s. = Laplace transform of C l ( t ) , defined lay eqn. (12), reel cm - a s. = Laplace transform of C2(t), reel cm -a s. = Laplace transform of Gout(t), defined by eqn. (17), mol cm -3 s. = column vector, first appearing in eqn (19). = exit-age probability density function of the tracer, defined by eqn. (8), s-1 . = y component of gravity, first appearing in eqn. (1), cm s-2. = first-order rate constant, first appeaxing in eqn. (17), s-1 . = mass flow rate of water that exits tile i-th stage, first appearing in eqn. (24), g s-1 . = slope of the linear equilibrium relationship that relates the aniline concentration in toluene to the aniline concentration in water appearing in eqn. (26), g of water/g of toluene. = moles of tracer injected in the input pulse, first appearing in eqn. (8), reel. = parameter define by eqn. (16), rm12 s-2.

Financial support for the second author to complete this work f,~om Autoclave Engineers Group is appreciated.

Typeset by .AA,~-TEX 53

54

W.R. THISSELL, P.L. MILLS

Q

= volumetric flow rate of the film, first appearing in eqn. (5), cm3 s-1 .

Qo

= volumetric flow rate of the incoming fluid, cmz s -1 .

Q~

-- volumetric flow rate of the exit fluid, ornz s-1 .

Q~

= volumetric flow rate of the incoming fluid to the i-th reactor, first appearing in eqn. (17),

cmz s-1.

8

---- Laplace transform variable, first appearing in eqn. (12), s-1.

t

= denotes time, s.

1)

-- film velocity, first appearing in eqn. (1), cm s-1 .

v~

= volume in the first well-mixed tank, first appearing in eqn. (9a), cm3.

Ѕ

= volume in the second welbmixed tank, first appearing in eqn. (9b), cm3.

= mass flow rate of toluene that exits the i-th stage, first appearing in eqn. (24), g s-1 .

w

= plate width in the z coordinate, first appearing in eqn. (5), cm.

:g

= denotes the x coordinate, cm.

xi

= concentration of aniline in the water phase on stage i, first appearing in eqn. (24), g

anifine/g water.

= column vector containing the xi, first appearing in eqn. (25).

Y

= denotes the y coordinate, cm.

Yi

= concentration of aniline in the toluene phase on stage i, first appearing in eqn. (24), g

aniline/g toluene.

Greek.

ot

:- fraction used to define the volumetric flow rate that enters the first well-mixed tank,

dimensionless.

= angle measured relative to the vertical, radiaus.

6

= film thickness first appearing in eqn. (3b), cm.

/t

= fluid viscosity, gcm -1 s-1 .

P

= fluid density, gcm -3.

INTRODUCTION The use of computers for performing symbolic manipulation was first introduced by Kahrimartian [1] and Nolan [2] in the early 1950's by i m p l e m e n t i n g various rules for differentiation of functions in closed-form. Advances in both Computer Science and computer technology, especially in the past decade, have resulted in the emergence of symbolic manipulation as a significant area for research and applications in mathematics, science, and engineering [3]. These have been described and summarized in the proceedings of various aspects of the subject. An important consequence of the above research has been the development of symbolic manipulation p r o g r a m s or so-called c o m p u t e r algebra software packages [3] which can execute a variety of mathematical operations associated with algebra, trigonometry, calculus, differential equations, vector and tensor analysis, and other related analytical techniques. These programs relieve the user of the traditional labor intensive and error-prone manual methods of mathematical analysis. In many cases, these programs provide other attractive benefits, such as computer graphics display of results, automatic generation of high-level language computer code, and implementation of numerical analysis techniques [3]. Examples of these programs include A L T R A N [18,19], C A M A L [20], F E S T E R [21], M A C S Y M A [16,22-26], M A P L E [27,28], M A T I t E M A T I C A [29], muMAWtt [30-32], REDUCE [33,34], SCRATCItPAD [35], SIMATtt [36], and SMP [37,38]. This represents a small sampling of more than sixty packages which were reviewed by Yun and Stoutemeyer [39]. Applications of symbolic manipulation programs are quite numerous and range from pure mathematics to more practical engineering type problems. Descriptions for some of these are scattered throughout the references cited above which also cite numerous application-oriented review papers. Excellent introductions to the subject are provided by Rand [16], Drinkard [22], I-Iussain and Noble [26], Calmet and van ttulzen [40], Noor and Andersen [41], Roach and Steinberg [42], and Pavelle [43,44]. Applications involving the use of MACSYMA, which is often cited as one of the most powerful and widely used computer algebra packages [45], are summarized in a bibliography [46] that contains several hundred literature citations.

Appllcations of computerized symbolic manipulation

55

An examination of the above literature reveals that the number of chemical engineering applications where MACSYMA and other computer algebra programs have been utilized is less than ten. Considering that fundamental aspects of this discipline are based upon various conservation laws and constitutive relations whose general forms are nonlinear partial differential equations, this would appear to be a surprising observation. Closer examination of the history of chemical engineering [47] suggests, however, that systematic use of advanced mathematical techniques to solve chemical engineering problems was not introduced until the mid 1950's and early 1960's with the pioneering work of Amundson and Aris [48]. Notable advances have occurred between then and the present such that a variety of sophisticated analytical and numerical techniques are routinely utilized to obtain solutions to mathematical models of chemical engineering systems. Use of computer algebra software to support the mathematical analysis of these models in both an academic and industrial setting would permit more detailed and rigorous investigations to be performed with obvious benefits. Incorporation of computer algebra into chemical engineering education on both an undergraduate and graduate level, such as in the teaching of transport phenomena, chemical reaction engineering, process dynamics and control, to name a few fundamental areas, would provide students with greater exposure to more realistic problem situations as one benefit. One objective of this paper is to demonstrate the utility of using a computer algebra package for solution of mathematical models that describe the performance of selected chemical engineering systems. It is surmised that if this approach was more widely advocated in education and these packages were readily available, use of computer algebra in mathematical analysis of chemical engineering systems would be more widespread. While anyone with sufficient stamina who is armed with a working knowledge of the mathematical techniques used and a computer could obtain the results reported here, the use of the computer algebra package represents a more general and modern approach to problem solving while providing benefits, such as greater efficiency and accuracy, when compared to the traditional approach. As might be expected, some tradeoffs and unexpected problems can occur when using a computer algebra package and these will be pointed out where applicable. The particular computer algebra package selected for this study was MACSYMA since it is widely available and was utilized in an earlier paper in this journal to analyze the dynamics of helicopter rotor blades [49]. It is expected that identical results could be obtained with one or more of the software packages cited above so the use of MACSYMA should not be viewed as a limitation.

USING COMPUTER-AIDED SYMBOLIC MANIPULATION AND MACSYMA Using computers to perform symbolic manipulation from both a theoretical and practical perspective has been the subject of much research. The recent monographs by Davenport et al. [3] and Rand [16] provide the necessary background information as well as other important details on how to use the MACSYMA package so these will be omitted here for brevity. A few comments about the basic usage are given to serve as an introduction to the example problems which follow. The MACSYMA user interface is structured into command lines and display lines which have the form ўi and di where i is an integer. Each command line ўi has a corresponding display line di where the former is input by the user and the latter is the MACSYMA response to a syntactically correct command line. For example, the mathematical function f ( x ) = 2x + sinhx would be entered in MACSYMA as ( c l ) f ( x ) ffi 2*x + s~nh (x); where the semicolon represents the command line terminator. To obtain the indefinite integral of f(x), one would enter a = integrate (f(x), x, 0, x), and b = diff (f(x), x,l)to obtain the first derivative of f(x) with respect to x. These and other commands are illustrated in the examples given below.

EXAMPLE MACSYMA PROBLEMS Several example problems which illustrate the use of MACSYMA to obtain solutions to typical chemical engineering systems are now introduced. The particular examples given here were selected, in terms of difficulty, to be relatively easy to more complex. Research type problems were intentionally omitted since their treatment lies outside the objectives of this work. Here, the goal is to provide some perspective of how the package can be used on various types of

56

W.R. THISS~.LL,P.L. MILLS

problems involving differential or difference equation models for selected chemical engineering systems. One possible long term goal is to maintain a library of these various computer algebra applications for general reference that users can share.

Example 1. Film Flow Down an lnclined Plane The hydrodynamics of thin liquid films has been widely studied by mathematicians and engineers in various contexts and leads to a number of challenging problems that have practical applications. The most fundamental problem of this type which is introduced in most chemical engineering texts on transport phenomena [e.g., 50] involves determining the velocity profile for flow of a viscous, isothermal liquid film along an inclined plate. A diagram of the system under consideration is presented in Fig. 2.2-2 of [50]. A cartesian coordinate system is chosen so that the origin of the (x, y, z) axes is at the top of the plate where x denotes the film thickness normal to the plate, y denotes the distance which the film has travelled down the plate, and z denotes the film width. The velocity vector "~ is assumed to have velocity components (u, v, w) corresponding to the x, y, and z directions, respectively. The inclined plane is assumed to make an angle ]~ relative to the vertical where ~ = 0. By making the usual assumption that the velocity component in the y direction v is a function of tim x and y coordinates with the remaining velocity components u and w being zero, the foUowing reduced form for the y momentum balance can be derived for constant fluid density p and viscosity p [50, Table 3.4-2]

P = \ ox2 + ou2 ) + pay.

(1)

The equation of continuity gives [50, Table 3.4-1]

COy -c-Oy ~ 0 .

(2)

The boundary conditions are no slip at the plate surface and zero shear stress at the film surface

v=0

atx=0,

(3a)

COy

/J ~xx = 0

at z = 6,

(35)

where the parameter ~f is the film thickness. Substitution of (2) into (1) gives the following second-order ordinary differential equation with constant coefficients whose solution is the velocity component v as a function of film thickness z

d2t)

+ Pgv = 0.

(4)

dx 2 I.t

Once the velocity profile v(x) is determined, the volumetric flow rate Q of the film per unit width W can be determined by integration of v(x) over the film cross-sectional area

wq- wif0' v(x)dx.

(5)

The manual solution of equation (4) is straightforward and requires two simple integrations whose constants of integration are determined from boundary conditions (3a) and (3b). Performing these operations and substituting for the y component of the gravity force gy=g cos(/~) yields the following expression for v(x)

,(,,,) pg6cos(#) pg cos(#) =2.

=

/.=

=

2/~

(6)

Applications of computerizedsymbolicmanipulation

57

Substitution of (6) into (5) and carrying out the integration yields the following expression for Q/w

Q p9 63 cos( )

=

3p

(7)

The computer algebra solution of equation (4) using MACSYMA is summarized in APPENDIX 1 which also serves as an introduction to the conversational mode of problem solving. A key step in the procedure involves using the MACSYMA package ODE2 to solve the second-order ordinary differential equation defined by equation (4). Inspection of the final results given by (dl0) and (rill) in APPENDIX 1 shows the total CPU time requirement was 10.9 seconds which is the typical order-of-magnitude for relatively easy problems of this type. The extension of the given procedure to solve more complicated problems in transport phenomena involves systematic application of these and other advanced MACSYMA commands.

Ezample 2. Ezit-Age Distribution for a Tracer Flow Model The use of tracer flow models to characterize the performance of chemical reactors and other related contactors has been widely used in chemical engineering and other applied sciences. These are introduced in nearly all texts on chemical reaction engineering [e.g., 51,52] and have been the subject of several monographs [53,54]. A recent review by Dudukovi6 [55] provides an excellent overview of this area. The theoretical portion of tracer methods typically involves mathematical models based upon ordinary or partial differential equations, integro-differential equations, or differential-difference equations. Application of symbolic manipulation software in the mathematical analysis of these models would appear to alleviate the manual labor associated with both analytical and numerical techniques. In this example, application of MACSYMA to solve the differential mass balance equations which occur in a tracer flow model for a nonideal continuous-flow stirred tank reactor is illustrated. A diagram of the flow model is given in Fig. 1 which consists of two well-mixed vessels whose volumes are Va and Vd, respectively. Fluid enters the system at a total volumetric flow rate of Qa where a fraction aQa enters the first vessel and achieves an instantaneous concentration Cl(t). The remaining fraction of the feed (1 - a ) Qa bypasses the reactor and is blended with an effluent stream whose volumetric flow rate is a Qa. The second tank receives fluid from the first tank at an inlet concentration Cl(t) and discharges fluid at an exit concentration C2(t) where the inlet and exit flow rate is constant at Qe. At time t = 0, a pulse of tracer M6(t), where M is the moles or mass of tracer, is injected in the flowing inlet stream. The quantity of interest is the exit-age density function of the tracer measured at the reactor outlet E(t). It is defined by [53,54]

E(t) = Co,,,(t)

M/Q. "

(8)

The significance of E(t) (it is that it represents the fraction of molecules in the reactor effluent which have residence times between t and t + dr. A knowledge of E(t) can be used to predict nonideal reactor performance for the case where first-order reactions are occurring and other troubleshooting applications which are summarized by Dudukovi6 [55]. To evaluate equation (8), a closed-form expression for Coup(t) is necessary. The tracer material balance equations for the system shown in Fig. 1 lead to the following coupled system of ordinary differential equations:

dC1

Va

= Q. C2 - C1 -

C1,

de2

Vd

= Q. Cl - Q. C2.

The initial conditions are

t = o,

= -aPMI '

t = O, C'~ = O.

(ga) (gb) (1oa) (lOb)

58

W . R . THISSELL, P . L . MILLS

M 5(0

aQ a

(1- a) Qa

r ::::::::::::::::::::::::.:::::::::: °ўQa

Cl Qe C2 :;iiiuiiiiiiiiiiiiiii Qe

Figure 1. Flow model for a nonideal continuous flow stirred tank reactor.

The tracer concentration at the system exit Coup(t) is related to the reactor exit concentration Cl(t) by a material balance around the intersection of the streams having volumetric flow rates o f ( 1 - a ) Q a , aQa, and Qa:

Cout(t) -- (1 - a) M 6(t) + . c l ( t ) .

(11)

One straightforward method for solution of (9a)-(ll) is to use Laplace transforms. The Laplace transform of Cl(t) is

-Cl(s) = Cl(t) e - " dr.

(12)

Analogous expressions for C2(t) and "Cout(s) follow directly from (12). Application of (12) to (9a)-(10b) results in two linear algebraic equations whose unknowns are ~l(S) and i22(s). Since Cl(t) appears in (11), we solve for Ul(s) by elimination of U2(s) and substitute the resulting expression for Cl(s) into the Laplace transform of (11) corresponding to "Coup(s). Inversion of Cout(s) into the time domain will yield an explicit expression for Cout(t) and hence the exit-age density function E(t) defined earlier by (8). APPENDIX 2 gives the details associated with performing the above operations using MACSYMA. The key commands are LAPLACE ( e x p r e s s i o n , t r a n s f o r m e d v a r i a b l e ) and ILT (exprusion, transformed variable, transform variable) which perform the Laplace transform and the inverse Laplace transform operations, respectively. The MACSYMA output at (d16) in APPENDIX 2 shows that the Laplace transformed outlet concentration of tracer is given by

c~(a M s Va + a M Qe)

(1-a)M

-Cout(s)= [s2Va+(Qe+aQa)s]Vd+QesV~+aQaQe + Q~

(13)

When the ILT (inverse Laplace transform) command is applied to (13), MACSYMA queries the user to determine if the following group of terms is positive, negative or zero

= 4aQa Q~ va Vd -- [(Q~ q" a Qa) Vd + Q~ V~]2.

(14)

It can be readily shown by inspection that (14) will always be negative for any situation of practical interest. For completeness, the inverse Laplace transform of ~ou~ (s) is determined for

Applications o( computerized symbolic manipulation

59

all three cases in (d17), (d20), and (d22) of APPENDIX 2 to illustrate the ease associated with determining the explicit expression of Co,z(t). The final expression for E(t) can be obtained directly from (all7):

E(t)= (2a2QoQ, V, Va-a~Q, Vd[(Qe+aQo)Va+QeVo])

x/-fi

+a' Oo V,zcoshk2Vo Vd] j V--~+(1-a)6(t). (15)

The parameter P which appears in (15) is defined by

P - (Qe + aQa) 2 V~ + 2Qe (Qe - aQa) V° Vd + Q,2 V~2.

(16)

The above flow model of a non-ideal stirred tank is just one among many that can be systematically analyzed using a symbolic manipulation program. This type of analysis provides useful insight into diagnosing the performance of real systems and is commonly used in many practical applications. Other tracer probability density functions, such as the cumulative age probability-density function F(t) and the intensity probability function A(t), could be readily derived from the above expression for E(t). Definitions for these functions in terms of E(t) and their significance is given by Dudukovi6 [55]

Example 3. Steady-State Performance of a Stirred Reactor Train Chemical engineering systems often consist of stagewise operations where various fluid streams at a given set of process conditions exit or enter a given control volume where they are contacted with other fluid streams and undergo physico-chemical changes. Some examples of these systems include distillation colunms, absorption towers, liquid-liquid extractors, liquid-solid leaching systems and stagewise reactor trains. Mathematical models for these various fluid-contactorstypically lead to difference-typeequations which can often be cast into matrix form [56]. In this example, M A C S Y M A is used to solve the linearsystem of algebraicequations which describe the product concentrations in a process system containing multiple continuous stirred tank reactors (CSTR's) connected in serieswith recyclestreams. A diagram of the system is given in Fig. 2 which consistsof four CSTR's of unequal volume. The fluidin each reactor is assumed to be well-mixed and contains species A that undergoes an irreversiblefirst-orderreaction. Each reactor is maintained at a constant temperature, but the temperature varies from one reactor to the next. The reaction rate is defined on a liquid volume basis and given by R / = kiVi where k is the first-orderrate constant, V is the reaction volume, and the subscript i can assume values of/= 1,2,3 or 4 corresponding to the reactor number. The volumetric flow rate into the i-th reactor from the (i- 1)-th reactor is Qi,ln,while the volumetric flow rate that exits the i-th reactor is Qi,out.If a portion of the effluentfrom reactor i + 1 is recycled back to the previous one, then the volumetric flow rate of this stream is Q~+1,r. The concentration of species A in the i-th reactor is CAi which is fixed at a known value CAo in the system feed. A material balance for conservation of species A around the i-th reactor at steady state gives the following differenceequation:

Qi,n CA,i-1 "~"Qi+l,r CA,i+l -- Qi,ou, CA,i -" ]gi CAi Yi.

(17)

A material balance for conservation of total mass around each reactor at steady state relates the inlet, exit and recycle stream volumetric flow rates and results in the following equation where constant fluid density is assumed:

Qi,ou, -- Qi+l,i. q" Qi.r.

(18)

Application of (17) and (18) to the reactor train shown in Fig. 2 gives a linear system of algebraic equations which can be cast into standard matrix form

e - d.

(19)

60 Q2, r l CA2 Q i, in

W . R . THISSELL, P.L. MILLS

l

~

Q4._.2r

c[ ·

k4 CA2, J[ ~ - -| k~ 3 ў:~1 CA4 I

Q1, out = Q2, in

Q2, out = Q3, in

Q3, out = Q4, in Q4, out

QI, in

CA4 Figure 2. Sequence of four continuous flow stirred tank reactors with recycle streams. P a r a m e t e r values: I,~ = 1000, V2 = 1500, V3 = 100, V4 = 500, kl = 0.1, k2 = 0.2, kz = 0.4, k4 = 0.3; Q l , i n = 1000, CAo = 1, Q2,r = o, Q z , r = Q4,r = 100.

m The matrix B and column vectors ~ and d are

[Ql,ou~q-kl V1 -Ql,r

0

0

1

[o ~ = I - Q 2 , , ,

Q2,,u,+k~V2 -Q3,r

0

-Q3,i.

Q3,o., + ks Va -Q4,r

] (20)

0

-Q4,1,,

Q4,o,,, + k4 V3

-~=[CA1 CA2 CA3 CA4]T,

(21)

----[ql,inCAo 0 0 0]T .

(22)

The tridiagonal matrix defined by (20) commonly occurs for mathematical models for staged chemical engineering systems and other applications, such as finite-difference approximations to differential equations. Although the number of unknowns for this example is small when compared to that encountered in more sophisticated applications, it provides a basis for illustrating the use of MACSYMA for solutions of linear systems of equations. APPENDIX 3 gives the details of solving the linear system (19) for the unknown species concentrations using MACSYMA. One key command used in (ў2) and (e3) in APPENDIX 3 is ENTEIOIhTRIX (rous, columns) which permits the user to input the numerical values assigned

to the array elements of B and d. Another important command is INVF.JtTwhich is used in (c4)

to perform a matrix inversion of B from which the unknown vector of species concentrations

is determined using ~ - B MACSYMA in (c4) is

d. For the process parameters used here, the vector ~ given by

10 28s00 2500 220o]T

~ = H 4132r 3757 3--~j

(23)

Appllcatlons of computerizedsymbolicmanlpulation

61

It is interesting to note that the solutions produced by MACSYMA are integers which represent the exact results. Use of standard methods for solving linear systems of equations give floating point results whose accuracy is limited by the computer precision and rounding errors. The linear system of equations defined by (19) was also solved by using FORTRAN subroutines from the IMSL library for comparison. The particular ones used include DLFCRG and DLFSRG which implement LU decomposition, DLSARB which solves banded linear systems, and a special purpose written subroutine that implements the Gauss-Seidel method. All three produced results that agreed within single precision accuracy (6 to 7 digits) with those given above in (23) obtained from MACSYMA. The effort required to write and debug the calling programs was considerably more when compared to that required to obtain a MACSYMA solution, however. This suggests that a good working knowlegde of symbolic manipulation software can lead to more efficient use of time for obtaining results to certain classes of problems.

Ezample 4. Steady-State Performance of an Aniline Eztraction Column

This particular example is given as problem 3-2 in the text by Constantinides [57]. Here,

it was selected to illustrate some possible disadvantages of using MACSYMA to obtain purely

numerical results. A diagram of the system is given in Fig. 3 which is s ten-stage countercurrent

tower for removal of aniline from water by solvent extraction. The water-aniline mixture, which

is the heavy continuous phase, enters the top of the column and flows downward. Fresh toluene

solvent, which is the light dispersed phase, enters the bottom of the column and flows upward

where it contacts the water-aniline mixture. Due to the greater solubility of aniline in toluene

when compared to water in toluene, the aniline is extracted from the aniline-water mixture

which results in a bottom product that contains predominantly water. The general problem is to

determine the concentrations in both the aqueous and organic phases on each stage of the tower

as well as at the tower exit. Additional details on extraction processes are available in standard

texts [e.g., 56] on the subject.

A material balance for aniline around the i-th stage in Fig. 4 at steady-state gives the following

difference equation:

Li zi + V~Yi - Li-1 z i - i -F ~+I Yi+1 = 0.

(24)

In the above equation, Li is the mass flow rate of water that exits the i-th stage, ~ is the mass flow rate of toluene that exits the i-th stage, zi is the concentration of aniline in the water phase in mass units on a aniline-free basis, and Yi is the concentration of aniline in the toluene phase where the toluene is defined on a aniline-frce basis. The stages are numbered from the top down so that definitions for the same variables with subscripts i - 1 and i-I- 1 is straightforward and omitted for brevity. Recycled toluene solvent is introduced above stage 7 so that the term F z! is added to the right-hand side of (24) for i = 6 where F and z! have the same units as Li and zi, respectively. If a linear equilibrium relationship is assumed, then yi = m zi so that expansion of (24) for i = 1, 2 , . . . , n gives a linear system of algebraic equations of the form

A ~= b,

(25)

where

[-ўL+, Irav1) my2 m~ -(~mv2)-(Ls+mvs) mV~

~-----

"'" L,-I -(Li'Jr'mV~) mV~+1

o

o ·

"'"

L',:_'2 -(L.-,'+'.,V._,)

Ln-*

] (26) ,,,V. -(L. + m Vn)

(27)

b=[L0z0 0 0 0 0 Y z! 0 0 0 0]T.

(28)

62

W.R. THISSELL, P.L. MILLS

The matrix A in (26) is given for an arbitrary number of stages for generality, but n = 10 is used here. T h e variable assignments are Li = 100 for i = 0, 1 , . . . , 10, 1,] = 23 for i = 0, 1 , . . . , 6 and = 10 for i = 7,8,... ,10, m = 9, x0 = 0.05 and Xl = 0.003.

Water Phase With Aniline L0 = 100 kg/hr water x0- 5kg/hranlline

~

[

Extract ~ Aniline-richToluene ~ V0 ,, 23 kg/hr toluene

·

yo=?

1

2

3

4

5.

Recycled Solvent

6

F =13 kg/hr toluene xf = 0.003 kg aniline/kg toluene

7

.8

9

Raffinate L10 = 100 kg/hr water Xl0 = ? Aniline-leanwater

10 Fresh Solvent V l l = 10 kg/hr toluene Yll = 1

Figure 3. Countercurrent flow liquid-liquid extraction column with ten ideal stages. Parameter values: Li -- 100, for i = 0,1 .... ,10; ~ -- 23 for i = 0,1 ..... 6 and I f / = 10 for i = 7 , 8 , . . . ,11; F -- 13; x0 = 0.05, xy = 0.003, yll = 1.

0.0

t

i

t

t

t

i

t

t

t

0.0

-0.5 --

-- -0.5

-1.0 -- -1.5 -- o- N -2.0 - ў.0 o -2.5 -- 1 -3.0 --

-1.0 -1.5 -2.0 ~ M -2.5 -3.0

-3.5 --

-3.5

-4.0 0

-4.0

1

2

3

4

5

6

7

8

9

10

STAGE NUMBER, N

Figure 4. Calculated stage compositions of aldline in the water (xi) and toluene (yi) phases. Parameter values are given on Fig. 3.

Applications of computerizedsymbollcmanipulation

63

APPENDIX 4 gives the M A C S Y M A solution to (25) for the unknown aniline stage concentra, tions in the water phase. The elements o f ~ and b are entered using the EFrERMATItIX command. Since ~ is a 10 x 10 matrix, all 100 elements must be manually entered which is a tedious task and represents one drawback of using this approach. A user written program would reduce the effort by accounting for the banded structure of the matrix and eliminate the need to input the zero elements. Additional inspection of APPENDIX 4 at (c7) shows that C P U time required to solve for · by matrix inversion and multiplication is 168 seconds. The same operation using the I M S L subroutine in F O R T R A N required less than 1 second so that the algorithm used in M A C S Y M A is inefficient. Use of M A C S Y M A as a device for extensive numerical calculations is not recommended and alternate approaches based upon conventional approaches should be used. Figure 4 shows the aniline concentrations in the water and toluene phases on each tray which follow the expected behavior for continuous countercurrent liquid-liquidextraction columns (Smith, 1963). The next extension of this problem involves the use of an equilibrium relationship of the form Yi = a + b zi where a and b are known constants. The system of equations defined above by (25) can be rewritten so that the column vector b is nonlinear in the zi. Successive substitution can be used to solve the system for an assumed initialestimate of the z~. This problem can be readily solved by M A C S Y M A also, although conventional F O R T R A N - b a s e d software requires two orders-of-magnitude less C P U time.

SUMMARY AND CONCLUDING REMARKS The use of symbolic manipulation software to obtain both analytical and numerical solutions to mathematical models for chemical engineering systems has been discussed and introduced through several example problems using MACSYMA. The emphasis here has been placed upon demonstrating the translation of the formal mathematical equations into a form which can be recognized by MACSYMA, and the associated conversational nature of using the software to obtain intermediate and final results. The particular examples selected here represent basic problems in fluid mechanics, reaction engineering, and staged operations where the solution of differential equations and matrices are involved. The solution of problems having a greater complexity can be attempted using advanced commands and other features of MACSYMA or other packages once the user has accumulated some additional experience. Examples in the literature where this has been performed on chemical engineering research problems are lacking, and tedious manual methods are still widely used. MACSYMA and possibly other symbolic manipulation software does not represent a replacement for most of the popular numerical software libraries since it lacks the necessary sophistication and efficiency in certain instances. In the last example, the solution of a 10 x 10 system of linear equations required significantly greater CPU time than a conventional numerical analysis approach based on Gaussian elimination. Development of integrated software which utilizes the latest technology in algorithms and computing hardware for both analytical and numerical application represents one of many challenges in this field.

REFERENCES

1. H.G. Kahrimanlan,Analptic Differentiation by a Digital Computer, M.A. Thesis, Temple Univ., Philadelphia, Pennsylvania,(May 1953). 2. J. Nolan, Analytic Differentiation on a Digital Computer, M.A.Thesis,M.I.T., Cambridge,MMsachusetts, (May 1953). 3. J.H. Davenport, Y. Siret and E. Tourmer, Competer Algebra. Systems and AIgorithme ]or Algebraic Computation, AcademicPress, London,(1988). 4. E.W. Ng (editor), ProceedingJ of the 1979 E*tropcan Symposiem on Symbolic and Algebraic Compwtation, Springer Lecture Notesin ComputerScience72, Springer-Verhg, Berlin,Heidelberg,New York, (1979). 5. J. Calmet (editor), Proceedings of tire 198~ E~ropcan Computer Algebra Congre,s, Springer Lecture Notes in Computer Science144, Springer-Verlag,Berlin,Heidelberg,New York, (1982). 6. J.A. van Hu]zen(editor), Proceedings oy the 1983 Eeropcan Comp,tLerAlgebra Congress, Springer Lecture Notes in Computer Science162, Sprlnger-Verlag,Berlin,Heidelberg,New York, (1983). 7. J. Fitch (editor), Proceedings o] the I984 European Computer Algebra Congreu, Springer Lecture Notes in Computer Science174, Springer-Verlag,Berlin, Heidelberg,New York,Tokyo,(1984).

C~IMA 21:10-E

64

W.R. THISSELL, P.L. MILLS

8. B. Buchberger and B.F. Cavine=m (editors), Proceedings of the 1985 European Conference on Symbolic and Algebraic Computation, Springer Lecture Notes in Computer Science 203 and 204, Springer-Verlag, Berlin, Heidelberg, New York, Tokyo, (1985). 9. N. Inada and T. Soma (editors), Proceedings of the fnd RIKEN Symposium on Symbolic and Algebraic Computation, World Scientific, New York, (1985). 10. J.H. Davenport (editor), Proceedings of the 1987 Eur6pean Conference on Symbolic and Algebraic Computation, Springer Lecture Notes in Computer Science 378, Springer-Verlag, Berlin, Heidelberg, New York, (1987). 11. Anonymous, Proceedings of the Second Sipmposium on Symbolic and Algebraic Manipulation, ACM Inc., New York, (1971). 12. Anonymous, Proceedings of the 1976 ACM Symposium on Symbolic and Algebraic Computations, ACM Inc., New York, (1976). 13. Anonymous, Proceedings of the 1981 ACM Symposium on Symbolic and Algebraic Computations, ACM Inc., New York, (1981). 14. Anonymous, Proceedings of the 1986 A CM Symposium on Symbolic and Algebraic Manipulation, ACM Inc., New York, (1986). 15. D.E. Knuth, The Art of Computer Programming, VoL II, Semi-numerical Algorithms, Addison-Wesley, 2nd. Edn., New York, (1981). 16. R.H. Rand, Computer Algebra in Applied Mathematics: An Introduction to MA CSYMA., Pitman, London, (1984). 17. D.S. Arnon and B. Buchberger (editors), Journal of Symbolic Computation, Academic Press, New York, (1983). 18. W.S. Brown, A L T R A N User's Manual, 4th. Edition, Bell Laboratories, Murray Hill, New Jersey, (1977). 19. A.D. Hall, Jr., The ALTRAN system for rational function manipulation - A survey, Commun. ACM 14, 517-521 (1971). 20. J.P. Fitch, CAMAL User's Manual, University of Cambridge. Computer Laboratory, Cambridge, England, (1974). 21. J.P. Fitch and R. Hall, Symbolic Computation and the Finite Element Method, Proceedings of the 1987 European Conference on Symbolic and Algebraic Computation (J.H. Davenport, editor). Springer-Verlag Lecture Notes in Computer Science 378, Springer-Verlag, Berlin, Heidelberg, New York, pp. 95-96, (1987). 22. R.D. Drinkard, Jr., MACSYMA: A Program for Computer Algebraic Manipulation (Demonstrations and Analysis), Technical Document 6401, Naval Underwater Systems Center, (1981). 23. V.E.Golden and the Matldab Group, Introductory MACSYMA Documentation: A Collection of Papers, Mass. Inst. Technology, Laboratory for Computer Science, Cambridge, Mass., (1983). 24. The Matldab Group, MACSYMA Reference Manual, Vols. I and II, Version 10, MIT. Laboratory for Computer Science, (1983). 25. Anonymous, A Brief Overview of MACSYMA, The Mathlab Group, Laboratory for Computer Science, MIT, (1983). 26. M.A. Hussain and B. Noble, Applications of MA CSYMA to Calculations in Applied Mathematics, Report No. 83CRD054, Technical Information Series, General Electric Co., (1983). 27. R.W. Char, K.O. Geddes, W.M. Gentleman and G.H. Gonnet, The Design of MAPLE: A Compact, Portable and Powerful Computer Algebra System, Proceedings of the 1983 European Conference on Symbolic and Algebraic Computation (J.A. van Hulzen, editor). Springer Lecture Notes in Computer Science 162, Springer-Verlag, Berlin, Heidelberg, New York, pp. 101-115, (1983). 28. H.W. Char, K.O. Geddes, G.H. Gonnet and S.M. Watt, MAPLE User's Guide, Watcom, Waterloo, (1985). 29. S.P. Wolfram, Mathcmatica--A System for Doing Mathematics by Computer, Addison-Wesley, New York, (1988). 30. A. Rich and D.R. Stoutemeyer, Capabilities of the muMATH-79 Computer Algebra System for the INTEL-8080 microprocessor, Proceedings of the 1979 European Conference on Symbolic and Algebraic Computation, (E.W. Ng, editor). Springer-Verlag Lecture Notes in Computer Science 72, Springer-Verlag, Berlin, Heidelberg, New York, pp. 241-248, (1979). 31. H.S. Wilf, The disk with the college education, Amer. Math. Monthly 89, 4-8 (1982). 32. The Soft Warehouse Group, The muMATH-83 Symbolic Mathematics System Reference Manual, Soft Warehouse, Honolulu, Hawaii, (1983). 33. J.P. Fitch, Solving algebraic problems with REDUCE, J. Symbolic Comp. 1,211-227 (1985). 34. A.C. Hearn, REDUCE-$ User's Mannal Version 3.3, Hand Corporation Publication CP78 (7/87), (1987). 35. J.H. Griesmen, R.D. Jenks and D.Y.Y. Yun, SCRATCHPAD User's Manual, IBM Research Publication HA70 (June 1975). 36. R. BSffgen and M.A. Reichert, The Computer Algebra System SIMATH, Proceedings of the 1987 European Conference on Symbolic and Algebraic Computation (J.H. Davenport, editor). Springer-Verlag Lecture Notes in Computer Science 378, Springer-Verlag, Berlin, Heidelberg, New York, pp. 48-49, (1987). 37. S. Wolfram, Symbolic mathematical computation, Commun. ACM 28,390--394 (1985). 38. S. Wolfram, SMP Reference Manual, Inference Corp., Los Angeles, Calif., (1983). 39. D.Y.Y. Yun and D.R. Stoutemeyer, Symbolic Mathematical Computation, Encyclopedia of Computer Science and Technology (J. Belzer, A.G. Holzman, A. Kent, editors). Volmne 15 Supplement, Marcel Dekker, New York, Basel, (1980).

Applications of computerized symbolic manipulation

65

40. J. Calmet and J.A. van Hulzen, Compnter Algebra Applications, Computer Algebra: Symbolic and Algebraic Computation (B. Buchberger, G.E. Collins, R. Loos, editors). Sec~d Edition, SpringerVerlsg, Wlan, New York, pp. 245-248, (1983). 41. A.K. Noor and C.M. Andersen, Computarized symbolic manipulation in structural mechanics progress and potential, Comp,t. Str~et. 1O, 95-118 (1979). 42. P.J. Roach and S. Steinberg, Symbolic manipulation and computational fluid dynandce, AIAA J. 22, 1390--1394 (1984). 43. R. Pavelle, M. Rothstein and J. Fitch, Computer algebra, Sei. Amer. 245(6), 136-152 (1981). 44. R. Pavelle (editor), Applications of Computer Algebra, Khiwer-Nijhoff. to appear 45. N.J.A. Sloane, My friend MACSYMA, Notices, ATe/T Bell Labs., 40--43 (January 1986). 46. Anonymous, Bibliography of Papers Referencing MACSYMA, Symbolics, Inc., Cambridge, Mass., (March, 1985). 47. W.F. Furter (editor),Histor~ ol Chemical Engineering,Advances in Chemistry SeriesNo. 190,American Chemical Society, Washington, D.C., (1980). 48. N.R. Amundson, RA of Minnesota, Chem. Engng. Sci. 44(9), 1761-1762 (1989). 49. M.R.M. Crespo da Silva and D.H. Hodges, The role of computerized symbolic manipulation in rotorcraft dynamic anA]ysis, Comp. ~ Maths. toith Appls. 12A, 161-172 (1986). 50. R.B. Bird, W.E. Stewart and E.N. Lightfoot, Transport Phenomene, Wiley, New York, (1960). 51. O. Levenspiel, Chemical Reaction Engineering, Wiley, New York, Second Edition, (1972). 52. G.F. Froment and K.B. Bischoff, Chemical Reactor Analysis and [email protected], Wiley, New York, (1979). 53. C.Y. Wen and L.T. Fan, Models for Flow Systems and Chemical Reactors, Marcel-Dekker, New York, (197s). 54. E.B. Nauman and B.A. Buflham, Mizing in Continuous Flow Systems, Wiley, New York, (1983). 55. M.P. Dudukovid, Tracer Methods in Chemical Engineering: Technif,es and Applications, ProceecYmgs of the NATO Advanced Study Institute on Chemical Reactor Design and Technology (H.I. deLasa, editor). NATO Advanced Study Institute Series, Series E: Applied Series No. 110, Martinus Nijhoff, Dordrecht, The Netherlands, pp. 107-190, (1986). 56. B.D. Smith, Design of Efuilibri~m Stage Processes, McGraw-Hill, New York, (1963). 57. A. Constantlnides, Applied Nnmerieal Methods with Personal Compnters, McGraw-Hill, New York, (1987).

APPENDIX 1

MACSYMA Solution for Ezample 1 Access the MACSYMA program after receiving the UNIX prompt which is a ~, sign for this particular VAX computer. Y. macsyaa

This is UNIX MLCSYML Release 309. l .

(c) t976, 1984 Massachusetts Institute of Technology.

All Rights Reserved.

Enhancements (c) 1984 Symbolics, Inc. All Rights Reserved.

Type describe(trade_secret);

to see Trade Secret notice.

Type exec("man macsyaa"); for help.

Specify that all execution times will be displayed after a g~ven input command.

(ct) showtime: all;

Time= 0 msec. MACSYMA verifies that aU execution times will be evaluated.

(dl)

all

Define the y-momentum balance defined in the text by (4) where ~(x) denotes the velocity distribution function, d i f f ( v ( x ) , x , 2 ) denotes the second derivative of ~ with respect to x, r h o = p, g y = g , , and mu = ~,.

( c 2 ) 0 ffi r h o * gy / a u + d i f f ( v ( x ) , x , 2 ) ;

Time= 66 msec. MACSYMA evaluates the second derivative of v(x) and outputs the second-order ordinary differential equation.

(d2)

2 d 0 =--- (v(x)) 2 dx

gy rho + ...... mu

Solve the previous expression, denoted in MACSYMA by ~,, using the O.D.E. package ODE2 (EQ, Y, X) where EQ is a first or second order O.D.E. to be solved for Y as a function of X. Here ~, = EQ, v ( x ) = Y, and x = X.

66

W.R. THISSELL, P.L. MILLS

(c3) ode2(~, v(·), x); Batching the file /src/local/pkg/macsyma/share/ode2.mac MACSYMAsolv~ ~ e O.D.E. and ~ v ~ ~eexp~cit ~ I u t i o n f o r v ( z ) w i t h two ~ t u n d e t e r m i n e d c o n s t a n t s V ~ l and Y~2.

Batching done.

Totaltime= (d3)

33266 msec. GCtime= 7333 msec.

2 gy rho ·

v(x) =

+ 7~2 ·

2 mu

+ Y~1

Substitutegcos(~)forgv m the pre~onsexpre~onusingtheSUBST (new expre~ion,~dexpre~ion, m equ~ion) COlTlln~d.

(c4) subst(g * cos(beta), gy, ~);

Time= 60 msec. MACSYMArespondswiththeupdated expr~sionforv(z).

2

cos(beta) g rho ·

(d4)

v(x) = -

+ Vk, 2 · + Yk, l

2 mu

The boundary con~tion ~ = 0 y~|c]s kl = 0 by inspection, so incorporate t ~ s into the previous expression for ~(z) using the S U B S T command.

(c5) 8ubet(O, Y.kl, ~);

Time= 33 msec. M A C S Y M A ~spondswiththeupdated expressionforv(z).

(dS)

2 cos(beta) g rho · v(x) = Y.k2 · - 2 mu

EvMuate dv/dz as a p~Hminary step ~ ap~ymg the boundary con~tion dv/dz = 0 ~ z = 6 using the DIFF (ex- p ~ i o n or function, independent variable, order ~ de~vatlve) command.

(c7) dill(d6, x, 1);

Time= 66 msec. MACSYMA responds with the expr~sion for dv/dx.

(d7)

d

cos(beta) g rho x

-- (v(x)) = ~k2 -

dx

mu

E v ~ u a t e dv/dz a t x = 6 M & preliminary s t e p for deterndmng k2.

(C8) ev(~, · = delta, diff(v(x), x, 1) = 0);

Time= 83 msec. MACSYMAr~pondsby~vingdu/dz~ z = 6.

(d8)

cos(beta) delta g rho 0 = 7,k2 - mu

AppHcations~ computeri1~dsymboficmanlpulst~n

67

Now solvefor~ m the aboveexpres~onuslngtheSOLVE (expremlon, varia~ename) command.

(c9) solve(%, ~2);

Time= 433 msec. MACSYMArespondswiththeexpfi~texpre~fion ~rk2.

cosCbeta) delta g rho

(d9)

[~2 = .....................

]

mu

Subs~tutefork~ m boththefirstand ~condtenusof(d5).

(c10) s u b s t ( p a r t ( ~ , l , 2 ) , ~.k2, d 5 ) ;

Time= 33 msec. MACSYMA ~spondswiththefin~ e x p ~ d o n f o r t h e v d o d t y ~s~ibutionv(=).

(dl0)

cosCbeta) delta g rho x vCx) = mu

2 cosCbeta) g rho x 2 mu

Ev~uatethe v~umetrlc flow rate per unlt Q/W by integrating(d10) with respectto ~ overthefimlts(0,6) uslngthelNTEGRATE (expre~ion, in~gration variabk, ~w~limit, upp~limlt) command.

(c11) q = i n t e g r a t e ( p a r t C d l O , 2 ) , x, O, d e l t a ) ;

I s c o s C b e t a ) d e l t a g mu r h o p o s i t i v e , n e g a t i v e , o r z e r o ?

positive;

2 TS c o s ( b e t a ) d e l t a g mu rho p o s i t i v e , n e g a t i v e , o r z e r o ?

positive;

Totaltime= 10800 asec. GCtimeffi 3883 msec. MACSYMA responds with the integrated result.

(d11)

3 cos(beta) delta g rho q= 3mu

The calculations are now finished so stop using the QUIT command.

(ct2) quitC);

68

W.R. THISSELL, P.L. MILLS

APPENDIX 2 MACSYMA Solution for Ezample Access the MACSYMA program after receiving the UNIX prompt which is a ~ sign for this particular VAX computer.

7. macsyffia

This is UNIX NACSYNA Release 3 0 9 . I. (c) 1976, 1984 Massachusetts Institute of Technology. All Rights Reserved. Enhancements (c) 1984 Symbolics, Inc. All Rights Reserved. Type describe(trade_secret); to see Trade Secret notice. Type exec("man macsyma"); for help. Apply the Laplace transform to the material balance equation defined in the text by (9b) where vd = Vd, CtWO(t) = C2(t), qe = Qe, and c o n e ( t ) = C1(t) using the LAPLACE command.

(cl) laplace( vd * 'diff(ctwo(t), t) = qe * cone(t) - qe * ctwo(t), t, s);

/src/Iocal/pkg/macsyma/maxerc/laplac.o being loaded. MACSYMA responds with the Laplace transformed equation where CtWO(0) = C2(t = 0).

(dl) (s laplace(ctwo(t), t, s) - ctwo(O)) vd =

qe laplace(cons(t), t, s) - qe laplace(ctwo(t), Spedfy tl~t all execution times will be displayed after a given input command.

t, s)

(c2) showtime: all;

Time= 16 msec. MACSYMA verifies that all execution times will be evaluated.

(d2)

all

Apply the Laplace transformto the nmterialbalanceequationdefinedin the textby (9a) where va = V,, qa = Q,, and alpha = ~ using the L A P L A C E conuuand.

(c3) laplace(va*)diffCcone(t),t) =qe*ctwoCt)-qe* coneCt)-alpha*qa*coneCt), t,s);

Ti,.e= 216 msec. MACSYMA responds with the Laplace transformed equation where c o n e ( 0 ) = Cl(t = 0).

Computing the display now

(d3) (s laplace(cone(t), t. s) - cone(O)) va = qe laplace(ctwo(t), t , S) -qe laplace(cone(t), t , S) - alpha qa laplace(cone(t), t, s) Substitutethe initialconditionCa(O) = 0 in (dl) usingthe SUBST conunand.

(c9) subst(O, ctwo(O), dl);

Timer 16 msec. MACSYMA responds with the updated expression for the Laplace transformed equation.

(d9) s laplace(ctwo(t), t, s) vd = qe laplace(cone(t), t, s)

- qe laplace(ctwo(t), t, s)

Applications of computerized symbolic manipulation

69

Substitute the initial condltion CI(0) = c~Mi/~/a in (d3) using the SUBST comnumd.

(c10) eubet(a*alpha/va, cone(0), d3);

Time= 66 meec. MACSYMA responds with the updated expression for the Laplace transformed equation.

(dlO) (s l a p l a c e ( c o n e ( t ) , t , s)

alpha m ) va = Va

qe laplace(cteo(t), t, s) - qe laplace(cone(t), t, s) - alpha qa laplace(cone(t), t, s) S~ve (d9)fortheLaplacetransformofC2(t) re~rredtohere M~9(s) u~ngtheSOLVEcommand. (cll) solve(d9, laplace(cteo(t), t, s));

TJ~e= 416 msec. MACSYMAresponds withtheexpH~texpre~ion ~r~2(s).

(d11)

qe l a p l a c e ( c o n e ( t ) , t , s)

[laplace(ct.o(t), t, s) :

---]

s vd + qe

Extract the first term on the right-hand side of ( d l l ) as the initial step for combinlng ( d l 0 ) and (d11) into a single equation in ternm of the Laplace transform of Cl(t) referred to here as C'l(s). (c12) p~rt(?,l,2);

Time= 0 msec. MACSYMA responds with the desired ~rm.

(d12)

qe laplace(cone(t), t, s) s vd + qe

Substitute the above term into ( d l 0 ) to eliminate ~2(s) and obtain an hnpllcit single equation where Cl(S) is the only unknown.

(c13) subst(part(d11,l,2), part(d11,1,1), dlO);

Time= 50 msec. MACSYMA responds with the imp];cit expression for ~1 (8).

(d13) (s laplace(cone(t), t, s)

alpha va

) va =

2 qe laplace(cone(t), s vd + qe

t, s) - qe laplace(cone(t),

t, s)

- alpha qa laplace(cone(t), t, e)

70

W.~. THISSELL,P.L. MILLS

S~ve the abo~ impUcit expre~on to obt~n an explidt exp~s~on for ~1 (s) using the SOLVE conunand.

(c14) solve(~, laplace(cone(t), t, s));

Totaltime= 6633 msec. GCtime= 3750 msec. MACSYMA ~spondswiththeexpU~texp~ion ~r~l(s). (d14) [ l a p l a c e ( c o n e ( t ) , t , s) =

alpha m s vd + alpha m qe 2 (s va + (qe + alpha qa) s) vd + qe s va + alpha qa qe The Lapl~e trans~rm of the trac~ exit concentration C~t(t), re,fred to hem M Cout(s), is obtained using equation (11) in the text. (c16) laplace(cout(t), t, s)= (1-alpha)*m/qa + alpha*laplace(oone(t),t,s);

Time= 83 msec. MACSYMArespondsby ~arran~ngthe previousinputexpresslon. (1 - alpha) m (d18) laplace(cout(t), t, s)= alpha laplace(cone(t), t, s) + qa

Sutmtltutethe dgh~hand dde ~ (dl4) intothefirsttermonthe dght-hand ~deof(dl6)toobt~nanexpHdt expression ~ r ~ t ( s ) .

(c16) subst(part(d14,1,2),laplaoe(cone(t),t,s),

diS);

Time= 50 msec. MACSYMArespondswiththeexp~citexpre~ion~r Co.t(s). (dlS) laplace(cout(t), t, s) =

alpha (alpha m s vd + alpha m qe)

(1 - alpha) m +

2

qa

(s va + (qe + alpha qa) s) vd + qe s va + alpha qa qe

Determlnethelnveme Laplace trans~rrnof~out(s) ~ingtheILT ( e x p ~ i o n , Lapl~etrandormv~ia~e, inve~e Lapl~etrans~rmvariable) conuuand.

(c16) ilt(~,s,t);

2 Is 4 alpha qa qe va vd - ((qe + alpha qa) vd + qe va) positive, negative, or zero? negative; Totaltime= 8866 msec. GCtime= 3666 msec. t ((qe + alpha qa) vd + qe va) 2 va vd

Appllcations of computerized symbolic manipulation

71

Note that MACSYMA queries the user to determine whether the indicated ~ o u p o~ tcrnm k podtive, negative

or zero. Here the user has assumed the group of terms will always be ne[|at|ve. The followimg e ~

for Come(t)

is then output by MACSYMA.

(d16) cout(t) = ~e

2

2

((2 alpha m qe va vd - alpha m vd ((qe + alpha q&) vd ч qe v|))

2

22 2

s~-h(t sqrt((qe + 2 alpha qa qe + alpha qa ) vd

2

22

+ (2 qe - 2 alpha qa qe) vavd + qe va )/(2 vavd))

2

22 2

2

/sqrt((qe + 2 alpha qa qe + alpha qa ) vd +(2 qe -2 alpha qa qe) vavd

22

2

2

22 2

+ qe va ) + alpha a vd cosh(t sqrt((qe +2 alpha qa qe + alpha qa ) vd

2

22

+ (2 qe - 2 alpha qa qe) va vd + qe va )/(2 va vd)))/(va vd)

(1 - alpha) m

+ ilt(

, s, t)

qa

MACSYMA d o ~ not know h e inverse Laplace ~ansl'orm d (1 -- ~)~,/Q., .o it ~ ~ . d (1 - ~)M~6(t)/Qa using the SUBST command.

by the tmea-u

(c17) subst(delta*me(1-alpha)/qa,

ilt(m*(1-alpha)/qa,s,t),

~);

Time= 216 msec. t ((qe + alpha qa) vd + qe va) 2 va vd MACSYMA responds by incorporating the above result into the expresdoa for Ceut(t) where delta = 8(t).

(d17) cout(t) = ~e

2

2

( ( 2 a l p h a in qe v a v d - a l p h a m v d ( ( q e + a l p h a q a ) v d + q e v a ) )

2

22 2

s~nh(t sqrt((qe + 2 alpha qa qe + alpha qa ) vd

2

22

+ (2 qe - 2 alpha qa qe) va vd + qe va )/(2 va vd))

2

22 2

2

/sqrt((qe +2 alpha qa qe + alpha qa ) vd +(2 qe - 2 alpha qa qe) va vd

22

2

2

22 2

+qe va ) + alpha m vd cosh(t sqrt((qe + 2 alpha qa qe + alpha qa ) vd

2

22

+ (2 qe - 2 alpha qa qe) va vd + qe va )/(2 va vd)))/(va vd)

(1 - alpha) delta m + qa

72

W.R. THISSELL, P.L. MILLS

The inverse Lap]ace transform of ~out(a) is determined again, except the user now specifies that the indicated group of terms is positive.

(c18) ilt(dlS,s,t) ; 2 Is 4 alpha qa qe va vd - ((qe + alpha qa) vd + qe va) positive, negative,

or zero?

positive; Totaltime= 9580 meec. GCtime = 3800 msec. t ((qe + alpha qa) vd + qe va)

2 va vd

MACSYMA responds with the inverse Laplace transform, tout(t).

( d 1 8 ) c o u t ( t ) = 7~e

2

2

((2 alpha m qe va vd - alpha m vd ((qe + alpha qa)

2

22 2

sin(t eqrt(- (qe + 2 alpha qa qe + alpha qa ) vd

vd + qe

va))

2

22

- (2 qe - 2 alpha qa qe) va vd - qe va )/(2 va vd))

2

22 2

2

/sqrt(-(qe + 2 alpha qa qe + alpha qa ) vd -(2 qe -2 alpha qa qe) va vd

22

2

2

22 2

- qe va )+ alpha m vd cos(t sqrt(-(qe +2 alpha qa qe + alpha qa ) vd

2

22

- (2 qe - 2 alpha qa qe) va vd - qe va )/(2 va vd)))/(va vd)

(I - alpha) m

+ ilt(--

, s, t)

qa

The SUBST command used previously is repeated since the same group of terms appears.

(c19) c17;

TiLe= 33 msec.

(d19)

(I - alpha) delta m

(I - alpha) m

eubst(

. ilt(--

, s, t). 7.)

qa

qa

MACSYMA responds by giving the final expression for Co, t (t).

(c20) subst(delta * m *(1-alpha)/qa, ilt((1-alpha)*m/qa,s,t),dlS);

TiLe= 200 msec.

t ((qe + alpha qa) vd + qe va)

2 va vd

( d 2 0 ) c o u t ( t ) = 7.e

2

2

((2 alpha m qe va vd - alpha m vd ((qe + alpha qa) vd + qe va))

2

22 2

sin(t sqrt(- (qe + 2 alpha qa qe + alpha qa ) vd

2

22

- (2 qe - 2 alpha qa qe) va vd - qe va )/(2 va vd))

2

222

2

/eqrt(-(qe + 2 alpha qa qe + alpha qa )vd -(2 qe - 2 alpha qa qe) va vd

22

2

2

22 2

- qe va ) + alpha m vd cos(t sqrt(-(qe + 2 alpha qa qe + alpha qa ) vd

2

22

- (2 qe - 2 alpha qa qe) va vd - qe va )/(2 va vd)))/(va vd)

(1 - alpha) delta m +

qa

Appllcstions d ccmputerised symbolic manlpulat~n

73

The commands ~ven above in (C18) through ((L~O) are repeated for the case where the indicated group d terms is gem.

(c21) ilt(dlS,s,t);

2 Is 4 alpha qa qe va vd - ((qe + alpha qa) vd + qe va) positive, negative, zero; T o t a l t ~ e = 6983 -aec. GCtiae= 3700 asec.

or zero?

(d21) cout(t) =

2

2

t(2 alpha m qe va vd -alpha m vd((qe +alpha qa) vd +qe va))

2 alpha m

C 22 2 va vd

+ ........ ) va

t ((qe + alpha qa) vd + qe va)

2 va vd

(I - alpha) m

%e

+ ilt( .............

, s, t)

qa

(c22) subst(delta*m*(1-alpha)/qa,ilt((1-alpha)*m/qa, s,t),~); Time= 133 msec.

(d22) c o u t ( t ) ffi

2

2

2

t(2 alpha m qe va vd -alpha m vd((qe +alpha qa)vd +qe va)) alpha m

(

+

)

22

va

2 va vd

t ((qe + alpha qa) vd + qe va)

2 vavd ~.e

(1 - alpha) delta m ч qa

The session is then tennlnated with the QUIT conunand.

(c23) quit() ; ~. e x i t

APPENDIX 3 MACSYMA 8ol~tio# ]or Example $ Access the MACSYMA program after receiving the UNIX prompt ~ sign. ?. macsyma T h i s i s UNIX MACSYMAR e l e a s e 3 0 9 . 1 . (c) 1976, 1984 Massachusetts Institute of Technology. All Rights Reserved. Enhancements (c) 1984 Symbolics, Inc. 111 Rights Reserved. Type describe(trade_secret); to see Trade Secret notice. Type exec("man macsyma"); for help.

74

W.R. THISSELL,P.L. MILLS

Speci~ thAt allexecution times wiH be displayedaftera ~ven mputcommand. ( c l ) showtime: all;

Time= 0 msec. M A C S Y M A verifies that allexecut~ntimeswm beev~u~ed.

(dl)

all

Speci~thatthee~mentsofa 4 X 4 gener~ma~ix w~beentered, andinputthearrayv~ues for eachrowand column.

(ў2) entermatrix(4,4);

Is the matrix 1. Diagonal 2. Symmetric 3. A n t i s ~ e t r i c Answer I, 2, 3 or 4 4; Row 1 Column 1: 1100;

Row 1 Column 2: O; Row 1 Column 3: O;

Row 1 Column 4: O;

Row 2 Column 1: - 1 0 0 0 ; Row 2 Column 2: 1400; Row 2 Col-mn 3: - 1 0 0 ; Row 2 Colm.n 4: O;

Row 3 Column 1: O; Row 3 Column 2: - 1 1 0 0 ; Row 3 Column 3: 1240; Row 3 C o l u = ~ 4 : - 1 0 0 ; Row 4 coiu~m 1: o;

Row 4 Column 2: O; ROW 4 Col-mn 3: - 1 1 0 0 ;

Row 4 Column 4: 1250;

Matrix entered.

Time= 216 msec. MACSYMA responds by outputtlng the BII matrix.

(d2)

[ 1100 [ [ - 1000 C [o C Co

0 1400 - 11oo o

o - 10o 124o - 1100

o] ] 0] ] - lOO ] ] 1250 ]

4. General

Applications of computerized symbolic manipulation

75

Specify that the elements of the d matrix will be entered with 4 rows and 1 colmnn, and input the array values for each row and column.

(c3) enteraatrix(4,1) ;

Koe 1 Col--- 1: 1000;

Roe 2 Col-m- 1: O;

Row 3 Col-m- 1: O;

Roe 4 Col-m~ 1: O;

Matrix entered.

Time= 50 mite. MACSYMA responds with the input values. Cd3)

[ 1ooo ]

[

]

[o ]

[

]

[o ]

[

]

[o ]

Solve for the ,nk~own tpecles concentr&tlonl by inverting the ~ matrix defined by (d2) and multiplying the result by the d matrix defined by ( d 3 ) , i.e., U-- B d.

(c4) invert(d2) . d3;

MACSYMA responds with the exact values for the column vector U. Batching the file /src/local/pkg/macsyma/ehare/invert.mac Batching done.

Time= 2983 msec. (d4)

[ lO ] [ --

[ 11 ]

[

]

[ 28800 ] [ ..... ]

[ 41327 "]

[

]

[ 25oo ] [ .... ]

[ 3757 ]

[

]

[ 2200 ] [ .... ]

[ 3757 ]

76

W . R . THISSELL, P.L. MILLS

APPENDIX 4 MA CSYMA Sois~io. .for Ez4mple Substitution of the numerical values for Li, Vi, m, F , z0 and ~! given in the text below (28) into (26) and (28) give the following linear system of equations.

'--307 207 0 0 0 0 0 0 0 0

100 --307 207 0 0 0 0 0 0 0

X2

ro

0 100 --307 207 0 0 0 0 0 0

x3

o

0 0 100 - 3 0 7 207 0 0 0 0 0

X~

o

0 0 0 100 --307 207 0 0 0 0

x5

o

0 0 0 0 1 ~ - ~ 7 207 0 0 0

:(6

o

0 0 0 0 0 1~ -~7 ~ 0 0

X7

--0.039

0 0 0 0 0 0 1~ -1~ ~ 0

Xs

0

0 0 0 0 0 0 0 100 -190 90

X9

0

0 0 0 0 0 0 0 0 100 --190

After accessing MACSYMA, the dements of the 10 x 10 matrix corresponding to A are entered on a row-by-row basis.

(c5) entermatrix(lO,lO) ;

Is the matrix 1. Diagonal Anseer I, 2, 3 or 4 4; Row I Col-m- 1: -307; Row 1 Column 2: 207; Row I Col-m- 3: 0; Row 1 Column 4: O; Row 1 Column 5: O; Row 1 Column 6: O; Row 1 Col-m- 7: 0; Row 1 Column 8: 0; Row 1Colum- 9: O; Roe 1 Column 10: O; Roe 2 Column 1: 100; Row 2 Column 2: -307; Roe 2 Column 3: 207; Roe 2 Col-m- 4: 0; Ro. 2 Col-mn 5: 0; Row 2 Column 6: 0; Row 2 Column 7: 0; Row 2 Column 8: 0; Row 2 Column 9: 0; Row 2 Column 10: O; Row 3 Column 1: O; Eow 3 Col-mu 2: 1 0 0 ; Roe 3 Column 3: -307; Row 3 Coltuan 4: 207; Row 3 Column 5: O; Row 3 Column 6: O; Row 3 Column 7: O;

2. Sy,metric

3. Antisymnetric

4. General

Applications of computerized symbolic manipulation

77

Koe 3 Column 8: 0; Kow 3 Column 9: 0; Row 3 Column 10: 0 ; Roe 4 Col-un 1: 0; Thisp~cessisco~mued until all of the array elements have been entered as shown below.

Row 10 Col-m- 4: 0;

Row 10 Column 5: O; Row 10 Column 6: O; Row 10 Column 7: O; Rog 10 Column 8: O; Row 10 Column 9: 100; Row 10 Column 10: - 1 9 0 ; Matrix entered. Time= 1110 msec. MACSYMA then responds with the input values.

[- 307 [ [ 100 [ [0 [ [0 [ [o (d5) [ [o [ [o [ [o [ [o [ [o

207 - 307 100 0 o o o o o o

0 207 - 307 100 o o o o o o

0 0 207 - 307 1oo o o o o o

o 0 o 207 - 307 100 0 0 0 0

o

o

0

0

o

o

0

0

207

0

- 307 207

100 - 307

0

100

0

0

0

0

o

o

0

0

o

o

0

0

0

0

0

0

90

0

- 190 90

100 - 190

0

100

The el~rnents of the 10 by 1 matrix correspond to b are then entered.

(c6) entermatrix(lO,1);

Kow i Col-mn 1: - 5 ; Rog 2 Column 1: 0; Row 3 Column 1: 0; Row 4 Column 1: 0; Kow 5 C o l - m ~ 1: 0; Kow 6 Column 1: 0; Row 7 C o l - m - 1: - 0 . 0 3 9 ; Row 8 Column 1: 0; Rog 9 Col-m- 1: O; Row 10 C o l - m - 1: 0; Matrix entered. Time= 133 msec.

o] ] 0] ] o] ] 0] ] 0] ] 0] ] 0] ] 0] ] 90 ] ] - 190]

78

W.R,. THISSBLL, P.L. MILLS

MACSYMA responds with the input values. Cd6)

[ -5 ]

[

]

[o ]

[

]

[o ]

[

]

[o ]

[

]

[o ]

[

]

[o ]

[

]

[ - 0.039 ]

[

]

[o]

[

]

[o]

[

]

[o]

Next, solve for the unknown concentrations of aniline in the aqueous phase by inverting the A matrix defined by ( d 5 ) and multiplyins the result by the b matrix defined by ( d 6 ) , i.e., - = ~-I~.

(c7) invertCdS) . d6;

Totaltime= 1006516 msec. (dT)

GCtime= 167666 msec.

[ 0.02426897197245182 ]

[

]

[ 0.01183852364996478 ]

[

]

[ 0.005833476151178772 ]

[

]

[ 0.002932487021330458 ]

[

]

[ 0.001531042997249146 ]

[

]

[ 0.0008540168986591411 ]

[

]

[ 0.0005269511505480273 ]

[

]

[ 0.0004152479261370032 ]

[

]

[ 0.0002911332323469764 ]

[

]

[ 0.0001532280170247244 ]

The session is then tennlnated with the QUIT command.

Cc8) quitC) ;

WR Thissell, PL Mills

doc.uments.com

About Us :: Privacy Policies :: Terms of Service :: Feedback :: Copyright :: Contact Us :: DMCA Policy

Copyright © 2018 doc.uments.com