Introduction

The future of general purpose scientific and engineering computing belongs to software systems that provide exceptional integration (pun intended) of several key functions: symbolic mathematics, graphics, numerical computations, document creation, and an ability to link to other software tools (e.g. spreadsheets, programming languages, numerical libraries, databanks, etc.). While we do not yet have tools that perform all of these tasks with the required level of excellence, the class of software that at present comes closest to this ideal are computer algebra systems (CASs). There are several computer algebra systems in use today: Macsyma, Reduce, Derive, Mathematica, Maple and Scratchpad are some of the better known ones. Many modern CASs integrates (if we may be permitted to use the word) symbolic mathematical capabilities (integration and differentiation, for example) with numerical capabilities (e.g. integration of ODE systems and linear and nonlinear equation solving) and sophisticated graphics (including 3-dimensional plots of (parametric) surfaces, spacecurves, and much more). Other tools are useful for the numerical and graphical calculations but only a CAS can help us with the symbolic computation.

For the past few years the author has made extensive use of Maple to tackle a wide variety of problems in chemical engineering (primarily thermodynamics) [1-10], numerical computing [10,11] and cartography [12]. This article provides an overview of some of our work with Maple (for complete details readers should consult the references).

Thermodynamics with Maple

Thermodynamics deals with relationships between eight state variables normally assigned the letters P, v, T, U, H, S, A, and G (these letters stand for Pressure, Volume, Temperature, Internal Energy, Enthalpy, Entropy, Helmholtz Energy, and Gibbs Energy respectively). In order to be useful in the exploration of thermodynamics, a computer algebra system needs:

1. The ability to differentiate undefined functions of several variables including undefined functions of an undefined number of variables.

2. To be able to work with the subscripted (indexed) partial derivatives of thermodynamic functions found in many thermodynamic formulae.

3. The ability to differentiate arbitrary sums (a summation where the upper limit is a symbol rather than a number) with respect to an indexed variable.

4. The ability to display the specialized graphical images required in thermodynamics.

To the best of our knowledge none of the major commercial computer algebra systems can do all (or, in some cases, any) of the above right out of the box. However, most CASs can be taught new tricks and we have done this in a package for Maple called TDtools that can "do" all of the above, as illustrated below.

Derivatives of Undefined Functions

Consider an arbitrary thermodynamic property, X, and assume it to be a function of two other variables, Y and Z say.

> xdef:=X=X(Y,Z): xdef;

In order to explore the world of thermodynamic property relations we need to be able to obtain differentials of undefined functions like this one. Maple's built in differential operator, D, is unable to perform this operation. We have, therefore, created a new operator called TD that is able to take the total differential of such functions as shown below. In all other respects, TD is supposed to behave like D.

> alias(d=D):

> t0:=TD(xdef): t0;

Partial Derivatives in Thermodynamics

Note the indexed partial derivatives in the above expression. Such derivatives are important in thermodynamics and are returned by the TDiff procedure.

> TDiff(X,Y);

If TDiff is called with more than two arguments, the additional arguments are assumed to be constant variables and become (or are added to) the index.

> TDiff(X,Y,Z);

Second (and higher) derivatives can be obtained by repeated application of TDiff.

The functions TD and TDiff (as well as a few others) allow us easily to derive many other relations between partial derivatives: the triple product rule,

the chain rule,

the insertion rule,

and the following expression (see [8] for details).

This seemingly complicated result is of immense utility for deriving relationships between thermodynamic properties. We may encapsulate Maple results in procedures for later use.

> DXYZ := proc(X,Y,Z,K,L)

> Diff(X,Y)[Z] = (Diff(X,L)[K] * > Diff(Z,K)[L] - Diff(X,K)[L] *
> Diff(Z,L)[K]) /

> (Diff(Y,L)[K] * Diff(Z,K)[L]-
> Diff(Y,K)[L] * Diff(Z,L)[K]);

> end:

Some of the partial derivatives in the above expression can be either 0 or 1 and these simplifications are encoded into `simplify/TDiff` to which we include a call inside the DXYZ procedure shown below.

> DXYZ := proc(X,Y,Z)

> local expr, L, K;

> expr:=Diff(X,Y)[Z] = (Diff(X,T)[P] > * Diff(Z,P)[T] - Diff(X,P)[T] * > Diff(Z,T)[P]) /

> (Diff(Y,T)[P] * Diff(Z,P)[T]-
> Diff(Y,P)[T] * Diff(Z,T)[P]);

> expr := simplify(expr,TDiff);

> expr := simplify(Right(
> simplify,expr,Thermo));

> end:

We have cut down the number of arguments to three since the last two are T and P. We have also added a call to `simplify/Thermo` that knows the Maxwell relations, the definition of heat capacity and a few other thermodynamic property relations. The code is not shown here since it is lengthy; however, we should note that all of the simplifications contained within that procedure can be derived using Maple. We can use this procedure to express any thermodynamic partial derivative in terms of measurable properties (PvT data, heat capacities) and entropy. Some examples follow.

> DXYZ(G,A,U);

> DXYZ(S,v,A);

All of the more than 300 different partial derivatives can now be presented error-free and completely automatically in just a few seconds. Alternative forms (perhaps involving compressibility) can be obtained by adding to the set of relations known to `simplify/Thermo`.

Differentiation of Arbitrary Sums

The NRTL model for the (dimensionless) Gibb's free energy of a mixture is:

> Qdef:=Q=sum(x[i] * (sum(x[j] *
> tau[j,i]*G[j,i],j=1..c) /sum(x[j] * > G[j,i],j=1..c)),i=1..c): Qdef;

This expression is the starting point for working out expressions for the activity coefficients of each species in the mixture defined by:

> ln(gamma[i])=Diff(n[t]*Q,n[i]);

The differentiation required by the two prior expressions involves differentiation of an undefined or arbitrary sum, where the upper limit is a symbol that represents some presently unknown number of species. Maple's own differentiation functions do not know how to differentiate such objects with respect to one of the indexed variables; TDiff does.

The next step in the derivation is to replace the mole fractions by their mole number ratios (and to multiply the result by the total number of moles):

> Q1:=n[t]*subs(seq(x[v]=n[v]/n[t],
v=[i,j,k,l,m,n]),Qdef):

> Q2:=expand(%): Q2;

and invoke TDiff to obtain the derivative of Q with respect to the number of moles of the i-th species (the parameters are assumed constant).

> t1:=ln(gamma[i])=
TDiff(rhs(Q2),n[i],[tau,G]):

The result of the above operation is quite correct but is not shown for reasons of space. In what follows we use Maple to present the result in a somewhat more familiar form using mole fractions rather than mole numbers.

> subs(seq(n[v]=x[v]*n[t],

> v=[i,j,k,l,m,n]),t1):

> expand(expand(%));

All this can, of course, be automated and the same commands can be used for any activity coefficient model [6]. Similar procedures may also be developed for deriving expressions for fugacity coefficients in multicomponent mixtures [6,7].

Maple Graphics

Maple possesses an extensive set functions for creating a wide variety of graphical effects (the plots and plottools packages). There remain several plot types that are of interest in thermodynamics but which cannot be displayed by standard Maple (curves, surfaces, and vector fields in triangular or tetrahedral coordinates, for example). Fortunately, it is reasonably simple to create new plot types in Maple. In the accompanying figures we show some of the thermodynamic images we have created - often surprisingly easily - with Maple. Limitations of space prevent us from showing in detail exactly how these plots are created, interested readers should consult the cited references and our web site. This author has found the Programming Guide that is supplied with Maple and the book by Klimek and Klimek [13] to be particularly helpful in explaining how to dismantle and reassemble Maple plots.

The fundamental relationships between state variables were explored in three famous papers by J.W. Gibbs [14]. In the second of these papers Gibbs explores the geometry of the U(S,v) surface. It is interesting to note that there are very few illustrations in Gibbs' papers; it seems that Gibbs could "see" the relationships between state variables. There should be little doubt that visualization of thermodynamic relations can aid an understanding of what for many is a difficult subject. Not all of us are able to see thermodynamics in our heads as easily as did Gibbs, but few texts contain graphical visualizations of any of what we refer to here as Gibbs surfaces. An exception is the advanced text of Tester and Modell [15] which contains a few images created using computer software developed by Jolls and Coy [16-19]. Their 12,000 line Fortran program, initially created to run on a few specific platforms and later ported to others, provides some wonderful graphic images of fundamental thermodynamic functions.

Figures 1 to 3 show Gibbs’ surfaces created using Maple. Given the parametric equations (in terms of temperature and volume) Maple can create these surfaces with a one-line command. The black and red lines are the binodal and spinodal curves, respectively. These lines require some additional computational work. Several other examples and complete details of how to construct these figures are given by Baur et al. [9]. It is very instructive to view these surfaces from different vantage points, a task that is made considerably easier in the latest version of Maple which can rotate three dimensional images in real time.

Figure 1: Helmholtz energy surface as a function of (dimensionless) temperature and volume.

Additional examples, including some frames from an animation showing the tangent plane sliding over the U(S,v) surface are given by Baur et al. (1998).

Figures 4 and 5 show residue curves and a vector field diagram for a mixture of acetone, methanol, and chloroform. These figures require extensive numerical computation prior to their construction; details are available in [8]. Maple can create vector field plots out of the box, but not for the triangular field needed here. Creating ones own field plots is very simple given the tools available in Maple.

Figure 2: Internal energy surface as a function of (dimensionless) entropy and volume.

Figure 3: Gibbs energy as a function of temperature and pressure.

Figure 4: Residue curves for acetone, methanol, and chloroform.

Figure 5: Vector field diagram for acetone, methanol, and chloroform.

The location of the binary and ternary azeotropes is clearly visible in Figures 4 and 5. The computation of all azeotropes in a multicomponent mixture has been the subject of a number of recent papers and sophisticated algorithms for this purpose have been published (e.g. [20]). We have found it possible to compute all azeotropes using nothing more than a constrained Newton’s method (at least, we find all the azeotropes listed in [20], the difference is that we cannot prove that we have them all). The result of our computations (see [8] for the code) is a Maple table, an example for this system is:

Note that our code has determined the composition and temperature of the azeotropes as well as whether it is a stable node (SN), unstable node (UN) or saddle point (S). We can pass this table to another procedure that works out the distillation boundaries and draws an approximate distillation boundary map as shown below. The procedure DBMplot is an implementation in Maple of the algorithm of Foucher et al. [21] with the exception that we have the advantage of knowing the stability of the azeotropic points.

> DBMplot(Azeotable3);

Figure 6: Approximate distillation boundary map for acetone, methanol, chloroform.

Numerical Computing

The emphasis on symbolic manipulation in any definition of computer algebra should not be taken to imply that a CAS is unsuitable for the predominantly numerical calculations that are such an important part of modern engineering computing. It is true that numerical computing is not Maple's strongest suit, but the program does possess a number of functions designed specifically for such activities. The built in numerical equation solver, fsolve, is an excellent tool for solving single nonlinear equations; sadly it is much less successful on problems involving systems of nonlinear equations. In view of this weakness we have developed a fairly sophisticated implementation of Newton's method for Maple that allows the user considerable control over the computations. Early versions of Newton were supplied with Maple in the share library of user supplied functions. The latest and considerably improved version is available from the authors web site. We have also implemented homotopy continuation methods in Maple, but that code has not yet been released as part of the Newton package (although we do intend to do so in the near future).

Many engineering problems require the numerical solution of differential equations, some require the solution of mixed systems of differential and algebraic equations (DAEs). Maple comes with a number of built in numerical intergration methods including Runge-Kutta, Gear's method and the LSODE algorithms. It is not, however, capable of solving DAEs. To deal with this shortcoming we have developed BESIRK [10], an implementation in Maple of Michelsen's 3rd order SIRK method combined with a Bulirsch-Stoer extrapolation technique. The BESIRK code can handle stiff systems of ODEs significantly faster than Maple's own integration routines and was used in the computation of the residue curves shown in Figure 4. We have also found the BESIRK code to be of considerable value in solving PDEs by the Method of Lines[11].

Numerical problem solving using interval methods has been used to solve a number of chemical engineering problems that often are considered “difficult” [22]. Even out of the box, Maple is capable of some elementary interval computations. A more complete package for interval computing with Maple has been developed by Rob Corless and Amanda Connell of the University of Western Ontario and supplied with Maple as part of its share library. George Corliss from Marquette University and his students developed an interval Newton method that uses the Corless/Connell package. However, their code does not work in the more recent releases of Maple (3, 4 and 5). We have modified the Corliss code for Release 4; the code is available on our web site.

Those of us that still prefer to program in more traditional languages may like to know that Maple can also be used to turn expressions and procedures into Fortran or C code. Maple has its own built in functions for this purpose and there is a package called Transfor developed by Claude Gomez that can be used for this purpose. Gomez is also a developer of Scilab, a (free) scientific software package for numerical computing in a user friendly environment. Scilab (with which this author has no personal experience) has links to high powered numerical algorithms (e.g. DASSL) and can carry out symbolic computing via a link to Maple.

More Maple

Readers wishing to experiment for themselves with Maple may find sample worksheets to be of some value. Complete worksheets for all of the above (and many more) exmaples are available on the authors web site. Worksheets by the author that demonstrate the use of Maple for solving simple material balance problems are supplied with Maple as part of the share library (in /share/engineer/chemeng). The author has also contributed 10 Maple worksheets to the CD that accompanies a new book by Cutlip and Schacham [23] (and many more are under development).

Perhaps the single most useful additional capability of Release 4 of Maple V (the current version is Release 5) was its ability to read binary random access files of arbitrary format. Maple's new found abilities at handling binary data were, however, flawed since, other than conversion to and from text strings, Release 4 did not come with utilities to convert binary data into the integers and floating point numbers that could be used in other computations. Following a request by this author posted on the internet, electrical engineer Joe Riel provided the functions needed to make it very easy to manipulate binary data within Maple. We have, for example, used his package to develop code for reading the databank that comes with ChemSep, thus providing Maple with access to physical property data for use in many engineering computations. Riel has also written a units conversion package that this author uses in a great many worksheets. Riel's packages (ieee and mks) are available on the web site of Alex Walz. Walz is also the author of some of the best mathematical and utility packages that significantly enhance the ease of use of Maple.

Acknowledgements

Several students, both graduate and undergraduate, contributed to the development of the Maple codes used in this work. They are: Arnoud Higler, Richard Baur, Dorothy Scheckler, James Reese, and Warren Hoffmaster.

References

1. Taylor, R. and K. Atherley, Chemical Engineering with Maple, Chem. Eng. Ed., 56-61, Winter 1995.

2. Taylor, R., Thermodynamics with Maple. I – Equations of State, Maple Tech, Issue 10, 50-59 (1993).

3. Taylor, R., Thermodynamics with Maple. II – Phase Equilibria in Binary Systems, Maple Tech, 1(1), 83-92 (1994)

4. Adams, S. and R. Taylor, Thermodynamics with Maple. III - Thermodynamic Property Relations and the Maxwell equations, Maple Tech, 1(2), 68-81, 1994.

5. Taylor, R., Thermodynamics with Maple. IV – The Properties of Steam, Maple Tech, 3(2), 61-68, 1996.

6. Taylor, R., Automatic Derivation of Thermodynamic Property Functions using Computer Algebra, Fluid Phase Equilibria, 129, 37-47, 1997.

7. Taylor, R., Thermodynamics with Maple. I - Symbolic Computation, Mathematics and Computation in Simulation, 45, 101-119, 1998.

8. Taylor, R., Thermodynamics with Maple. II - Numerical and Graphical Applications, Mathematics and Computation in Simulation, 45, 121-146, 1998.

9. Baur, R. J. Bailey, B. Brol, A. Tatusko, and R. Taylor, Maple and the Art of Thermodynamics, Computer Applications in Engineering Education, 6, 223-234, 1998.

10. Schwalbe, D., H.A. Kooijman, and R. Taylor, Solving Stiff Differential Equations and Differential Algebraic Systems with Maple V Maple Tech, 3(2), 47-53, 1996.

11. Kooijman, H.A. and R. Taylor, Numerical Computing with Maple - Solving PDEs by the Method of Lines, under review, 1998.

12. Taylor, R., R. Baur, and J. Oprea, Maple Maps, under review, 1998.

13. Klimek, G. and M. Klimek, Discovering Curves and Surfaces with Maple, Springer, New York, 1997.

14. Gibbs, J.W., The Scientific Papers of J.W. Gibbs, Ph.D., LLD., Volume 1 - Thermodynamics, Dover, New York, 1961.

15. Testor, J.W. and M. Modell, Thermodynamics and its Applications, 3rd edition, Prentice-Hall, Inc, 1996.

16. Jolls, K.R. Gibbs and the Art of Thermodynamics, Proc. Gibbs Symposium (G.D. Mostow and D.G. Caldi, eds.), 293, American Mathematical Society, 1990.

17. Coy, D.C., Visualizing Thermodynamic Stability and Phase Equilibrium through Computer Graphics. Ph.D. thesis, Iowa State University, 1993.

18. Jolls, K.R. Understanding Thermodynamics through Interactive Computer Graphics, Chem. Eng. Progress, 64-69, Feb. 1989.

19. Jolls, K.R., M.C. Schmitz, and D.C. Coy, Seeing is Believing: A New Look at an Old Subject. The Chemical Engineer, 42-46, May 30, 1991.

20. Fidkowski, Z.T., M.F.Malone, and M.F. Doherty, Computing Azeotropes in Multicomponent Mixtures, Comput. Chem. Eng., 17, 1141-1155, 1993.

21. Foucher, E.R., M.F. Doherty, and M.F. Malone, Automatic Screening of Entrainers for Honogeneous Azeotropic Distillation, Ind. Eng. Chem. Research, 29, 760-772, 1991;

22. Hua, J.Z., J.F. Brenneke, and M.A. Stadtherr, Ind. Eng. Chem. Research, 37, 1519-1527, 1998; Fluid Phase Equilibria, 116, 52-59, 1996.

23. Cutlip, M.B., and M. Schacham, Problem Solving in Chemical Engineering with Numerical Methods, Prentice-Hall, 1998.

Web Sites

Claude Gomez:

www-rocq.inria.fr/scilab/gomez/transfor/transfor.html

R. Taylor: www.clarkson.edu/~chengweb/faculty/taylor/taylor.html

Alexander F. Walz: