Front Cover
 Editorial: Engineers Deserve a...
 Author Guidelines
 Table of Contents
 Jennifer Sinclair Curtis of the...
 Chemical Engineering at Tufts...
 Chemical Kinetics, Heat Transfer,...
 Active Problem Solving and Applied...
 Student Ratings of Teaching: Myths,...
 Introducing Stochastic Simulation...
 Applications of the Peng-Robinson...
 Celebrating ChE in Song
 Can I Trust This Software Package?...
 Back Cover

Chemical engineering education
http://cee.che.ufl.edu/ ( Journal Site )
Full Citation
Permanent Link: http://ufdc.ufl.edu/AA00000383/00174
 Material Information
Title: Chemical engineering education
Alternate Title: CEE
Abbreviated Title: Chem. eng. educ.
Physical Description: v. : ill. ; 22-28 cm.
Language: English
Creator: American Society for Engineering Education -- Chemical Engineering Division
Publisher: Chemical Engineering Division, American Society for Engineering Education
Publication Date: Winter 2008
Frequency: quarterly[1962-]
annual[ former 1960-1961]
Subjects / Keywords: Chemical engineering -- Study and teaching -- Periodicals   ( lcsh )
Genre: serial   ( sobekcm )
periodical   ( marcgt )
Citation/Reference: Chemical abstracts
Additional Physical Form: Also issued online.
Dates or Sequential Designation: 1960-June 1964 ; v. 1, no. 1 (Oct. 1965)-
Numbering Peculiarities: Publication suspended briefly: issue designated v. 1, no. 4 (June 1966) published Nov. 1967.
General Note: Title from cover.
General Note: Place of publication varies: Rochester, N.Y., 1965-1967; Gainesville, Fla., 1968-
 Record Information
Source Institution: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: oclc - 01151209
lccn - 70013732
issn - 0009-2479
sobekcm - AA00000383_00174
Classification: lcc - TP165 .C18
ddc - 660/.2/071
System ID: AA00000383:00174

Table of Contents
    Front Cover
        Page i
    Editorial: Engineers Deserve a Liberal Education, C. Judson King
        Page ii
    Author Guidelines
        Page iii
    Table of Contents
        Page 1
    Jennifer Sinclair Curtis of the University of Florida
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
    Chemical Engineering at Tufts University
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
    Chemical Kinetics, Heat Transfer, and Sensor Dynamics Revisited in a Simple Experiment
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
    Active Problem Solving and Applied Research Methods in a Graduate Course on Numerical Methods
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
    Student Ratings of Teaching: Myths, Facts, and Good Practices
        Page 33
        Page 34
    Introducing Stochastic Simulation of Chemical Reactions Using the Gillespie Algorithm and MATLAB: Revisited and Augmented
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
    Applications of the Peng-Robinson Equation of State Using Mathematica
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
    Celebrating ChE in Song
        Page 52
    Can I Trust This Software Package? An Exercise in Validation of Computational Results
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
    Back Cover
        Page 60
Full Text

chemical engineering education



on the



... of the University of Florida

Introducing Stochastic Simulation of Chemical Reactions Using the Gillespie Algorithm
and MATLAB: Revisited and Augmented (p. 35)
Argoti, Fan, Cruz, Chou
Applications of the Peng-Robinson Equation of State Using Mathematica (p. 47)
Random Thoughts: Student Ratings of Teaching: Myths, Facts, and Good Practices (p. 33)
Felder, Brent
Editorial: Engineers Deserve A Liberal Education (inside front cover)
Celebrating ChE In Song (p. 52)
Active Problem Solving and Applied Research Methods in a Graduate Course on Numerical Methods (p. 231
Maase, High
Chemical Kinetics. HeatTransfer, and Sensor Dynamics Revisited in a Simple Experiment (p. 17 1
Sad, Sad, Castro, Garetto
Can I TrustThis Software Package? An Exercise in Validation of Computational Results (p. 53)
Shacham. Brauner, Ashurst. Cutlip


Tufts University



Engineers Deserve a

A s much as or more than for any other career, en-
gineers need the benefits of a liberal bachelor's
education. Other professions have recognized
this need for many years. It is time for engineering
to change to the master's as the professional degree.
Bachelor's education should then be modeled after pre-
medical education-a set of generally required courses
within thecontextofa liberal arts bachelor's degree, with
a wide choice of majors. I have presented this proposal
and the rationale for it in more detail elsewhere.
Some of the drivers for this change are
Engineers need the flexibility to move among
careers-to management, public service, or
out of the engineering profession altogether.
SEngineers must interact effectively with people
of different backgrounds and with the public on
issues such as environment, energy, food sup-
ply, water supply, and national security.
SWith instant broadband communication and
access to information, these days "the world
is flat." Straightforward engineering jobs are
going overseas, leaving U.S. engineers with
more complex, interactive jobs and the need
to understand and work closely with other
* There is simply more to be had from life if
one has a broader background and interests.
All these factors call for a broader engineer, with
more understanding of society and the human nature
and condition.
Engineering is the only major profession with the
bachelor's degree as the first professional degree and
accreditation primarily at that level. Other professions,
such as medicine, law, business, architecture, planning,
journalism, and public health, build upon a liberal arts
bachelor's degree.
To its credit, ABET has recently sought more breadth
is that a nominally four-year bachelors degree is just too
small a boxin which to package both a professional engi-
neering education and an adequate liberal education.
Some will say that fewer students would enter engi-
neering becauseofthe prospectsoflongertimeand more

Liberal Education

C. Judson King is Provost and
Senior Vice President, Emeritus of
the University of California. He is
also Professor Emeritus of Chemi-
cal Engineering on the Berkeley
campus, where he was formerly \
Provost - Professional Schools and
Colleges and Dean of the College of
Chemistry. He presently directs the
Center for Studies in Higher Educa-
tion at Berkeley.

cost. But the change will also enable students to make a
much more informed choice as the end of the bachelor's
degreeapproaches. Breadth and theopportunityfora de-
layedcommitmentshould bring morestudentstoconsider
engineering seriously, especially women and minorities.
The community college transfer route will be much more
attractive and workable. The added time and cost are
excellent investments for a future career, and the time is
offset bythefactthat mostengineering students nowtake
well overfouryears to completetheir degree.The change
to model pre-medicaleducation mayactuallyincreasethe
number and diversity of engineering students.
Some pre-engineering bachelor's graduates will not
proceed onward to the professional master's.They will
still have good options, including other professional de-
grees and jobs linking engineering with otherfunctions.
They will contribute toward the goal of a technologically
literate population.
Some will observe that there is a strong industrial de-
mand for bachelor's engineers. But here the interests of
corporations and individuals diverge. Corporations will
gladly hire engineers for entry-level functions when the
openingsarethere. Individual engineers need to enhance
theirfiexibility and opportunities for moving to different
functions astheircareerdevelops and/orastheeconomic
cycle turns downward.
Still others will observe that it will be difficult to make
this change. But medicine and law did it in the early parts
ofthelastcentury.Surelyengineering-the problem-solv-
ing profession-can find the way to make the change.
My class at Yale had its 50th reunion this year. In the
reflective 50-yearclass book, oneofour numberobserved
and regretted that he had"missed hisYaleeducation"be-
causeof having majored in engineering.With thechange
that I am urging, he would have had that education. p

Author Guidelines for the



The laboratory experience in chemical engineering education has long been an integral part
of our curricula. CEE encourages the submission of manuscripts describing innovations in the
laboratory ranging from large-scale unitoperations experimentsto demonstrationsappropriate
for the classroom. The following guidelines are offered to assist authors in the preparation of
manuscripts that are informative to our readership. These are only suggestions, based on the
comments of previous reviewers; authors should use their own judgment in presenting their
experiences. A set of general guidelines and advice to the author can be found at ourWeb site:

c Manuscripts should describe the results of original and laboratory-tested ideas.
The ideas should be broadly applicable and described in sufficient detail to
allow and motivate others to adapt the ideas to their own curricula. It is noted
that the readership of CEE is largely faculty and instructors. Manuscripts must
contain an abstract and often include an Introduction, Laboratory Description,
Data Analysis, Summary of Experiences, Conclusions, and References.
* An Introduction should establish the context of the laboratory experi-
ence (e.g., relation to curriculum, review of literature), state the learning
objectives, and describe the rationale and approach.
* The Laboratory Description section should describe the experiment in
sufficient detail to allow the reader to judge the scope of effort required
to implement a similar experiment on his or her campus. Schematic dia-
grams or photos, cost information, and references to previous publica-
tions and Web sites, etc., are usually of benefit. Issues related to safety
should be addressed as well as any special operating procedures.
SIf appropriate, a Data Analysis section should be included that concisely
describes the method of data analysis. Recognizing that the audience
is primarily faculty, the description of the underlying theory should be
referenced or brief. The purpose of this section is to communicate to the
reader specific student-learning opportunities (e.g., treatment of reac-
tion-rate data in a temperature range that includes two mechanisms).
*The purpose of the Summary of Experiences section is to convey the
results of laboratory or classroom testing. The section can enumerate,
for example, best practices, pitfalls, student survey results, or anecdotal
* A concise statement of the Conclusions (as opposed to a summary) of
your experiences should be the last section of the paper prior to listing

Chemical Engineering Education
Department of Chemical Engineering
University of Florida * Gainesville, FL 32611
PHONE and FAX: 352-392-0861
e-mail: cee@che.ufl.edu

Tim Anderson

Phillip C. Wankat

Lynn Heasley

James O. Wilkes, U. Michigan

William J. Koros, Georgia Institute of Technology


John P. O'Connell
University of Virginia

C. Stewart Slater
Rowan University

University of Colorado
Jennifer Curtis
University of Florida
Rob Davis
University of Colorado
Pablo Debenedetti
Princeton University
Diane Dorland
Thomas F. Edgar
University of Texas at Austin
Richard M. Felder
North Carohna State University
H. Scott Fogler
University of Michgan
Carol K. Hall
North Carolina State University
Jim Henry
University of Tennessee, Chattanooga
Steve LeBlanc
University of Toledo
Ron Miller
Colorado School of Mines
Susan Montgomery
University of Michgan
Ronald W. Rousseau
Georgia Institute of Technology
Stan Sandler
University of Delaware
Donald R. Woods
McMaster University

Chemical Engineering Education
Volume 42 Number 1 Winter 2008

10 Chemical Engineering at Tufts University
Nakho Sung, and Daniel Ryder

2 Jennifer Sinclair Curtis of the University of Florida
Nicholas Peppas, Gintaras Reklaitis, and Erik Ydstie

33 Student Ratings of Teaching: Myths, Facts, and Good Practices
Richard M. Felder and Rebecca Brent

23 Active Problem Solving and Applied Research Methods in a Graduate
Course on Numerical Methods
Eric L. Maase and Karen A. High

35 Introducing Stochastic Simulation of Chemical Reactions Using the
Gillespie Algorithm and MATLAB: Revisited and Augmented
A. Argoti, L.T Fan, J. Cruz, and S.T Chou

17 Chemical Kinetics, Heat Transfer, and Sensor Dynamics Revisited in a
Simple Experiment
Maria. E. Sad, Mario R. Sad, Alberto A. Castro, and Teresita F Garetto

52 Celebrating ChE In Song
Alan M. Lane

47 Applications of the Peng-Robinson Equation of State Using Mathematica
Housam Binous

53 Can I Trust This Software Package? An Exercise in Validation Of Com-
putational Results
Mordechai Shacham, Neima Brauner, W Robert Ashurst, and
Michael B. Cutlip

inside front cover Editorial: Engineers Deserve A Liberal Education
C. Judson King

CHEMICAL ENGINEERING EDUCATION (ISSN 0009-2479) is published quarterly by the Chemical Engineering
Division, American Societyfor Engineering Education, and is edited at the University ofFlorida. Correspondence regarding
editorial matter, circulation, and changes of address should be sent to CEE, Chemical Engineering Department, University
of Florida, Gainesville, FL 32611-6005. Copyright @ 2005 by the Chemical Engineering Division, American Societyfor
Engineering Education. The statements and opinions expressed in this periodical are those of the writers and not necessarily
those of the ChE Division,ASEE, which body assumes no responsibility for them. Defective copies replaced if notified within
120 days of publication. Writefor information on subscription costs and for back copy costs and availability. POSTMASTER:
Send address changes to Chemical Engineering Education, Chemical Engineering Department., University of Florida,
Gainesville, FL 32611-6005. Periodicals Postage Paid at Gainesville, Florida, and additionalpost offices (USPS 101900).

Vol. 42, No. 1, Winter 2008

rij] educator
---- U s_____________________________________

Jennifer Sinclair Curtis

of the University of Florida

The University of Texas
Purdue University
Carnegie Mellon University
P professor Jennifer Sinclair Curtis
of the University of Florida cel-
ebrates 20 years of an academic
career that has made her not only an
internationally known researcher and
educator, but also a great role model
for younger generations of chemical
engineers. During her life, she has
faced family adversities and almost
insurmountable crises. But through
her perseverance and convictions,
as well as a deep Christian faith, she
has triumphed. She has touched many
others with her concern, enthusiasm,
and dedication. She is now recognized as a leading figure in
the chemical engineering profession.

Jennifer is a wonderful educator recognized with numerous
teaching awards, a student counselor par excellence, and a
leading administrator (with two chairmanships -one in Fresh-
man Engineering at Purdue University and another now in
Chemical Engineering at the University of Florida-as well
as service as an associate dean of Engineering at Purdue). She
is a national leader in numerous societies, including AIChE
where she has been a major motivating force of the Particle
Technology Forum, and as a recently elected AIChE director.
But Jennifer Curtis is also the internationally known leader
in particle technology, an innovative engineer, and a scholar.
Along with all these achievements she enjoys a wonderful fam-

l ily life with her husband, Barry Curtis,
and her children, Jennette and Derek.

Jennifer is a product of the rich Mid-
western tradition of hard work and ap-
preciation for education. She was born
Jennifer Lynn Doloresco, the daughter
of Jerry Doloresco and Carolyn Stortz
Doloresco, on Christmas Day 1960 in
Cincinnati, OH. Her father enlisted
in the Navy immediately after high
school to acquire enough money to
go to Miami University of Ohio. He
. Received a B.S. in accounting in three
years -a necessity because he didn't
have enough money to go for a lon-
ger period of time. Nicholas Peppas
remembers, "Jerry Doloresco was
especially happy when Jennifer went
on to graduate school and received her
Ph.D. from Princeton University. He was sad that he could
not afford more education himself but he was so happy his
children (Jennifer and her brother, Bryan Doloresco) could. I
guess this is a pretty common story from the generation before
ours that grew up in the post-Depression period. Jerry was a
loving father, always proud of Jennifer's achievements. We
kept frequent correspondence, and he was particularly proud
when Jennifer started her academic career at Carnegie Mellon
University." Jerry passed away from cancer in 1995, when
Jennifer was 35. Jennifer's mother, Carolyn, is of German
descent and still lives in Cincinnati. This marvelous mixture
of Italian and German background has given Jennifer her
wonderful enthusiastic spirit and optimism that everything
will go well, but also her dedication-bordering on stubborn-
ness-to achieve her goals

� Copyright ChE Division of ASEE 2008

Chemical Engineering Education

Jennifer graduated from Northwest High School in Cincin-
nati in May 1979 as valedictorian of her class. Having worked
in a summer job at Procter & Gamble, she was sure she wanted
to become a chemical engineer and enrolled in the School of
Chemical Engineering at Purdue in the fall of 1979. While
taking Honors Chemistry during her freshman year she met
fellow student Gavin Sinclair. They started dating in their
sophomore year and were married in 1981.
At that time, the School of Chemical Engineering at Pur-
due had instituted a new advising program for all students.
Each professor was in charge of 30 undergraduates with full
responsibility for advising and registration.
Nicholas Peppas remembers, "I was a young associate
professor, just in my third year at Purdue, and I had been as-
signed my first advising group of 30 incoming sophomores.
My first advisee was Jennifer. And to the door of my office
at CMET 210A came this shy 19-year-old with golden locks,
holding an already-filled-out departmental class registration
form. Without even stepping into the room she extended her
right hand and announced 'Here is my form, please sign it.' I
asked her to come in with the explanation that I'd like to learn
more about her family and educational background, as well as
her future goals. We ended up talking for an hour. I was truly
impressed with her intellect and enthusiasm. My intrigue was
similar to that of, say, a collector who finds a rare art object.
There and then we decided that Jennifer would start doing
research with me during her sophomore year."
Jennifer's performance at Purdue was stellar. She was se-
lected a top sophomore woman in engineering in 1981, a top
junior woman in engineering in 1982, and a top senior woman
in engineering in 1983. She received an Olin Fellowship for
undergraduate research, the first prize of the 1983 University
of Akron Polymer Engineering competition, and a best paper
award at the Mid-Western regional AIChE competition at
Dayton, OH, in 1983. In addition, she received the Best ChE
Senior awardin 1983. She graduated in May 1983 with "high-
est distinction," the Purdue equivalent of a summa cum laude.
Next to her was another graduate-her husband of two years
Gavin Sinclair, who had achieved the impossible, i.e., to be
awarded two difficult degrees (management and ChE) in four
years. While at Purdue, Gavin published his first paper, on the
analysis of drug release mechanisms from swellable polymers.
The paper was published in 1983 in the Journal of Membrane
Science and has been cited more than 200 times.
Not unexpectedly, the question of graduate education came
up very early in Jennifer's undergraduate studies, and the pos-
sibility of an academic career had appeal as early as 1980. So,
that fall Jennifer started working with Nicholas Peppas on a
project that led to a two-year collaboration, her first two pub-
lications, and her first exposure to fluid mechanics and mass
transport. The project involved the development of the first
advanced models of non-Fickian and anomalous transport in
glassy polymers with consideration of the relaxational effects
Vol. 42, No. 1, Winter 2008

Jennifer in third
grade, left, and with
her dad, Jerry Dolo-
resco, when she was
finishing high school
in 1979, below.
Opposite page: Jen-
nifer with fellow
researchers in her lab
at the University of

S 7------- -- S AN -- I' - 'I
of the macromolecular chains and the associated anomalous
transport. The problem had been addressed by others with
the addition of a pseudo-convective term (resulting from the
relaxational effects) until Jennifer's model set the physics
correctly, recognizing the importance of the moving bound-
ary value problem (glassy/rubbery transition and erosion
fronts). The paper, "Anomalous Transport of Penetrants in
Glassy Polymers," published in Colloid and Polymer Science
almost 25 years ago [261,404-409 (1983)] has received 180
citations and is one of the most cited articles in the history of
this journal. Another paper was published in the Proceedings
of the International Congress ofRheology of1984, where the
work had been presented.
Planning for graduate school included some long evening
discussions between Jennifer and Gavin. Jennifer had received
an NSF Fellowship and was clearly interested in a university
with a scholarly reputation where she could continue working
on transport phenomena. The choice came down to Princeton
University or the University of Pennsylvania-where Doug
Lauffenburger had put together a great academic and re-
search program and was particularly successful in attracting
a wonderful group of exceptional women chemical engineers.
Finally, Jennifer picked Princeton, but Doug has remained a
great supporter and mentor throughout her life.

Jennifer started at Princeton University in 1983, while
Gavin was employed a state away-at Air Products in Al-
lentown, PA! They decided to live in Flemington, NJ, which
meant that every day Jennifer had to drive 50 minutes each
way to school at Princeton. In the demanding classes there,
she met a number of other successful chemical engineers who
have remained good friends, including Paul Johnson (now
dean at Arizona State University), Harry Ploehn (now at the
University of South Carolina), and Spyros Anastasiadis (now
at the University of Thessaloniki, Greece). At Princeton, she
had many inspiring teachers such as Carol Hall, Bill Russell,
and Roy Jackson. Her rapport with the latter was exceptional;
he became the main advisor of her Ph.D. thesis, entitled "Ver-
tical Transport of Gas and Solids with Radial Solid Density
Variations." It was published-among others-in a seminal
issue of the AIChE Journal in 1989 ["Gas-Particle Flow in
a Vertical Pipe with Particle-Particle Interactions," AIChE J.
(1989) 35, 1473-1486]. It's one of the 50 most cited papers in
the history of the journal.
While she was in graduate school, Jennifer and Gavin
decided to start a family. Nicholas Peppas remembers, "In
October 1985 Jennifer called me with the good news she
was expecting her first child. And it was just two months
later that she called me back to share the terrible news that
Gavin had been diagnosed with cancer, probably in the lungs,
although that was not clear yet. We spent several days and
many hours talking about doctors and possibilities and we
ended up with the solution of Sloan-Kettering, where Gavin
was indeed saved after more than six months of treatment. I
will never forget that in all these months, with her graduate
research, her commuting almost two hours a day, and her
being pregnant, Jennifer did not sound weak even for one
minute. I never heard or saw her cry, complain, or give up.
This was an incredible lesson for me and for all of us around
her." Indeed, Jennifer was able to face these family crises
because of her deep faith and Christian beliefs that have given
her strength all her life.
Fortunately, 1986 was a much better year! A daughter, Jen-
nette, was born in April and the couple was reunited after the
completion of Gavin's surgery and many treatments. Gavin
went on to write a book, All Things Workfor Good (based
on Romans 8:28), about his cancer experiences, focusing on
the many lessons he learned and how his Christian faith was
strengthened. As for Jennifer, these were the days when she
began to master the art of juggling family, career, and life.
Her major professor in those days, Roy Jackson, a member
of the Royal Society, is now retired and lives in North Caro-
lina. He remembers, "As you may know Jennifer was
working for her Ph.D. in my group at Princeton at the
time her husband, Gavin, was diagnosed with cancer
and a very grim prognosis. At the same time she was
pregnant with their first child, so the pressures on her

Winning the 1989 Lafayette College Teaching Award
(pictured with students Paul Verderber
and Dave McVeigh).
were enormous. In addition she was working on the
problem of predicting, in detail, the behavior of gas-
particle flows in vertical tubes. The behavior of these
flows turned out to have an unexpectedly complex (and
interesting) structure, whose details were revealed only
by meticulously careful computation and interpretation.
In this she succeeded brilliantly, despite everything else
weighing on her mind at the time. Her resilience and
cheerfulness were simply amazing."
He continued, "After graduation the demands of her career
still had to be balanced against difficulties associated with
Gavin's health. But, in spite of this, her reputation grew steadi-
ly and I was proud to follow her progress. She is an amazing
woman who combines academic ability and administrative
competence with a warm and friendly personality that has
made her so many friends. Of my many academic 'children'
she is, perhaps, the one in whom I take the greatest pride."

In August 1987, while still a graduate student at Princeton,
Jennifer answered the calling to teach. As the family could
not leave the area because of Gavin's medical insurance (the
dreaded "pre-existing condition"), she accepted an offer from
Lafayette College and began teaching as an assistant profes-
sor in that school's Department of Chemical Engineering.
She stayed there until December 1989, one semester after
her formal graduation from Princeton. During her years at
Lafayette she became the master teacher and educator that
she is, winning numerous teaching awards.
As an educator, she taught a wide range of courses such as
Separations, Transport, Transport Lab, Fluid Mechanics, Heat
and Mass Transfer, Control, Fortran, and Differential Equa-
tions. To these courses she later added Kinetics, Advanced
Transport (grad level), Advanced Math (grad level), Fluid
Mechanics, and Heat and Mass Transfer at Carnegie Mellon;

Chemical Engineering Education

Jennette reads Daddy's book, All Things Work for Good,
to Derek (early 1996).

Kinetics, Fluid Mechanics, and Control at the University of
Arizona; and Materials and Energy Balances, Particulate
Systems, and Transport Phenomena at Purdue University.
By 1989, Gavin's medical condition had improved and he
was able to move to a new location-Pittsburgh-where he
was employed by PPG Industries. Jennifer interviewed and
was offered a position as an assistant professor at Carnegie
Mellon University. Nicholas Peppas remembers again: "She
was truly excited with the opportunities at CMU. In December
1989 she flew to Indianapolis and we spent three wonderful
days in our home-Jennifer, Lisa, and I, planning the future,
defining new problems, and talking about possible collabora-
tions on biomedical applications of micro- and nanoparticles.
I could see the future leader of particle technology in all her
intellectual glory."

At CMU, Jennifer found a truly exceptional faculty with
colleagues who provided the scholarly environment for her
professional development. During the interview process she
was particularly affected by then Head of the Department John
Anderson (now president of IIT in Chicago) as well as by
colleague Ignacio Grossmann. Anderson remembers, "When
I first met Jennifer during her interview at Carnegie Mellon,
I knew she was not only an excellent scholar and teacher but
also a future academic leader. She has proven me correct.
She leads by example, and her enthusiasm is infectious. The
way Jennifer has handled her personal and professional lives
provides a role model for all of us."
In a short amount of time Jennifer developed new courses
and attracted first-rate research students, such as Eduardo
Bolio, who is currently working for McKinsey in Siberia,
Vol. 42, No. 1, Winter 2008

and Christine Hrenya, who is an associate professor at the
University of Colorado in Boulder. The students helped her
develop a strong research program supported by industry and
government. Ignacio Grossman remembers that "Jennifer was
the strongest and most popular teacher in the department, and
received the highest scores on student evaluations. This is
especially noteworthy because she was teaching core courses
like Mathematics for Chemical Engineers and Heat and Mass
Transfer." She also started to give invited seminars and re-
search presentations outside the department while engaging
herself in university and professional service. While chair
of the undergraduate committee, she conceived a plan and
raised funds for a new computer cluster for undergraduates.
The computer cluster is still a very popular and well used
place for students to work and meet. While the computers
have been upgraded many times since, credit for the initial
concept goes back to Jennifer's ingenuity.
"Jennifer managed to do all of this, not only extremely
well, but always with a big smile, and while placing the ma-
jority of her attention on her husband and children," writes
Erik and Genevieve Ydstie moved to Pittsburgh in 1992
when he joined the CMU faculty, and they made a point of
moving into the same neighborhood as Jennifer and Gavin.
"We have a daughter in the same age group as Jennette and
we wanted to spend time together," Erik Ydstie remembers.
"We knew we could not go wrong as far as location was
concerned-Gavin and Jennifer had done extensive research
and had even written a computer program to take as many
considerations into account as possible to help in their hous-
ing-decision process! We have many fond memories of good
family times spent together after church and late evenings
discussing politics and all manner of other things."
Yet in the mid-1990s, after Jennifer received promotion to
associate professor at CMU, happy times gave way to worries
once again. There was some deterioration in Gavin's health,
especially due to the cold winters in Pittsburgh. Grossmann
notes, "I still recall the day that Jennifer walked into my office,
closed the door, and let me know she was leaving for Arizona
because Gavin's health required that he live in an environ-
ment with dry air. Despite the fact that her departure was a
big setback to our department, I admired how she was able to
put her family's priorities above her professional ambitions.
This to me speaks volumes of Jennifer's integrity, and her
great dedication to her family and her profession."
In 1995, the family moved to Tucson, where Gavin con-
tinued to work for PPG remotely and Jennifer became an
associate professor at the University of Arizona-whose
ChE department was under the expert leadership of Tom
Peterson, now dean of Engineering. The relocation rounded
out happily with the birth of Jennifer's and Gavin's son Derek
in early 1996.

But while Tucson was a great place for Gavin's health, and
Jennifer's career continued flourishing, the remote work with
PPG could not continue and medical insurance problems re-
surfaced. It was extremely difficult for Gavin to get insurance
from a new company due to his pre-existing condition. And
that's when Gavin and Jennifer remembered their roots.

In summer 1997, Jennifer and Gavin returned to Purdue to
launch the next phase
of their careers - she
as a tenured associate
professor and he as an
assistant professor in
organizational leader-
ship and supervision
as well as in agricul-
tural economics. Al-
though he was physi-
cally limited due to the
many health problems
he faced, Gavin had -w -
completed a master's
degree in economics
at Lehigh University
while working full time
at Air Products, and Above, the Sinclair family (Jenn
then completed a Ph.D. 1998 Caribbean cruise. An annu
in social and decision Jennifer and her family. Initially
sciences at Carnegie- then they became Alaskan and r
Mellon while working as well. Below, Gavin Sinclair's
full time at PPr which was featured in Newswe'
full time at PPG.
front page of The Wall S
This dual-academic-
career couple immedi-
ately set about having an impact on their respect ti i
departments and students. In spring 1998, Jennil .
launched a survey course in Particle Technolo,_'
(described in a CEE article in fall 1999). Begin-
ning in 2000, this course was video broadcast iti
include industrial audiences, and it was co-taughi t
with Carl Wassgren (M.E.) from 2002 on. The ,
initiative has continued to flourish, with Mike 0~i
Harris taking up the course when Jennifer no
longer could.
In 1999, Jennifer was elected trustee of
CACHE (Computer Aids for Chemical Engineel -
ing), where she championed the introduction of compuiduioial
fluid dynamics (CFD) into chemical engineering teaching (see
her CEE article in spring 1998). It was through her initia-
tive that Rodney Fox and Dick LaRoche were elected to the
CACHE board to further promote this domain in chemical
engineering research and education.
In fall 2000, while Gavin launched his book, Life, Love and

Economics, based on his very successful introductory econom-
ics course, Jennifer launched the administrative component
of her rapidly accelerating career: She accepted the position
of head of Freshman Engineering, leading a department that
was home to some 2,200 beginning engineering students.
She took on this major challenge with her usual optimism
and can-do style, re-energizing the staff, recruiting new fac-
ulty, and rallying the first-year students to overcome that big
transition to college. She was remarkably adept at soothing

. . . --. . .- _
ette, Jennifer, Gavin, and Derek) on a
al cruise has remained a tradition for
/ they were all Caribbean cruises, but
nore recently Mediterranean journeys
textbook, Life, Love and Economics,
ek magazine in May 2000 and on the
street Journal in October 2000.

the concerns of parents
while encouraging their
sons and daughters to
work hard but also find
time for campus activi-
ties. She was also very
effective at both attract-
ing new resources and
capitalizing on those
the department had.
She worked closely
with an Industrial Advi-
sory Council for Fresh-
man Engineering that
helped to refocus that
department's strategic
directions to its central
role of preparing stu-
dents for success in the
professional schools.
The chair of that coun-
cil, Bob Davis (Air
Products & Chemicals)
comments on her effec-
tive leadership style:

"Jennifer added energy to everything she
touched. Her real gift was empowering
all the members of her department so that
the resources of the department were seem-
ingly multiplied. The impact on results and
morale was remarkable!"
Sadly, into this whirlwind of positive pur-
'ntits came a devastating downturn: Gavin's
h .alth deteriorated again, and he was moved
Iick to Tucson because of the excellent
im.dical facilities there. He passed away during
S'linstmas 2000.

Jennifer's optimistic outlook prevailed. In spite of her many
administrative duties and the never-ending meetings, she
continued growing her research group to 10 solely orjointly
advised graduate students. In Fall 2001 she was promoted to
full professor and shortly thereafter (March 2002) she was
named the associate dean of the College of Engineering with

Chemical Engineering Education

responsibility for Undergraduate Education, while continuing
her duties as head of Freshman Engineering. And, of course,
the particle technology course had to go on as well! In her
associate dean role, Jennifer was instrumental in setting the
stage for the reorganization of the Freshman Engineering
Department and the subsequent launch of a new Department
of Engineering Education, the first in the United States. Other
challenges included restructuring and restaffing the Minority
Engineering Program, the Women in Engineering Program,
and the Cooperative Education Program. Her highly effective
leadership in engineering education was recognized through
the 2003 Sharon Keillor ASEE Award for Women in Engi-
neering Education. By fall 2003, it was clear that Jennifer
was ready to move on to further challenges-to the good
fortune of the University of Florida. Without question, her
Purdue years spanned a very productive phase of her profes-
sional life as teacher, researcher, and academic leader, but
also involved very low as well as high points in her personal
life. She remains a valued colleague and has certainly left her
impact on Purdue Engineering.

Throughout her academic career, the impetus for Jennifer's
energy and creativity has been her attraction to, and thus
seminal contributions toward, key research issues. In this area,

she has contributed on manifold fronts both on a national and
international scale.

Jennifer has an internationally recognized research program.
Her work in the development and validation of numerical
models for the prediction of gas-solid flow phenomena is in
an interdisciplinary research area of both great technologi-
cal and commercial importance. Particle flow processes are
encountered in the pharmaceutical, chemical, mining, agri-
cultural, food processing, and petroleum industries. Jennifer
has received various awards for her pioneering work in this
area. One of her most notable successes is the adoption of
her multiphase flow models by the two key commercial CFD
software-package vendors. Both companies-Fluent and
CFX-have adapted her models for dilute and dense-phase
gas-solid flow in software that is used internationally by both
academic and industrial experts.

Jennifer received the NSF Presidential Young Investigator
Award and an engineering college research award when she
was a professor at Carnegie Mellon University. Also early
in her career, she was awarded NSF equipment grants that
funded the purchase of important equipment such as a laser
Doppler velocimeter-a purchase that helped her propel her
early research program further toward the leading edge of


1. John Doney ChE, M.S., Carnegie Mellon Uni-
versity, "Heat Transfer in Dilute Particle-Laden
Turbulent Flows," 1993
2. Eduardo Bolio, ChE, Ph.D., Carnegie Mel-
lon University, "Dilute Turbulent Gas-Solid
Flows with Particle-Particle Interactions,"
3. Christopher Alexander, ChE, M.S., Carnegie
Mellon University, "Pneumatic Transport with
a Bimodal Particle Size Distribution," 1995
4. John Doney, ChE, Ph.D., Carnegie Mel-
lon University, "Electrostatic Behavior of
Charged Particle Mixtures," 1996
5. Christine Hrenya, ChE, Ph.D., Carnegie Mel-
lon University, "Dense Gas-Solid Suspension
Flow,"1996; 1997 recipient of Best Ph.D.
in Particle Technology by the International
Particle Technology Forum (Dr. Hrenya is
now associate professor of chemical and
biomolecular engineering at the University
of Colorado)
6. Pawan Agarwal, ChE, M.S., The University
of Arizona, "Modeling Gas-Solid Flow with
a Bimodal PSD," 1997
7. Michael Schabel, Materials Science, M.S.,
The University ofArizona, "Characterization
of Particle Traps in Plasmas," 1997
8. Tim Mallo, ChE, Ph.D., Carnegie Mellon
University, "Modeling Heat Transfer in Dilute
and Dense Phase Gas-Solid Fow," 1997
9. BerendvanWachem, ChE, Ph.D., DelftUniver-
sity, "Flow in a Fluidized Bed with a Bimodal
PSD," 2000 (Dr. van Wachem is now associate

professor at Chalmers University, Gothenburg,
Sweden, Department of Applied Mechanics)
10. Cynthia Blake-Powell, ChE, M.S., Purdue
University, "Comparison ofAd-Hoc Theory
Predictions and CFD Predictions for Con-
fined, Recirculating Flows," 2000
11. Agus Sumantri, ChE, M.S., Purdue Universi-
ty, "Comparative Analysis of Computational
Fluid Dynamic Models for Gas-Solid Flow
in Risers," 2000
12. Edward (Nick) Jones, ChE, Ph.D., Purdue
University, "The Effect of PSD on Gas-Solid
Fows," 2001
13. Chie Min Chung, ChE, M.S., Purdue Univer-
sity, "CFD Modeling of a Fluidized Bed with
a Wide Particle Size Distribution," 2002
14. Stephanus Budilarto, ChE, Ph.D., Purdue
University, "Experimental Study of Velocity
Ratio, Particle Size and Size Distribution
Effects in Particle-Laden Jet Fow," 2003
15. Kim Hayden Henthorn, ChE, Ph.D., Purdue
University, "The Effect of Particle Rough-
ness on the Fow Behavior of Fine Powders,"
NSF Graduate Fellow, 2004 (Dr. Henthorn
is now assistant professor in chemical and
biological engineering at University of Mis-
16. Kunn Hadinoto, ChE, Ph.D., Purdue Uni-
versity, "Dense-Phase Fow Measurements
using LDV," 2004 (Dr. Hadinoto is now
assistant professor in chemical and biomo-
lecular engineering at Nanyang Technologi-
cal University)

17. Michael Lasinski, ChE, Ph.D., Purdue
University, "Particle Clustering in Granular
Flows," 2004 (co-advised student)
18. Bill Ketterhagen, ChE, Ph.D., Purdue
University, "DEM Simulations of Granular
Discharge from Hoppers," 2006 (co-advised
19. Robert Hamilton, ChE, Ph.D., Purdue Uni-
versity, "Evolving Particle Size Distribution
in the Flow of Gas-Solid Mixtures," (1999 -)
(co-advised student)
20. Benjamin James, ChE, Ph.D., University of
Florida, "Effect of Particle Shape on Particle-
Phase Stress," (2004 -)
21. Mark Pepple, ChE, Ph.D., University of
Florida, "LDV Measurements of Liquid-Solid
Flow," (2004 - )
22. Anshu Anand, ChE, Ph.D., University of
Florida, "Effects of Particle Cohesion and
Shape on Particle Segregation in Hopper
Flows," (2005 - )
23. Julio Cesar Castro Vazquez, ChE, Ph.D., Uni-
versity of Florida, "Effect of Wall Roughness
on Wall Friction Angle," (2006 - )
24. Juan Pedro Marval, Departamento de Ener-
getica, Universidad Nacional Experimental
Francisco de Miranda, Punto Fijo, Edo.
Falc6n, Venezuela, "Numerical Simulation
of Saltating Particles using Granular Kinetic
Theory," (2006 - ) (co-advised student)
25. Byung Chu, ChE, M.S., University of Florida,
"Experimental Measurements of Sifting
Segregation in Hoppers," (2005 -)

Vol. 42, No. 1, Winter 2008

Right: At Carnegie Mellon University with
John Anderson and Ignacio Grossman in 1995.
Below: In Australia with a friend in 2005.
Facing page: Jennifer, Derek, and Barry celebrating
Jennette's graduation from her parents'
alma mater, Purdue University, in May 2007.

The importance of her research has been validated through
multi-year funding from major corporations such as Chevron,
DuPont, Alcoa, Dow, and Pfizer. Jennifer consults for a wide
range of organizations. She has also served as one of five aca-
demic participants in a U.S. Department of Energy-Office of
Industrial Technology Consortium on Multiphase Phase Flow.
This consortium involved collaboration between industry,
academia, and national labs. It received national recognition,
and as part of it Jennifer gave a panel presentation before
Congress in 1999.
Jennifer is a popular speaker not only for scientific but also
inspirational talks. In her academic career, she has given more
than 180 invited lectures.
Over the years, she has served the scientific community
in various positions. For example, she served on the Board
of Directors for the American Chemical Society-Petroleum
Research Fund (ACS-PRF), as a member of an NSF-orga-
nized delegation to South Africa promoting collaboration
in particle technology, as a member of an external review
panel to review the entire Chemical and Transport Division
of the NSF, as a member of the editorial advisory board of
the Journal of Powder Technology and the Journal of Phar-
maceutical Development and Technology, and since 2007 as
associate editor of the AIChE Journal. Within AIChE, she
served as the national meeting co-chair, along with Gintaras
Reklaitis, for the 2002 Annual AIChE Meeting that was held
in Indianapolis. She was elected to the executive committee
of the Particle Technology Forum for eight years and to the
AIChE Fluid Mechanics Programming Committee for four
years, and was chair of the Transport and Energy Processes
Division in 2006. A unique contribution to the annual chemi-
cal engineering meetings has been her continuous leadership
and support (with Rosanne Ford of the University of Virginia)
of the Christian ChE Fellowship meetings that are held on
Tuesday or Wednesday mornings of the AIChE meetings.
These meetings were initiated in 1989.


Jennifer is a highly effective teacher and mentor. She re-
ceived departmental and university teaching awards while at
Lafayette College and the University of Arizona. At Purdue,
she received the departmental Kimberly Clark Mentoring
Award and was named to the university-wide Teaching for
Tomorrow Program. In 2001, Jennifer served, along with
then Purdue Engineering Dean Linda Katehi, on a 10-mem-
ber invited NSF panel on the Future of Engineering Educa-

tion. From 2003-2006,
she also served on the
NAE's Engineering
Education Committee.
She is very active in
training new chemical
engineering professors.
She has participated as
an instructor for new
chemical engineering
educators at the Ameri-
can Society for Engi-
neering Education's
Chemical Engineer-
ing Summer School,
held every five years.
In 1997, she served as
a lecturer at this work-
shop. In 2002 and 2007, she organized and conducted a very
well attended workshop on CFD.
In August 2000, when she was named head of Purdue's
Freshman Engineering Department, she embarked upon
attending a series of workshops that better prepared her
for administration. For example, in 2001 she attended the
two-week Management Development Program at Harvard
University and in 2002 the American Council on Education's
Chairs and Deans Workshop.
As a Ph.D.-student advisor, Jennifer has excelled by trans-
ferring to her students the same style and the same warm
concern that Roy Jackson, her Ph.D. advisor at Princeton,
had transferred to her. Her former Ph.D. student Kimberly
Hayden Henthorn-now an assistant professor in chemical
engineering at the University of Missouri-Rolla-remarks,
"I think the thing that impressed me most about Jennifer
was her leadership style--she led by example. She kept her
graduate students motivated and working hard by showing us
that she was working just as hard as (and sometimes harder
than) we were."

Chemical Engineering Education

- -1
Naturally, Jennifer has not only been an exceptional
engineer, researcher, and educator, but also a wonder-
ful role model for the young generation of women
chemical engineers. More than 20 years ago, she had
the fortune to find a great role model in distinguished
ChE and NAE member Carol Hall, then at Princeton.
And many years later, Jennifer is now recognized as
an example for others. Again, Kim Hayden Henthorn
notes: "Having a strong, female role model like Jen-
nifer has made a tremendous impact on me. She has
been through so much, but has kept an upbeat, positive
attitude throughout. When I feel overwhelmed trying to
balance everything on my plate, I always remind myself
that Jennifer has been through all of this before, and
I give her a call. She is still willing to listen and give
good, thoughtful advice about everything from family
life to being a woman in academia."
Christine Hrenya, another former Ph.D. student of
Jennifer's and now an associate professor at the Uni-
versity of Colorado, admires Jennifer's ability to focus
and be successful. "Jennifer is incredibly productive,
but she has an unfair advantage over the rest of us -she
only needs four to five hours of sleep per night, and
this is her longtime average! Along the same lines, she
is an expert at multitasking. Once when I was visiting
her house for dinner, she promptly pulled out the iron-
ing board when dinner had finished (though perhaps
I was overstaying my welcome!)." This is the love
and admiration shown by her students. And Hrenya
concludes: "On a more serious note, Jennifer was an
outstanding advisor, and continues to be a role model,
friend, and sounding board for all topics under the sun.
I am extremely fortunate to have spent so many of my
formative years under her guidance."


In 2001 Jennifer met Barry Curtis, another Purdue chemical
engineer, who had graduated two years before her (in 1981), hav-
ing done the five-year co-op program. Jennifer comments, "I met
Barry for the first time at a Purdue Chemical Engineering Industrial
Advisory Board dinner in Houston. Jeff Hemmer then hooked us
up on a strategic planning project for Freshman Engineering, when
I was head of that department."
The rest is history. They were married in 2002 and moved to
Gainesville in 2005. Barry has been a wonderful father to Jennette
and Derek and a most supportive and admiring husband and partner
to Jennifer.
Just before the move to Florida though, the family took yet another
trip-to Australia! Jennifer was a visiting professor at the Chemi-
cal Engineering Department of Monash University in Melbourne
in 2005, received the Eminent Overseas Lectureship Award from
the Institution of Engineers in Australia, and visited a number of
Australian universities for seminars.
In 2005, Jennifer Curtis assumed a new responsibility as the chair
of the Department of Chemical Engineering at the University of
Florida. In a very short time, she has taken an already nationally
recognized program (achieved under the continuous and imaginative
leadership of the editor of this journal, Tim Anderson) and further
improved its national ranking.
Here is how Anuj Chauhan, associate professor in this depart-
ment, sums up the thoughts of the younger generation of professors
there: "Jennifer runs the department in a very fair, transparent, and
efficient manner. She is helpful to the entire faculty, particularly the
juniors, and she deserves a lot of credit for the rapid growth of the
department and the improvement in rankings."
At the University of Florida Jennifer keeps an active research
program and collaborates with the nationally known ERC program
on particle technology.
Meanwhile, the Sinclair Curtis family is thriving. Derek is now 11
years old and in sixth grade. Jennette attended Purdue in the School
of Management-the same school from which her father graduated
24 years earlier-and received her B.S. in management in May
2007 with highest distinction (a perfect 4.00). She is now a graduate
student at Michigan State University, specializing in organizational
behavior. Like her parents, she wants to become a professor.
Looking back at Jennifer Curtis's life and accomplishments, one
cannot help but admire her concern for others, her high standards,
and her dedication to education and advising. But one cannot forget
the adversities she has faced in her life. As Nicholas Peppas likes
to point out: "At some time or another all of us have had crises or
adversities in our lives. How Jennifer was able to cope with them
has been an inspiration for all of us. I think Jennifer, Barry, and
the children have found the golden rule of a good, inspiring, and
productive life. They are a wonderful example that all ;lii,,o work
for good." 0

Vol. 42, No. 1, Winter 2008

rj[ =department
---- U s_____________________________________

ChE at

Tufts University

The atrium of the Tufts Science and Technology Center, home to Chemical and Biological Engineering.

Tufts College was chartered in 1852 as the 163rd
institution of higher learning in United States. Two
years later, its first class matriculated. Today, Tufts
University is a diverse ensemble of undergraduate, graduate,
and professional schools that includes the Schools of Arts and
Sciences, Engineering, the Cummings School of Veterinary

Medicine, the Schools of Dental Science and Medicine,
the Sackler School of Graduate Biomedical Sciences, the
Friedman School of Nutrition Science and Policy, the Tissue
Engineering Resource Center, and the Fletcher School of
International Law and Diplomacy. Full-time students total
8,830, of which 4,950 are undergraduates.

� Copyright ChE Division ofASEE 2008
Chemical Engineering Education

The School of Engineering includes the Department of
Chemical and Biological Engineering as well as the Depart-
ments of Biomedical Engineering, Civil and Environmental
Engineering, Electrical Engineering and Computer Science,
and Mechanical Engineering.
Chemical Engineering debuted at Tufts in 1898 to become
only the fifth such program in this country, and accepted its
first students in 1900. The department was officially part of the
Chemistry Department until February 1949; its undergraduate
program was initially accredited three years later. An M.S.
program was added in 1962 and a Ph.D. program in 1965.
The first M.S. degrees were granted in 1964 to Dewey Ryu
(now a professor in the department of Chemical Engineer-
ing and Materials Science at the University of California,
Davis) and Wayne Sanborn (who worked for Matheson Gas
Products). The first Ph.D. degrees were awarded in 1971 to
Donald Merchant, Joseph DelPico (who worked for Polaroid
Corporation), and Edward Denk, Jr. (now a consultant at
Evergreen Solar).
This year the department will award some 15 graduate
diplomas, including five Ph.D.s Though it will award only
20 B.S.ChE degrees, enrollment in each of the continuing
undergraduate classes exceeds 30 students.
A perennial high-tech incubation zone, the greater Boston
area has long been one of the nation's hubs for biotechnology
R&D. Within this community, numerous professionals are
Tufts Chemical Engineering alumni. Many remain actively
involved with the department, which offers teaching and
research programs in biological engineering and has contrib-
uted nationally to development of the relevant undergraduate
For the past 15 years, the department has been assembling
a core group of biological engineering faculty. Graduate
research is divided within three main areas-biological en-
gineering, energy and environment, and advanced materials.
The current plan is to focus on expanding clean energy and
environmental sustainability along with systems engineering
research and teaching.

The department is currently housed in the Tufts Science
and Technology (SciTech) Center, which was constructed
in 1989 with major financial support from the U.S. Depart-
ment of Energy. At that time, the Department of Chemical
Engineering was the only academic department housed at the
facility. The Tufts High Energy Physics Laboratory, Biotech-
nology Engineering Center, Electro-Optic Research Center,
and Laboratory for Materials and Interfaces (LMI) were also
housed at the site. Both the Biotechnology Center and LMI
were associated research facilities of the ChBE Department.
Later, through the generous support of Raytheon Corporation,
the Pollution Prevention Projects Laboratory was established

within the department's facilities. In addition to the Depart-
ment of Chemical and Biological Engineering, the complex
currently houses the recently established Department of
Biomedical Engineering and the Gordon Institute for Engi-
neering Management.

The department has a long-standing history relating to the
development of educational programs in bioengineering. In
1986, the then Chemical Engineering Department initiated
a Massachusetts-sponsored continuing education program
to train technical staff for biotechnology business. The Bio-
technology Engineering Certificate program, which began
as a series of hands-on, introductory, technical courses and
remains a vital link to local industry, has led to the develop-
ment of a graduate-degree program. Today, the department
offers graduate degrees in both chemical engineering and
biotechnology engineering, including course-only M. Eng.
degrees attractive to working professionals, as well as tradi-
tional thesis-based M.S. and Ph.D. programs.
In 2004 the department hosted a series of NSF-sponsored
workshops focused on the development of chemical and bio-
logical engineering curricula. Participants from academia and
local bio-industries helped to shape our current undergraduate
program (see Table 1), which has evolved to include intro-
ductory courses on biology, statistical thermodynamics and
molecular structure (Physical Chemistry II), and biochemistry.
Simultaneously, the chemical engineering science core has
incorporated biological systems and processes, while the
upper-level course offerings have been supplemented with
applied biology-based electives (see Table 2). We are currently
redesigning the senior laboratory courses to comprise indus-
try-sponsored research projects in biotechnology, metabolic
engineering, energy sustainability, systems engineering, and
advanced materials.


Current Faculty
The current department faculty includes nine tenure-tracked
members and two jointly appointed members:
Professor Linda M. Abriola, currently dean of Engineering,
does research in the mathematical modeling of the transport
and fate of organic chemical contaminants in porous media.
She developed one of the first mathematical models to de-
scribe the interphase mass partitioning and nonaqueous phase
migration of organic liquid contaminants in the subsurface.
Her recent research involves the use of models and laboratory
experiments to examine abiotic and biotic processes influenc-
ing the persistence of organic and controlling the effective-
ness of aquifer remediation technologies. She is also a member
of the Department of Civil and Environmental Engineering,
and a member of the National Academy of Engineering.

Vol. 42, No. 1, Winter 2008

Professor Maria
nopoulos spear- 4
heads clean energy
technology and
catalysis research
in the School of
Engineering. Her
work focuses on Faculty member
catalytic fuel pro- Martin Sussman.
cessing, including
hydrogen generation, for fuel cells;
and new materials for air pollution
abatement and fuel gas desulfuriza-
tion. Her group's fundamental research
on catalyst structure and function has
led to the identification of a family of
very active oxidation catalysts based Some of the
on Au, Pt, and Cu on cerium oxide. Kyongt
An exciting new development in the Bottom
hot gas cleanup area is the group's recent identification of rare
earth oxides that remove hydrogen sulfide to sub-ppm levels
via reversible adsorption. Flytzani-Stephanopoulos's research
group presently comprises one postdoctoral fellow, seven Ph.D.
and four M.S. research assistants, and several undergraduates.
Her research is supported by the NSF, the DOE, the Army
Research Laboratory, and industry. Since 2002, she has served
as the editor of Applied Catalysis B: Environmental.
Professor Christos Georgakis established the Systems
Research Institute for Chemical and Biological Processes in
2005. Modeled on the Chemical Process Modeling and Con-
trol Research Center at Lehigh University, the Institute fosters
interdisciplinary research in modeling, optimization, and
operation of chemical and biological systems and processes
in collaboration with leading pharmaceutical and chemical
companies. The institute's first initiative-the Industrial
Consortium on Batch Manufacturing of Pharmaceuticals and
Chemicals (BMPC)-began operation in 2006. The consor-
tium is a cooperative research program that will focus on
improving the development, design, and operation of batch
pharmaceutical and chemical processes, with the ultimate
goal of halfing the time it takes to design new processes.
Initial research areas encompass designing batch reactors
for multistep synthesis; enantiomeric and polymorphic crys-
tallizations; plant-wide simulations, scale-up, and process
sensitivities; and statistical process monitoring and process
control that will effectively contribute to the implementa-
tion of the FDA's Process Analytical Technology Initiative.
The consortium officially launched its research program on
Oct. 25, 2006.
Professor David L. Kaplan's group applies biopolymer
engineering to elucidate structure-function relationships,
including those of fibrous proteins (collagens, silks) and
polysaccharides. They apply physiological and genetic engi-

:hBE faculty. Top, L-R: Christos Georgakis, Blaine Pfeifer,
bum Lee, Jerry Meldon, Daniel Ryder, Hyunmin Yi.
L-R: Ramin Haghgooie, Nakho Sung, Eric Anderson.
neering approaches to control biopolymer chemistry in order
to then explore self-assembly and function from a materials
science and engineering perspective. The biopolymer sys-
tems are studied for potential biomaterials applications, such
as matrices or scaffolds to support cell and tissue growth,
for controlled release of pharmaceuticals, and as biosen-
sor platforms. Selective processing, surface chemistry, or
genetic engineering redesigns are employed to functionalize
biopolymer material surfaces or control morphology and
structure, in order to enhance biomaterial interactions with
inorganic and biologic interfaces. In combination with human
stem cells and complex bioreactor systems, the biopolymer
biomaterial systems are employed in the generation of func-
tional human tissues including bone, cartilage, and adipose
systems. These systems are studied both in vitro and in vivo
for tissue regeneration, as well as to serve as human disease
models with which to understand disease progression, and
subsequently for pharmaceutical discovery. Transport issues
are addressed via vascularization strategies in vitro and in vivo
to promote improved functional tissue outcomes. Kaplan's
lab collaborates with many colleagues within and outside
of Tufts. He also serves as the chair of the Department of
Biomedical Engineering.
Professor Kyongbum Lee's research in systems biology,
metabolic engineering research in systems biology, and meta-
bolic engineering and tissue engineering is aimed at better
understanding the chemical design of living systems through
the development and application of analytical platforms, mod-
eling tools, and other enabling technologies. He is especially
interested in characterizing and manipulating cellular metabo-
lism for a variety of basic and applied studies. His research is
also motivated by the growing need to find and evaluate new
therapies for obesity and type 2 diabetes. Current work focuses
on the metabolism of liver and white fat cells.
Chemical Engineering Education

Left, one of several students in
the senior projects lab who got
hands-on experience when the
department's distillation column
needed repair. Below, Danny
Pierre, a graduate student in Pro-
fessor Flytzani-Stephanopoulos's
rntnlxi. lanh

There is increasing evidence that the functions and struc-
tural features of a biochemical network are intertwined. Thus,
a useful approach to studying the design of these complex
systems has been to characterize their layout, or "wiring."
Lee's research aims to develop a modeling framework to
integrate the analysis of the structural and functional layout
of biochemical networks. For example, he has recently de-
veloped and applied a series of algorithms for the rational
decomposition of metabolic networks based on experimen-
tally derived reaction activity data.
In the context of metabolic engineering, Lee's group uses
mathematical models of metabolic pathways in conjunction
with experimental biochemistry techniques to quantitatively
study the metabolic basis of adipose tissue development and
growth. In addition to fundamental insights into adipocyte
metabolic regulation, they seek to obtain: 1) novel drug tar-
gets, for example enzymes that catalyze the key controlling
steps in adipocyte lipid accumulation; and 2) robust, sensitive,
and easily measured metabolite biomarkers for diagnosis of
obesity-related diseases and efficacy assessment of thera-
peutics. To these ends his group has performed in silico and
in vitro studies characterizing the metabolic flux profiles of
de novo adipocyte formation, identified target pathways for
metabolic intervention, and demonstrated the feasibility of
achieving a global shift in lipid accumulation through the
forced expression of a single protein.
Lee also uses tissue engineering to design, build, and
optimize model systems whereby important cellular func-
tions may be studied under well-controlled conditions while
Vol. 42, No. 1, Winter 2008

Above, Gregory Botsaris gives a
presentation at the 2005 sympo-
sium in his honor.

maintaining a high degree of physi-
ological relevance. The primary
focus here is on engineering the
cellular component and its medium
(as opposed to the materials or
reactor housing). Current projects leverage developments in
micro-fluidics and stem cell biology to generate flow-through
incubators that afford spatial control over the chemical en-
vironment of the cultured cells. This incubator is used for
a number of on-going studies involving heterogeneous cell
components (e.g., adipocyte-preadipocyte co-culture) and
solution gradients of exogenously added chemicals (e.g., drug
transformation by hepatocytes).
Professor Jerry H. Meldon's research interests are in mass
transfer and separation processes, particularly when revers-
ible chemical reactions control the rate and selectivity of a
diffusional separation. His primary focus is on membrane
processes, gas scrubbing, and mathematical modeling of mass
transfer with chemical reaction. Recent work includes theoreti-
cal and experimental investigations of hydrogen transport in
palladium membranes and its role in catalytic membrane reac-
tors for the production of pure hydrogen. One Ph.D. student is
building and modeling a novel spinning disk system for reaction
and separation; another is modeling NOx scrubbing.
Professor Blaine A. Pfeifer's research group seeks to
influence genetic, metabolic, cellular, and process events in
the synthesis of therapeutic products. For example, they ap-
ply molecular and process engineering in conjunction with
molecular biology, microbiology, analytical chemistry, and
bacterial genetics to develop microbial bio-processes and
products that can potentially target cancer, bacterial infections,
diabetes, and other diseases.
Pfeifer's group also seeks more efficient and economical
ways to generate biological products. One approach is to

Left, Zheng Zhou, a graduate student in Professor Flytzani-Stephanopoulos's catalysis lab. Right, graduate student Lisa
Schupmann works with a multi-reactor system in Professor Georgakis's Systems Research Institute.

transplant genetic material responsible for an important thera-
peutic product into a convenient and process-friendly bacterial
microorganism (such as Escherichia coli) for eventual product
scale-up and development. Current projects include: 1) cellu-
lar and metabolic optimization for production of the antibiotic
erythromycin; 2) production of new and established anticancer
agents; 3) development of new biological vaccine systems;
and 4) modeling and experiments to identify and elucidate
cellular bottlenecks in natural product biosynthesis.
Professor Daniel F. Ryder's primary research lies in the
processing and characterization of ceramic and organically
modified ceramic materials for electronic and optical appli-
cations. Specific areas of study include: Sol-gel processing
of PLZT ceramic thin films for applications in ferroelectric
memories; chemical processing of hexagonal ferrite thin films
for application in high frequency devices; chemical process-
ing of superconducting oxide materials; the development of
organically modified ceramic materials for abrasion resistant
and antireflective coatings on optical polymeric substrates;
and chemical vapor deposition methods for III-V nitride
Professor and Chair Nakho Sung's research in the area
of structure-property-processing relationships of polymers
and composites includes development of in situ monitoring
techniques for characterization of reactions in polymers and
composites via fiber-optic fluorescence and UV reflection.
Both extrinsic and intrinsic reactive fluorophores are devel-
oped for various matrix polymers for advanced composites to
follow cure reactions as well as water uptake using distal-end
and evanescent fiber optic probes. UV reflection spectroscopy,
coupled with bifurcated fiber-optic probes, is also developed
to characterize reactions in composites that are fluoresce-
insensitive, such as those based on polyimides, epoxies, and
unsaturated polyesters; this technique complements fluores-
cence monitoring.

Another area of research has been to develop a methodol-
ogy for tailor designing composite interphases. Using a low-
temperature plasma of various monomers and inert gases,
interfacial structure is modified to enhance such properties
as adhesion strength and long-term stability against tempera-
ture and moisture-induced degradation. Sorption/diffusion
behavior of organic molecules in multiphase polymers is
investigated using a model composite consisting of glassy
and rubbery domains where morphology is systematically
varied. These model composites facilitate elucidation of
complex diffusion behavior and development of a predictive
diffusion model. Recently, Sung's group has been developing
low-cost membranes for fuel cell application and hydrogen
purification. Current focus is on anionic membranes suitable
for low-temperature direct methanol fuel cells.
Professor Hyunmin Yi's research interests lie in biochemi-
cally driven nanometer scale fabrication (nanobiofabrication)
of high throughput biosensors, biophotonic devices, and
nanocatalysts for biomedical, environmental, and energy ap-
plications using smart biopolymers and viral nanotemplates.
His primary research areas are nanobiofabrication with geneti-
cally modified viral nanotemplates, and biophotonic device
fabrication with smart biopolymers.
Yi's group has developed nucleic acid hybridization-based
surface assembly strategies of genetically modified tobacco
mosaic viruses through their own genomic mRNA. These
biologically derived nanotubes - with the precise dimensions
of 300 nanometer (nm) length,18 nm outer diameter and 4
nm inner channel-can serve as nanotemplates for covalent
coupling of functional nanoparticles. Building upon these
facile assembly strategies of these potent nanotemplates, his
group is working toward the development of high-throughput
biosensors and BioMEMS devices for environmental and
biological threat detection. In closely related collaboration
with Flytzani-Stephanopoulos, they are also developing na-

Chemical Engineering Education

A close-up of the Tufts Chemical and Biological Engine
X-ray diffractometer.

noscale gold particle-based catalysts for energy applications
using the viruses as nanotemplates with high capacity and
precise nanoscale spacing.
In his second research area, Yi seeks to develop biocom-
patible high-throughput biosensing platforms and implant-
able biophotonic devices for biomedical and environmental
applications. They exploit the stimuli-responsive properties
of smart biopolymers to fabricate nanoscale patterns and
waveguides with high spatial, temporal, and orientational
control. Example biopolymers include structural proteins
such as gelatin and silk as well as polysaccharides such as
agarose and chitosan. Recently his research group has shown
that the thermo-responsive morphology transition of common
biopolymers such as gelatin and agarose can lead to efficient
means for consistent manufacturing of nanometer-scale sur-
face diffraction gratings under mild processing conditions.
They are currently working on building all-biopolymeric
implantable waveguides, photonic bandgap crystal fiber-based
biosensors and nanoimprinted biopolymeric gratings.
Research Professor Aur6lie Edwards's research falls broadly
in the domain of biological transport phenomena, developing
theoretical models of transport in the renal medullary micro-
circulation, both at the macroscale and the cellular scale. The
microcirculation in the renal medulla plays an essential role in
the regulation of fluid and electrolyte excretion and in the long-
term maintenance of arterial blood pressure. Understanding
how medullary blood flow is regulated and how it influences

in turn arterial blood pressure is
the main objective.
Specifically, a detailed math-
ematical model of transport in
the tubular and microcirculatory
systems of the renal medulla is
in development in order to ex-
amine how interactions between
oxygen, nitric oxide (a potent
vasodilator), and superoxide
affect tubular sodium reabsorp-
tion, medullary blood flow, and
blood pressure.
Mathematical simulations
of ion concentration changes
Within the bulk cytoplasm and
microdomains are under devel-
opment to elucidate contractil-
ity of resistance vessels, such
as medullary descending vasa
recta, regulated by variation of
intracellular Ca2+ concentration
([Ca2+]cyt) in both endothelium
ring Department's and smooth muscle. This model
accounts for the characteristics
of the channels and transporters that exchange ions between
the plasma membrane, SR stores, and cytosol. Edwards's
group is investigating which factors are most likely to serve
as principal determinants of [Ca2+]cyt and delineate the
mechanisms through which the bulk cytosol, microdomains,
and SR communicate. Their goal is to ultimately relate local
medullary NO concentrations to [Ca2+]cyt changes and to
blood vessel diameter variations.
Emeritus Faculty
The early transition years of any department present a series
of challenges that require the selfless dedication, vision,
and determination of its faculty. The accomplishments
of these early efforts are often taken for granted with the
passage of time, but whatever level of success that this
department presently enjoys is due in large part to a small
group of faculty that has served Tufts throughout their
professional career.
Martin V. Sussman joined the department in 1961 and
served on the faculty for 37 years before retiring due to ill-
ness in 1998. His popular "Technology as Culture" course,
which offered liberal art students insights into the world of
technology, made him a household name across the Tufts
community. He earned his doctorate in chemical engineer-
ing from Columbia University and was a Fellow of both the
American Institute of Chemical Engineers and the American
Chemical Society. He was the author of two books on ther-
modynamics as well as numerous technical articles. Martin's

Vol. 42, No. 1, Winter 2008

research interest spanned the fields
of materials science and biotechnol-
ogy. He was author to some 20 pat-
ents that included the "Incremental
Draw Process" for the processing of
synthetic fibers and the "Fibra-Cel"
tissue culture matrix.
The department was saddened by
his passing in 2005, but his legacy
and connection with Tufts survive
through the generous funding of the
Jeanne and Martin Sussman Endowed
Gregory D. Botsaris joined the
department in 1965 immediately after
earning his Ph.D. in chemical engi-
neering at M.I.T. He officially retired
in 2004 after 39 years of service, but
remains on campus to continue his
research in crystallization processes.
He is the author of numerous technical
articles in that field and other sepa-
ration processes. His contributions
to crystallization technology were
recognized by his peers through the
establishment of the Gregory Botsaris
Lecture series.
Ken Van Wormer is the new-
est emeritus faculty member. Van
Wormer joined the department as an
instructor in 1954 while completing
his doctoral studies at M.I.T. His
continuous service to the department
for an impressive 53 years-including
10 years as chair-is a testament to his
dedication to the profession.

Chemical & Biological Engineering Curriculum
at Tufts University
First Year
Fall Term Spring Term
EN 1 (Intro. to Computers in Engineering) EN 2 (Intro. to Engineering Communication)
Math 11 (Calculus 1) Math 12 (Calculus 2)
Chem 1 (Chem. Fundamentals w/ Lab) Chem 2 (Chem. Principles w/ Lab)
English 1 (Expository Writing) Humanities/Social Sciences Elective
EN Elective (Intro. to Engineering) Physics 11 (General Physics w/ Lab)
Second Year
Fall Term Spring Term
ChBE 10 (Thermo & Process Calc I) ChBE 11 (Thermo & Process Calc II)
Math 13 (Calculus 3) Math 38 (Differential Equations)
ES 11 (Intro. to Biology) ChBE 39 (Appl. Math and Numerical Methods)
ES 10 (Intro. to Materials Science) Chem 32 (Physical Chemistry II)
Chem 33 (Physical Chemistry Lab I) Humanities/Social Sciences Elective
Third Year
Fall Term Spring Term
ChBE 21 (Fluid Dynamics & Heat Transfer) ChBE 22 (Mass Transfer)
Chem 51 (Organic Chem I) ChBE 102 (Reactor Design)
Chem 53 (Organic Chem Lab I) Bio 152 (Biochemistry & Cellular Metabolism)
ES 3 (Intro to Electrical Engineering) Advanced Chemistry Elective
Humanities/Social Sciences Elective ChBE Concentration Elective
Free Elective
Fourth Year
Fall Term Spring Term
ChBE 45 (Separation Processes) ChBE 60 (Chemical Process Design)
ChBE 51 (ChBE Unit Ops Lab) ChBE 52 (ChBE Projects Lab)
ChBE 109 (Process Dynamics & Control) ChBE Concentration Elective
ChBE Concentration Elective Advanced Chemistry Elective
Humanities/Social Sciences Elective Humanities/Social Sciences Elective


Sources for this feature article came primarily from the
ChBE Department archives, the Tufts Office of Institutional
Research, and the Tufts Office of Alumni Relations. Professor
Kenneth Van Wormer's "Speech on the Occasion of the 85th
Anniversary of Chemical Engineering at Tufts, April 1986"
was an especially valuable source for the department's early
history. His speech contained references to Russell F. Miller's
Light on the Hill: A History of Tufts College (Boston: Beacon
Press, 1966). Authors acknowledge the help from Joanna
Adele Huckins in drafting this article. 7

Biological Engineering Electives
ChBE 160 Biochemical Engineering
ChBE 161 Protein Purification
ChBE 162 Introduction to Biotechnology
ChBE 163 Recombinant DNA Techniques Laboratory
ChBE 164 Biomaterials & Tissue Engineering
ChBE 166 Principles of Cell and Microbial Cultures
ChBE 167 Metabolic and Cellular Engineering
ChBE 168 Biotechnology Processing Projects Laboratory

Chemical Engineering Education

M] =1laboratory



in a Simple Experiment

Universidad Nacional del Litoral * Santa Fe, Argentina

Since its beginnings, chemical engineering education
has incorporated practical activities into the curricula
(known as laboratory practices) in varying degrees
(number of teaching hours, number of credits) according
to the type of course or level of advancement in a given
program. The historical evolution of this type of activity in
chemical engineering programs has been recently discussed
in literature.[1,2] It seems quite clear that, with the exception of
certain courses and subjects directly involved with standard
practical experiments (e.g., general and instrumental analytic
chemistry), there is ample choice in their design and imple-
mentation. On the other hand, new developments both in the
processing of experimental data and in process control and
automation techniques make it possible to design and imple-
ment laboratory practices in practically every area at any
level of sophistication and scientific accuracy. This requires
important investments in time and money that, unfortunately,
many schools are not willing to undertake.
A valid alternative, then, is to redesign these practical activi-
ties to maximize the use of resources. This work describes a
simple, inexpensive lab experiment that was originally de-
signed for qualitative demonstration on the existence of "hot
spots" in discontinuous reactors. It has been deeply modified
from the viewpoint of the processing and interpretation of
the data collected through the use of a relatively sophisti-
cated mathematical model and the independent validation
� Copyright ChE Division of ASEE 2008
Vol. 42, No. 1, Winter 2008

of certain parameters. With these modifications, the field of
knowledge has been broadened (development and validation
of mathematical models, heat transfer and sensor dynamics)
at practically null costs (considering time and economic
In what follows, both the experimental procedure and the
development and validation of the model are described. Re-
sults illustrate this type of activities performed by students.
The subject to which this activity belongs is Chemical Reac-

Maria Eugenia Sad is a Ph.D. student in chemical engineering. She
received her bachelors degree from Universidad Nacional del Litoral in
chemical and food engineering. Her research interests are in chemical
reaction engineering and heterogeneous catalysis by zeolites.
Mario Ricardo Sad is an assistant professor of chemical reaction
engineering. He received his bachelors degree from Universidad Na-
cional del Litoral in chemical engineering. His research interests are in
chemical reaction engineering and catalytic reactions of hydrocarbons
on metallic catalysts.
Alberto Antonio Castro is a full professor in chemical reaction engineer-
ing. He is dean of chemical engineering faculty and director of Institute
of Research in Catalysis and Petrochemistry (INCAPE). His research
interests are in chemical reaction engineering and catalytic reaction of
hydrocarbons by metallic catalysts.
Teresita Francisca Garetto is an associate professor of chemical reac-
tion engineering. She received her bachelors degree from Universidad
National del Literal and her Ph.D. degree from University of Zaragoza,
both in chemical engineering. Her research interest are in environmental
catalysis and catalytic reactions of hydrocarbons.

tion Engineering I, which also has other practical activities
of different design that range from pilot scale experiments
with semi-industrial equipment to computational simulations
of real processes.

Several authors have reported studies with adiabatic reac-
tors using the strongly exothermic reaction between sodium
thiosulfate and hydrogen peroxide in aqueous media. At first,
this reaction was used to experimentally assess theoretical
predictions about hot spots in discontinuous or plug flow reac-
tors and steady state multiplicity and hysteresis or oscillations
in continuous stirred tank or recycle plug flow reactors.[37]
More recently, this reacting system was used to evaluate some
new techniques in control and optimization engineering and
in risk management.[818]
Heat transfer in chemical reactors is an important topic
both from theoretical and practical points of view. A common
practice in chemical reaction engineering courses (which
may be taken in parallel with energy transfer courses by the
students) is to assume very simple models of this phenom-
enon to avoid further complications in the classical chemical
reactor models.
Another aspect added to thermal kinetic data processing is
the incorporation of the "time lag" concept in the temperature
measurement. This concept naturally arises in the course of
process control, which students take normally afterward, but
may be easily modeled as a simple, first order heat transfer
lag to help introduce students to the idea of non-instanta-
neous measurement. In this case, chemical reaction is very
rapid, and we have chosen a sensor (digital thermometer)
with a deliberately high thermal mass so the lag effect is
more noticeable.

The reaction between sodium thiosulphate and hydrogen
peroxide in aqueous solutions is highly exothermic (AHR ~
147 kcal/mol Na2S203) and practically irreversible. The cur-
rently accepted stoichiometry for the reaction is:
2Na2S203 +4H202 - Na2S30 + Na2SO4 + 4H20 (1)
2A+4B- C+ D +4E
The kinetic behavior of this system has been well studied
by several authors[3 6] and is first order with respect to each
of the reagents:
r = k.CA.CB [r]= gmolNa2S203 / cm3.s (2)
The dependence of k with T will be taken as the usual Ar-
rhenius form:

k =A.eE/RT (3)

even though some slightly different form was reported in
earlier studies.
For a batch reactor the mathematical model is made up of
two differential equations resulting from mass and energy
balances. Concentration of reagents sodium thiosulphate (C,),
hydrogen peroxide (CB), and temperature (T) can be described
by using these balances, the stoichiometry ratio between both
reagents, and the initial conditions:
dCA / dO = -k(T.CA .CB (4)
V.Cp,.(dT / dO) = k(T).CA.CB.(-AHR)- A,.U.(T - t) (5)
CB = CB0 -2.(CA- CA) (6)
T(0) = 0 CA(O) =AO CB(O) = CB (7)
The thermal mass of the beaker used as a reactor has been
neglected in the energy balance (Eq. (5)), due to its low mass
(less than 4g), to avoid further complications in the proposed
This set of differential equations can be solved using an
appropriate numerical method (Euler, Runge-Kutta, etc.) in
order to obtain the evolution of products concentration and
temperature with time. It must be noted that this system may
exhibit a "stiff'"17I18] behavior, so the numerical integration
method should be a robust one.[191
For exothermic reactions such as this one, with heat ex-
change resulting from an external cooling media, the tem-
perature vs. time profile has a maximum value ("hot spot").
At this point:
k(T).CA.CB.(-AHR) = A.U.(T - t) = dT / dO = 0 (8)
After this point the temperature will decrease with a di-
minishing slope, as the reaction rate tends asymptotically to
zero. At infinite time, the temperature within the reactor will
be the same as that of the surrounding media:
[k(T).C.CB .(-AHR)1- = 0 (9)
[dT /d '.. =(A,.U/V.Cp,).(T-t)= Toe = t (10)
This behavior at sufficiently high times is useful to "inde-
pendently" estimate the value of U (from At.U/V.Cpv) using
values of temperature/time in this time domain. For the vessel
used as a reactor, values of V and At were 100 cm3 and ca.
120 cm2. The value of Cpv was 1 cal/g.K.81 These values were
used to estimate the heat transfer coefficient.
The "time lag" of the sensor (digital thermometer) was taken
into account through the simple first-order model:[22 24]
dTR/ d = (1/ T).(T- TR) (11)
which may result from the finite heat transfer rate between the
fluid, which is at the "real" temperature T, and some internal
point in the sensor probe, which is at the read temperature TR.
Both the high chemical reaction rate and the deliberate choice
of a sensor with high thermal mass justifies the incorporation
of the new dependent variable TR in the model with its ad-
ditional differential equation and a new parameter (T).
Chemical Engineering Education

The equipment used is simple, inexpensive and is comprised
of the following elements: a 100 mL disposable, lightweight
polystyrene beaker, a magnetic stirrer, a digital thermometer
of the type used for frozen meat control at cold-storage plants
(273 K/423 K range, 0.1K resolution) and a 500 mL crystal-
lizer to be used as an ice bath. A few millimeters of clearance
was left between the bottom of the beaker and the stirrer plate
to achieve a more homogeneous external heat transfer envi-
ronment. The chemical reactants employed are of analytical
grade, but commercial hydrogen peroxide (10-20 vol, with
stabilizer) or commercial thiosulfate have been used with no
substantial result alteration.
A series of experiments is synthetically described below,
for typical work performed by students that can be easily
modified from operating conditions or sequencing.


Chemical reaction involving heat exchange with
the environment (kinetic experiment)
This experiment consists of loading the reactor with 50
mL of a Na, S, ~ 0.8M solution, placing the thermometer
in an adequate fixed position, and stabilizing the stirring
rate to finally add 50 mL of HO, ~1.2M as rapidly as pos-
sible. Then, temperature readings should immediately start
at appropriate time intervals. Figure 1 shows typical results
obtained by students conducting this experiment at 303 K
(room temperature).



. 340

- 330




0 10 20 30 40 50 60
1 Time (sec)
Figure 1. Reaction between Na2S20, and H202 in a batch reactor with
heat exchange by natural convection (room temperature: 300 K).
Vol. 42, No. 1, Winter 2008

It was observed that the temperature increased rapidly at
the beginning, went through a well-marked maximum for
about one minute, and then decreased increasingly at a slower
pace. This results were obtained because when the limiting
reactant is practically exhausted, the reaction rate becomes
negligible and the system behaves as a simple vessel with
fluid at a temperature close to room temperature, where it
exchanges heat by convection (internally forced, externally
natural) with low heat transfer coefficients.[2021]


Natural cooling of the vessel with no chemical
This experiment consists of following the liquid cooling
(preferably using the same solution used in Experiment #1
after its completion) under the same conditions previously
analyzed, from a temperature of about 353-363 K. The final
liquid from the previous experiment, preheated up to that
temperature level, is left to evolve spontaneously under the
same stirring conditions. Many time/temperature values
are then registered (10 - 15 values) to estimate a value of
the U global exchange coefficient under conditions similar
to those of the kinetic experiment. To do so, it suffices to
use the asymptotic integrated form of the energy balance
differential Eq. (5):
ln[(T0 -t) (T- t)]= (A,.U V.Cpv). = ca.0 (12)
The slope of a straight line in the graph of ln[(T - t)/(T - t)]
vs.6 is the value of a. Figure 2 (next page) shows the results
of an experiment of this type together with statistical
analysis of the data.


Determination of the thermometer time lag
(time constant)
This experiment should also be performed with
the liquid sample of the first experiment, preheated
to a convenient temperature level (343-353 K) and
existing under the same conditions (stirring, etc.).
In this case, the thermometer (previously in thermal
equilibrium at room temperature) is introduced and
fixed as quickly and as smoothly as possible and the
registered temperature evolution is read as a function
of time (given the ordinary response times of tools
such as the thermometer employed, the readings
should be carried out at a frequency of about 2-5 s).
An essential condition for this experiment is that the
temperature of the liquid mass be constant during the
experiment. Since the time consumed is short, this
condition is practically accomplished. The integrated
form of differential Eq. (11) yields the classical step
input time response for a first order system:

In[(TR - T,) /(T R T) = ( ] ( r).O

from which the value of time constant, i
by regression of the experimental data. Th
"final" value shown in the display when
(T, = TRI >). Figure 3 shows typical
results as well as the estimation of T.
Analysis of the data: Models and
interpretation of results
The results obtained with a typical
kinetic experiment may be correlated
using a simple model (Model I) inte-
grated by the mass balance (Eq. (4)) and
energy (Eq. (5)) differential equations
with proper boundary conditions and
parameters A, E and U (or a). They
can also be correlated by a more so-
phisticated model (Model II) including
the differential equation accounting for
thermometer time response with the ad-
ditional T parameter. On the other hand,
the number of parameters (and hence
the error associated with its determina-
tion) are reduced to 2 in both models
if values of a and T, obtained from the
analysis of the results of Experiments
# 2 and 3, are taken as known. Figure
1 shows the adjustment of both types
of models (with its associated floating
parameters), and Table 1 summarizes the
results obtained as well as the compari-
son between the parameters obtained
and those from independent experi-
ments or reported in the literature. No-
tice the better adjustment of the model
in Figure 1, including the time lag of
the thermometer readings. Likewise, in
Table 1, a good correspondence can be
observed between the parameters de-
termined from the kinetic experiment,
those from independent experiments,
and those reported in the literature.
This correspondence is more marked
in the case of the model including the
thermometer time lag. The interpreta-
tion and analysis of the results is usually
performed using any available software
capable of performing the regression
of the temperature-time experimental
values, by way of numerical integra-
tion of the differential equations (with
methods adequate for "stiff' systems)
using any parameter search engine


, can be estimated
ie value ofT is the
it varies no longer

such as simplex, steepest descent, Monte Carlo, etc. There
is a lot of general software available for these ends (Matlab,
Matematica, MathCAD, to mention just a few) and some
direct-oriented free applications (such as Scientist from
Micromath) to solve the problem of statistical determination

370 2.5

In [(T-Tmo)/(T-Tm)] = 0.001009.0 + 0.0034
360 , R2 0.9987
A.U/V.Cpv= 0.001009 - 2

350 -

-1.5 F7

" 330
E . 1
l- _ =


, * * - 0.5
310 *

300 . . 0
0 400 800 1200 1600 2000
Time (s)
Figure 2. Determination of parameter a defined in Eq. (12).

360 8

.--0 ..... 0 ..-... 4-- ...- ....-- - .
350 - "" .- . ..

-- 6
340 -
. 330 - . .
S- 4 ^
) : c -
: In[(TR - TF)/(TR -TF)]= 0.1259.08
310- R2 = 0.9973
0 1/ = 0.9973s, 2


290 q 0
0 10 20 30 40 50 60
Time (s)
Figure 3. Determination of the thermometer time lag.

Chemical Engineering Education

of the parameters and/or to simulate the system. It is not the
aim of this work to approach the wide field of mathemati-
cal and statistical methods employed in the correlation and
modeling of experimental data, so only the results are shown
as long as alumni are not compelled to use any particular
software. The analysis of data presented in this paper was
done using software with Episode (stiff) integrator, which
includes the analytical calculus of the Jacobian matrix and
a least square method (minimizing relative error) for the
parameter adjustments.
Testing the predictive power of the model: Ice
bath experience
An additional experiment to assess the usefulness of the
model to predict the behavior of the system can be done by
modifying the heat transfer scenario (i.e., using an ice bath as
external transfer medium). In this case the values of kinetics
parameters (A and E) and the sensor time constant are taken as
the previously determined ones, and the heat transfer param-
eter is separately determined in an individual cooling experi-
ence as described before (as heat exchange coefficient value
is expected to change for the ice bath surrounding media)
or estimated from temperature-time values for higher times.
Typical results of such experiences are not shown here, but in
Fig. 4 a good concordance can be seen between predicted and

0 20 40 60 80 100
Time (min)
Figure 4. Testing the predictive model in ice bath experience.

measured temperature vs. time profiles in this case, confirming
the usefulness of the model for predictive purposes.

Some very simple laboratory experiments can be inexpen-
sively modified to obtain higher pedagogical output through
a more open analysis and correlation of obtained data and
some minor operative modifications, including concepts and
methods that strictly pertain to another (but closely related)
area of knowledge. Modeling and simulation, two fundamen-
tal tools in chemical engineering, should be used despite the
modest quality of equipment and obtained data to illustrate
important aspects of the formulation and validation of mod-
els, design of complementary experiences for estimation of
parameters, and the utilization of current numerical analysis
and statistical procedures and software. Rather surprisingly
in this case, as long as only one variable is measured, the
results of parameter estimation and simulation are good and
constitute an additional stimulus for students.

A Arrhenius preexponential factor [cm3 mol /s 1]
A, Effective heat transfer area [cm2]
C Molar concentration [mol/cm3]
Cpv Volumetric heat capacity [cal/cm3/K 1]
E Activation energy [cal/mol 1]
k Kinetic constant [cm3 mol /s 1]
r Reaction rate [mol/cm 3s 1]
R Gas constant [cal/mol 1/K 1]
T Temperature of system, real [K]
t Temperature of surrounding media [K]
T Initial temperature [K]
T Constant temperature of medium in Experiment #3 [K]
TR Digital thermometer temperature reading [K]
TRO Initial digital thermometer temperature reading [K]
U Global heat transfer coefficient [cal/s 1/cm 2/K 1]
V Reactor volume [cm3]
0 Time [s]
T Thermometer time constant [s]
a Parameter defined as At.U/V.Cp [s 1]
AH, Enthalpy change of reaction [cal.mol 1]

1. Feisel, L.D., and G.D. Peterson, "The Challenge of the Laboratory
in Engineering Education," J. Eng. Ed., 91(4), 367 (2002)

Comparison of Parameter Values

A 2.1010 2.1010 - 2.1010 - 6.85.10"
E (Cal/mol) 8120 8238 - 8238-9200
U (KCal/s m2 K) 0.0031 0.0035 0.0036 0.002-0.02
T (s 1) - 0.1259 0.1259

Vol. 42, No. 1, Winter 2008




E 280

- Simulation
o Experimental

2. Feisel, L.D., andA.J. Rosa, "The Role of the Laboratory in Undergradu-
ate Engineering Education," J. Eng. Ed., 94(1), 121 (2005)
3. Root, R.B., and R.A. Schmitz, "An Experimental Study of Steady State
Multiplicity in a Loop Reactor," AIChE J., 15(5), 670 (1969)
4. Vetjasa, S.A., and R.A. Schmitz, "An Experimental Study of Steady
State Multiplicity and Stability in an Adiabatic Stirred Reactor,"AIChE
J., 16(3), 410 (1970)
5. Chang, M., and R.A. Schmitz, "An Experimental Study of Oscillatory
States in a Stirred Reactor," Chem. Eng. Sci., 30, 21 (1975)
6. Guha, B.K., G. Narsimhan, and J.B. Agnew, "An Experimental Study
of Transient Behaviour of anAdiabatic Continuous-Flow Stirred Tank
Reactor," Ind. Eng. Chem. Process Des. Dev., 14(2), 146 (1975)
7. Schmitz, R.A., R.R. Bautz, WH. Ray, and A. Uppal, "The Dynamic
Behavior of a CSTR: Some Comparisons of Theory and Experiment,"
AIChE Journal, 25(2), 289 (1979)
8. Grau, M.D., J.M. Nougues, and L. Puigjaner, "Batch and Semibatch
Reactor Performance for an Exothermic Reaction," Chem. Eng. and
Proc., 39(2), 141 (2000)
9. Ferrouillat, S., P Tochon, D. Della Valle, and H. Peerhossaini, "Open
Loop Thermal Control of Exothermal Chemical Reactions in Multi-
functional Heat Exchangers," Int. J. of Heat and Mass Transfer, 49
(15-16), 2479 (2006)
10. Bhanot, S., S. Kumar, and M.K. Vasantha, "On-line Instantaneous
and Delayed Control of a Fast Tubular Reactor," TENCON '98. 1998
IEEE Region 10 International Conference on Global Connectivity in
Energy, Computer, Communication and Control, Vol. 2, New Delhi,
India, 448 (1998)
11. Palanki, S., and J. Vemuri, "End-Point Optimization of Batch Chemi-
cal Processes," Proceedings 42nd IEEE Conference on Decision and
Control, 5, 4759 (2003)
12. Kulkarni, A.J., S.V. Patil, and V.K. Jayaraman, k!..... Based Local
Learning: Application to Process Engineering Problems," Int. Journal

of Chem. Reactor Eng., 1, A23 (2003)
13. Prat, L., A. Devatine, P Cognet, M. Cabassud, C. Gourdon, S. Elgue,
and E Chopard, "Performance Evaluation of a Novel Concept Open
Plate Reactor Applied to Highly Exothermic Reactions," Chem. Eng.
& Tech., 28(9), 1028 (2005)
14. Ausikaitis, J., andA.J. Engel, "The Behavior of the Iron(III)-Catalyzed
Oxidation of Ethanol by Hydrogen Peroxide in a Fed-Batch Reactor,"
AIChE J., 20, 256 (1974)
15. Ruiz, D., "Fault Diagnosis in Chemical Plants Integrated to the Infor-
mation System," Thesis, Universitat Politecnica de Catalunya (2001)
16. Shukla, PK., and S. Pushpavanam, "Parametric Sensitivity, Runaway,
and Safety in Batch Reactors - Experiments and Models," Ind. & Eng.
( I ..... . , Research, 33(12), 3202 (1994)
17. Gao, Q.Y., Y.L. An, and J.C. Wang, "A Transition from Propagating
Fronts to Target Patterns in the Hydrogen Peroxide-Sulfite-Thiosulfate
Medium," Phys. Chem., 6(23), 5389 (2004)
18. Karer G., G. Musi, and B. Zupan, "Predictive Control of Temperature
in a Batch Reactor with Discrete Inputs," Proceedings of the 13th
Mediterranean Conference on Control and Automation Limassol,
Cyprus, 855, June 27-29 (2005)
19. Solodov A., and V. Ochkov, "Stiff Differential Equations" in Differen-
tial Models. An Introduction with Mathcad, Springer Berlin Heidelberg
21. Perry R.H., and C.H. Chilton (Eds.), Chemical Engineers'Handbook,
5th Ed., McGraw-Hill Book Company, Chapt. 10 (1973)
22. Sutton, C.M., "In Situ Measurement of Resistance Thermometer Self-
Heating and Response Time," Meas. Sci. Technol., 5, 896 (1994)
23. Dally, J.W., WE Riley, and K.G. McConnell (Eds.), "Instrumentation
for Engineering Measurements," 2nd Ed., Iowa State Univ. (1993),
11.6 Dynamic Response of Temperature Sensors.
24. 1

Chemical Engineering Education

--- - ^ K.___________________________-



In a Graduate Course on Numerical Methods

Oklahoma State University * Stillwater, OK 74078
For the past 20 years, educators in the engineering field
have implemented ways of better engaging students,
including active and cooperative learning, learning
communities, service learning, cooperative education, inquiry
and problem-based learning (PBL), and team projects.[1 Peda-
gogies of these types have been shown to clearly be effective
in allowing students to improve their learning. It is no longer
important to consider simply the content covered in a course,
but how the material is going to be taught.
Researchers at McMaster University[21 found that to solve
problems successfully students need both comprehension
of chemical engineering and what are called general prob-
lem-solving skills. The researchers developed workshops to
explicitly develop skills deemed appropriate for improving
students' problem-solving skill and confidence.
Prince,31] in his article on active learning, provides additional
motivation for incorporating problem-based learning in the
classroom. One benefit that he shows is that PBL typically
involves significant amounts of self-directed learning on the
part of the students. The professor is no longer solely respon-
sible for delivering course content. Prince provides extensive
and credible evidence suggesting that faculty should consider
nontraditional models for promoting academic achievement
and positive student attitudes.

Based on the above evidence of the value of PBL and active
problem solving, the instructors reconstructed their graduate
course in Chemical Engineering Modeling. In Fall 2006, the
instructors introduced components to the course exposing

Eric L. Maase is a visiting assistant profes-
sor at Lafayette College and a research
associate in the Laboratory of Molecular
Bioengineering at Oklahoma State Univer-
sity. He received his B.S. degree from the
University of Maryland, his M.S. degree from
the Colorado School of Mines, and his Ph.D.
from Oklahoma State University. His research
interests include biochemical/biologicalmod-
eling and simulation, polymer science and
tissue engineering, thermodynamic modeling
of electrolytes, computational modeling and
design, and engineering education.
Karen A. High earned her B.S. and M.S.
degrees from the University of Michigan and
her Ph.D. from Pennsylvania State University.
She is an associate professor in the School
of Chemical Engineering at Oklahoma State
University. Her research interests are Sus-
tainable Process Design, Industrial Catalysis,
and Multicriteria Decision Making. She is
a trainer for Project Lead the Way pre-En-
gineering curriculum. She is also involved
with the development of an entrepreneurship
minor at Oklahoma State University.

� Copyright ChE Division of ASEE 2008

Vol. 42, No. 1, Winter 2008

students to case studies, active problem solving, teamwork,
and experimentation with an the aim of promoting creative
and critical thinking as students identify and practice the art
of modeling.

The Chemical Engineering Modeling course has been
structured to give the students experiences to meet the fol-
lowing objectives. At the end of the semester, the students
should be able to:
(a) Develop mathematical models describing chemical
engineering phenomena.
(b) Evaluate the assumptions, limitations, and restrictions
necessary to solve practical problems by mathematics.
(c) Use classical numerical techniques to solve the equa-
tions that resultfrom model formulation (ordinary and
partial differential equations, and linear and nonlinear
simultaneous algebraic equations).
(d) Become familiar with available computational tools
that incorporate these numerical techniques (specifi-
cally in MATLAB and VBA/Excel).

The course topics include:
(a) Classification of Equations
(b) Matrices
(c) Model Classifications/Concepts
(d) Nondimensionalization
(e) Taylor's Series

(f) Ordinary Differential Equations
(g) Linearization
(h) Linear/Nonlinear Equations
(i) Finite Differences

(j) Initial Value Problems (IVP-ODEs)
(k) Boundary Value Problems (BVP-ODEs)
(1) Partial Differential Equations (PDEs)

(m) Applications of .1. i,,,, .' " i, . I / Analysis
All students in the class are graduate status and, more im-
portantly, are also in their first semester of graduate school.
The course is a requirement for all chemical engineering
graduate students at Oklahoma State University. Students are
concurrently enrolled in graduate-level thermodynamics and
diffusional operations courses during this semester as well.
The diffusional course covers the equations of continuity and
focuses on analytical solutions to equations.
Activities and instructional materials teaching program-
ming and problem-solving principles using Visual Basic for
Applications (VBA) linked through Microsoft Excel and
MATLAB have developed during previous semesters as

course materials. These computer tools are employed by the
students in solving algebraic, ordinary differential equations
(initial value and boundary value) and partial differential
equations (parabolic and elliptic).
The instructors predominantly use two references during the
semester that help students develop and set up mathematical
models of physical phenomena. The Smith, Pike, and Murrill
(SPM) 41 methodology employs an "IOGA" (input, output,
generation, and accumulation) approach of building models
using a control volume. The Himmelblau and Bischoff (HB))51
methodology relies on beginning with equations of continuity
and "crossing" off terms that do not apply for any specific
situation. The course includes homework assignments in
mathematical theory (paper solutions), computer exercises
in numerical methods-both programming with VBA and
MATLAB -and exams evaluating student's understanding.
In past years students have worked on projects that have
them evaluate the effectiveness of the ODE boundary value
problem (BVP) strategies of matrix method, shooting method,
and the successive overrelaxation method (see Hornbeck&61
and Ri._' ~ I. The students are given the opportunity to pick
chemical phenomena that give rise to ODE BVPs. Prior
students have stated that they learned a great deal from such
projects but also have wondered if there were anyway to be
given a project that allowed them to model actual experimental

In modification to the course the instructors introduced
a semester-long group project where the students perform
experimental work in a lab, report and discuss their efforts
via memoranda to the instructors, and apply mathematical
theory, numerical methods, and developing modeling skills
to reach a specific goal.
The students were presented with the following problem
statement at the beginning of the semester:
"Design an encapsulated drug (ofyour group's choice) that
effectively delivers an appropriate dose for an appropriate
number of hours. Select the most -..-, i . f.. in.' method."
This problem statement is intentionally both inadequately
defined and open-ended. Defining the problem and the over-
all goals of the group projects is an aspect of the students'
Other educational objectives inherent in the problem state-
ment include:
* Consider different types/scales of phenomena (micro-
scopic, multiple gradient, maximum gradient, macro-
* Understand what physical models and mathematical
models are.
* Consider unsteady phenomena.

Chemical Engineering Education

* Develop models that cover as many types of phenom-
ena as are appropriate.
* Develop the maximum gradient modelfor spherical
* Develop experimental protocol around ..... i,,,i , drug
delivery with Tootsie Pops.
* Consider numerical techniques to solve the problem
(ODE vs. PDE).

The class meets twice weekly (Tuesdays and Thursdays)
and was structured so that Tuesday classes during the first
third of the semester are in computer labs and students are
provided instruction in VBA/Excel and MATLAB. Thursday
lectures provided instruction on mathematical principles and
modeling methodology. Over the remainder of the semester,
Tuesday lecture periods served as times for students to be
engaged either in group work or with instructors providing
support and mentoring in sessions where specific content
could be presented (lab safety discussion, review sessions,

Figure 1. Drug (lollipop) dissolution.

Vol. 42, No. 1, Winter 2008

presentation skills, etc). The students are initially daunted
by the possibilities inherent in the project but also are will-
ing-and in some cases, quite eager-to rise to the challenge
that the project represents. Two individual case studies, which
are introduced at separate times, provide the basic impetus for
the students' project development and direction. Discussions
of projects during class periods and project status reporting
(memos) provide additional mentoring moments for the in-
structors over the course of the semester.
To accommodate the extra time required in class and to
minimize lecturing, a variety of reading protests (administered
as online quizzes in Blackboard) were given to ensure stu-
dents comprehended material prior to class. Also, it allowed
the instructors to shorten "lecture" time and increase active
learning experiences.
The overall approach taken is intended to encourage critical
thinking and creativity in the process of problem solving by con-
necting applications of the mathematical and numerical methods
taught in lectures more closely to other engineering research
methods. The restructured course also mimics experiences that
students should expect to encounter while working towards a
thesis (M.S.) or dissertation (Ph.D.) in a graduate program.
Homework assignments and exams were administered to
students to ensure that they were able to individually learn
and demonstrate knowledge in the content. The overall course
grade was from exams (60%), individual homework (20%)
and the group project (20%). This ensured that the students
were learning all of the course objectives.

At the start of the semester, the students were put in groups.
The questions in Table 1 were used to ensure heterogeneity in
student groups based on experience, gender (when determin-
able from name), and ethnicity.

Once the students were in groups, the instructors led an ex-
perimental session where the students began their project work
by examining a dissolution case study. A laboratory was set up
that was designated for the students use (see Figure 1).

Questions Used to Group Students
Encapsulated Drug Project, Aug. 31, 2006
On the card, write:
Your name
Your home (i.e., China, Michigan, Antarctica)
Your degree (i.e., M.S. Chemical Engineering, Ph.D. Chemistry)
Your computer experience (High/Medium/Low)
Your experimentation experience (High/Medium/Low)
Your transport equation experience (High/Medium/Low)
Your numerical method experience (High/Medium/Low)

The experiments are driven by the goals and basis outlined
in an initial case study-providing a starting point for the
student projects. The intent of the initial lab period is also to
provide instruction on lab notebooks and experimental con-
cerns such as safety and experimental measurement, as well
as to outline the challenges that experimental work entails.
The class would later discuss and determine the safety rules
for the lab with all students signing the resulting safety agree-
ment as shown in Table 2. A discussion of this nature resulted
in a more open-ended perspective on safety along with a dif-
ferent rational for and understanding of other users' feelings
in a laboratory environment and their sense of safety.


Case studies served as the primary creative impetus
for students' projects. While the students may be eager
to undertake the challenge posed by the project they also
need a sense of the structure that allows such studies
to be successful. The initial case study provides their
starting point.
Case Study #1 -Experimental Investigation
and Modeling
This case study, from Jensen,[81 is intended for freshman
students. The goal is for students to investigate the factors
affecting the dissolution of a proposed coating for a new
pharmaceutical. The story is written from the perspective
of two chemical engineering students in "ChemE 101," a
fictional course. The learning objectives for the students
reading and working on the case study are: 1) describe
the factors affecting mass transfer between phases; 2)
collect and evaluate data on mass transfer; and 3) evalu-
ate a model to describe mass transfer.
These objectives remain appropriate for the graduate
students' projects as well. The intent of the case study is
simply to first determine the appropriateness of the linear
model for a lollipop (the students used Tootsie Pops) in
order to experimentally investigate a system having an
encapsulating material over an inner "drug" material. The
students were then tasked with developing and outlining
their own project: developing a better model of encapsu-
lated drug delivery, as can be seen in Table 3.
As students were completing the initial experimental
aspects of the project, they were allowed to ask for
materials to be provided for them to help solve various
problems. This exposed the students to a just-in-time
teaching[91 or lecture-by-demand strategy.
Once students' initial experimental and numerical
techniques have been applied to Case Study #1 and they
have begun further experiments of their own designs,
a second case study is introduced. The purpose of the
second case study is to introduce the concepts of com-
putational modeling and refocus students on the overall


aims of encapsulated drug delivery. This case study is dis-
cussed after the necessary computer instruction (VBA/Excel
and MATLAB) is complete.
Case Study #2-Modeling and Numerical
The first step in any modeling process is to simplify the
system of interest into essential aspects that can be math-
ematically described while still retaining essential features
of the larger real-world system of interest. This case study
provides an example of a model of a human body interacting
with drugs delivered orally.
The second case study arises from an assignment previ-
ously developed and used in an undergraduate Computer

Laboratory Safety/Practices for 307 EN CHE 5743 Class Lab
1. Safety glasses (eye protection) are encouraged
2. Move deliberately-no running
3. No more than 2 groups and 1 individual in 307 at one time
4. No food/drink/gum
5. Closed shoes and full pants (or skirt) will be worn all times in 307
6. Keep lab clean and equipment stored in a neat fashion
7. Unplug all equipment when done
8. Close door and lock when the last person leaves
Name (print)

First Project Memo
To: Students of Chemical Engineering 5743
From: The Instructor
Date: Aug. 31, 2006
Subject: 1st group memo for encapsulated drug class project.
A 1-2 page memo is due Sept. 7 at the beginning of class that answers/consid-
ers the following.
1) Initial experimental protocol for your physical model
a. Equipment/supplies needed-your group will be
responsible for procurement
b. Data collection procedure
c. Other information
2) Initial computer modeling/numerical methods for mathematical model
a. Computer software to be used
b. Numerical methods to be used
c. Other information
3) Simplifications made to start the initial work
4) Future considerations that more appropriately model the real system
5) Questions for the instructor.
6) A reference list to help in the project. Include at least 2 references per
group member.

Chemical Engineering Education

Programming for Engineers course taught at Oklahoma State
University by the co-author. The case study outlines a model
of the human body as a set of three interacting compartments
where each represents a broader aspect of human physiology
deemed important in drug delivery. The three interacting com-
partments are referred to as the oral, central, and peripheral
The oral compartment represents primarily the mouth and
stomach where a drug is gradually dissolving. The central
compartment represents both the blood stream (circulatory
system) and those organs that are highly perfused (liver, kid-
ney, lungs, etc.). The peripheral compartment includes muscle
and tissues where diffusion of medications is often slower and
where drugs may build up in concentration over time.
Case Study #2 also defines the specific interactions among
the compartments with governing equations for the transfer
of a drug to and from compartments, as shown in Figure 2.
These equations include rates/pathways for a drug leaving or
entering the central compartment, where k is the rate drugs are
eliminated from the body, k, is the drug metabolizing, and k
is the drug being further absorbed (to the peripheral compart-
ment). All of the rates are dependent on the concentration of
the drug currently in the central compartment, Co. The central
compartment balance also includes terms describing the drug
entering from the oral compartment, kCo, and returning from
the peripheral compartment, kpCp.
The case study tasks students with creating a computer
program capable of solving the drug delivery model, inves-
tigating simulated behaviors, and discussing the results of
their efforts.
The overall perspectives and intentions of Case Study #2
prove valuable for the graduate course as well. In extending
the basic princi-
ples outlined by
Oral the second case
dC study, graduate
d =-kAo students will ex-
I, & plicitly consider
multiple drug
L doses, statisti-
' cal uncertainties
Central inherent in hu-
cc /man reactions
cdk=- E +, )C+C to medications,
the limitations
Sand assumptions
necessary for

c, A4 Figure 2.
= -kpCp+kcC Case Study #2:
d1 Drug delivery

Vol. 42, No. 1, Winter 2008

complex system modeling, and applications of numerical
methods discussed in the course as tools for solving systems
of equations.

Case studies served as the primary impetus for students'
projects in the course. Through use of the first case study stu-
dents were provided an initial direction from which they could
develop their own experimental goals and were presented
with the realities inherent in both mathematical modeling
and experimental work.
The second case study provided a different perspective for
the students where the challenge became more open ended
and constructing better models-in tractable mathematical
terms able to capture essential behaviors of the actual system
of interest-was a goal. These models and related simulations
provided additional paths for examination of complex systems
where important factors and influences can be estimated and
examined computationally when experimental data or meth-
ods are difficult at best.
The modeling and numerical methods in computer
simulation complemented the experimental laboratory
work initiated with Case Study #1. The second case also
served to remind the students of the overall project goals)
as outlined in the original problem statement. Translat-
ing mathematical and conceptual models into computer
simulations also tasks the students with understanding and
integrating their previous experimental approaches (those
completed as well as those under way) into a more com-
plex system description. This aspect of applied numerical
investigations asks students to determine a better overall
"picture" and to hopefully find a more effective resolution
for the project goalss.
The case studies used in the course are complementary
tools in engineering research. The implantation of the cases
in the course provided introductions for students to the role
of laboratory experiments along with computational numeri-
cal simulations.

Over the course of the semester, project status, learning,
and student assessment (grading) are addressed by asking the
student groups to prepare memoranda at appropriate points
during the project. Table 4 (next page) lists the additional proj-
ect status memos that student teams were asked to prepare.

The written projects were due the second to last day of the
semester. Presentations were on the last day of the semester.
Four days before the final presentations were given, the in-
structors provided the students with instruction on creating
appropriate Microsoft PowerPoint documents and delivering
an effective presentation.


Student Project Examples

One group undertook an investigation of the design of a
cost effective controlled-release drug delivery system for
Digoxin.1101 In their efforts they proposed a design - a four-
layer delayed-release pill-intended to sustain optimal blood-
plasma concentration over a 24-hour period; the proposed

Project Status Memos
2nd memo due September 25th at the beginning of class.
1) Revise your 1st memo considering that your first experiment
should take 1 hour and 15 minutes (class time). You should keep
this very simple and consider using the physical and mathemati-
cal model in case study #1. You need a very detailed plan for
your experiment. Where will you get the equipment for the
2) Develop the mathematical model in the case study using the:
a. Himmelblau and Bischoff approach
b. Smith, Pike, and Murrill approach
3rd memo due October 6th (4:30 p.m. 423 EN)
1) Develop the mathematical model in the case study using the:
a. Himmelblau and Bischoff approach
b. Smith, Pike, and Murrill approach
Follow Example A approach.
2) Discuss your data and numerical model from your experi-
ments on Sept. 28. Use the in-class discussions from Oct. 3 as
your guide.
4th memo due October 20th (4:30 p.m. 423 EN)
1) Find your initial experimental data from Sept 28. Plot the
data r vs. t. Does equation 20.1 of the case study "fit" your data?
2) Rerun your experiment with a solid lollipop. Plot the data r
vs. t. Does equation 20.1 "fit" your data?
3) Develop a next-level model (mathematical and physical)
that gets your group toward the ultimate goal of designing an
encapsulated drug. Clearly plot your data from your experiment
and show how well it "fits" your data from your experiment
(physical model). Discuss the appropriateness of the models
(physical and mathematical). Don't get too fancy yet!
5th memo due Friday December 1 at noon in 423 EN
1) Tell us how your group is considering determining "cost
2) What are your aims (i.e., which drug are you using, how are
you delivering the drug, where are you delivering it to)?
3) How are you focusing the problem statement into specific
4) Progress on physical model and experiments.
5) Progress on mathematical model.
6) Back up your claims with references to literature. What are
you using for your main references? Include an accurate and
complete bibliography that you have been compiling since the
start of the project.

four-layer coating dissolution is shown in Table 5.

The students generated a conceptual model, determined
simplifications and assumptions, and created numerical simu-
lations of their system. An example is shown in Figure 3 (the
graph is unedited student work). The work shown represents
only one example of the students' thinking and submitted
work during the course.

In their final report, the economics of the proposed drug
delivery method were outlined along with their current in-
vestigations to date.

Overall the work submitted during the course, along with
the interactions the instructors had with the students, showed
students to be very actively engaged with course materials
and found to be meeting or in many instances exceeding the
instructors' goals for the reconstructed course.

Survey Assessment
An end-of-course assessment was given to the 13 graduate
students in the course. University Institutional Review Board
approval was given for the survey. The students gave consent
for the authors to use the results of the survey in presentations
and publications. The survey was anonymous and examined
after the conclusion of the semester. The survey had three

Section 1

Q1. I can appropriately use a lab book.

Q2. I am confident with my experimental skills.

Q3. I am comfortable working in groups/on teams.

Q4. I am able to give effective presentations using Power-

Q5. I feel confident that I can use VBA/Excel to solve
problems in this course.

Q6. My numerical methods skills are appropriate to solve
problems in this course.

Q7. I feel confident that I can use MATLAB to solve
problems in this course.

Q8. I am able to develop mathematical models that ap-
propriately consider relevant phenomena.

Q9. I understand the need to follow appropriate lab safety

Student Designs: Four-Layer Pill
Layer Coating (hrs dissolve) Digoxin (pg)
1 .167 2.35
2 5.83 0.9
3 4.67 0.9
4 4.67 0.9

Chemical Engineering Education

Q10. I enjoy a problem-based learning environment.

Q11. I can search the literature for relevant background
material for projects.

Q12. I understand the components of a research report /
Q13. I know what goes into a developing a good presenta-

The students were given a 5-point Likert scale to evaluate
their agreement with the above statements. Table 6 shows the
results for Section 1 of the assessment.

All scores listed were above neutrality, indicating a positive
agreement of all questions. The students felt the most agree-
ment with the statement "I enjoy a problem-based learning
environment." In the authors' opinion the students gained a
lot from the active learning experience and were engaged
throughout the entire project and course. The next highest
agreement was for the statement "I know what goes into
developing a good presentation." The instructors felt that this
is principally due to the class period four days prior to class
end, in which presentation skills were discussed.

The students had the lowest level of agreement for the
statement "I feel confident that I can use MATLAB to solve
problems in this course." This is not surprising since the
co-authors spent a majority of the
class discussing solution strategies
in VBA/Excel. The next lowest level 2.5
of agreement was for the statement
"I am able to give effective pre-
sentations using PowerPoint." This 2 - ""
was interesting given their strong
agreement to knowing what goes
into a good presentation. One pos- 1 5
sibility is that this result may be due E
to language challenges that many of
foreign students experienced. The i 1
students also had not given their pre-
sentation yet. The project was turned
in on the same day as the assessment o 5

D Figure 3. Digoxin concentra- _
tion in blood plasma. o

and the presentation was the next day. Foreign students prob-
ably had anxiety about giving a presentation based on a project
they just turned in and delivering a presentation the next day
in a second language.

Section 2
Q1. How would you rate your knowledge of MATLAB?

Q2. Did you use MATLAB to solve problems in the
course? Why or why not?

Q3. How would you rate your knowledge of VBA/Excel?

Q4. Did you use VBA/Excel to solve problems in the
course? Why or why not?

Q5. How would you rate your math skills, particularly
eigenvalues and matrix calculations?

Q6. How would you rate your knowledge of the transport
equations (the equations of continuity for energy,
mass, momentum)?

Q7. How would you rate your knowledge of numerical
methods (finite difference.)?

Q8. How helpful were the MATLAB help sessions that
Dr. High provided?

Q9. How helpful were the VBA help sessions that Dr.
Maase provided?

5 = Strongly agree, 4 = agree, 3 = neutral, 2 = disagree, 1 = strongly disagree.
*Indicates where there is a tie in the ranking.

Vol. 42, No. 1, Winter 2008

20 25

Time (hrs)

Q10. How helpful was Blackboard to you?

All of these questions had room for student comments.
Table 7 shows assessment results.
The highest rating for the five questions where students rat-
ed themselves was in their knowledge of numerical methods
(Q7). They also felt confident with their math skills (Q5). The
lowest was for their knowledge of MATLAB (Q1), again this
wasn't surprising considering that VBA was emphasized.
For the questions where the students responded to how
helpful various resources were, they gave highest score for the
VBA sessions and for the Blackboard classroom management
system and lower for the MATLAB sessions.
For Q2 (did you use MATLAB?), eight students said yes,
three said no and two didn't answer. Student comments for
why they used MATLAB included:
"Only when we were assigned to"

"Sometime because it is easier to use"

"Its easy to use"

And for comments why they did not use MATLAB:
"Used VBA (more comfortable with)"

"I don't have enough experience to work on it com-

For Q4 (did you use VBA?), 12 said yes and one didn't an-
swer. Student comments for why they used VBA included:
"I love to use it."

"For numerical methods and solve diff. equations"

"Because I was trying to learn VBA."

"It's very simple."

"Easy to use and easily accessible."

"(For) all problems because of ease of interface and
prior background."

Section 3

This section was where students could write free-form
comments to the questions. Student comments are added after
each statement. See the Appendix for comments.
According to the assessment, it appears that the students
appropriately developed their
expertise in four areas: math-
ematical modeling; experimental
strategies; numerical methods; Q1
and computer strategies on MAT- Value 2.46
LAB and VBA/Excel. It is also Rank 5

clear that they enjoyed the prob-
lem-based environment in which
the course was administered.

For Q1, Q3, Q5, Q6, Q7: 4=1
For Q8, Q9, Q10: 4=Extrem
*Indicates where there is a ti

This project proved invaluable for beginning graduate
students to identify the connection between modeling, math-
ematical descriptions, and experimental reality. It was also
intriguing for the students to not only consider mathematical
but physical models in which experimentation with the real
system is, at best, impractical-as is the case in biological
systems. The students gained valuable teamwork skills not
common for many graduate students as well as improvement
in communication skills. The assessments show that the stu-
dents were particularly engaged in the course.

The approach taken encourages critical thinking and creativ-
ity in the process of problem solving, while connecting appli-
cations of mathematical and numerical methods more closely
with other engineering research methods. The restructured
course also mimics experiences that students can expect to
encounter while working toward a thesis (M.S.) or dissertation
(Ph.D.) in a chemical engineering graduate program.


1. Smith, K., S. Sheppard, D. Johnson, and R. Johnson, "Pedagogies
of Engagement: Classroom-Based Practices," J. Eng. Ed., 94(1), 87
2. Woods, D., A. Hrymak, R. Marshall, P Wood, C. Crowe, T. Hoffman,
J. Wright, P Taylor, K. Woodhouse, and C. Bouchard, "Developing
Problem Solving Skills: The McMaster Problem Solving Program,"
J. Eng. Ed., 86(2), 75 (1997)
3. Prince, M., "Does Active Learning Work?A Review of the Research,"
J. Eng. Ed., 93(3), 223 (2004)
4. Smith, C., R. Pike, and P Murrill, Formulation and Optimization of
Mathematical Models, International Textbook Company, Scranton, PA
5. Himmelblau, D., and K. Bischoff, Process Analysis and Simulation,
Wiley, NY (1968)
6. Hornbeck, R., Numerical Methods, Prentice Hall, Upper Saddle River,
NJ (1982)
7. Riggs, J, An Introduction to Numerical Methods for Chemical Engi-
neers, Texas Tech University Press, Lubbock, TX (1994)
8. Jensen, J., "Chapter 20- Dissolution Case Study," in A Users Guide to
Engineering, Prentice Hall Source, Upper Saddle River, NJ (2004)
9. Novak, G., A. Gavrin, W. Christian, and E. Patteson, Just-In-Time
Teaching: Blending Active Learning with Web Technology, Benjamin
Cummings, San Francisco (1999)
10. Weiss, M., and W Kang, "Inotropic Effect of Digoxin in Humans:
Mechanistic Pharmacokinetics Pharmacodynamic Model Based on Slow
Receptor Binding," Pharmaceutical Research, 21, 235 (February 2004)

Assessment Results for Section Two
Q3 Q5 Q6 Q7 Q8 Q9 Q10
2.85 3.00 2.77 3.08 2.77 3.08 3.08
3 2 4 1* 3 1* 1*
Extensive, 3=Pretty good, 2=Limited, 1= None
ely, 3=Very, 2=Somewhat, 1= Not helpful
e in the ranking.

Chemical Engineering Education


The following are free-form responses to the survey ques-
tions (comments are verbatim).
Q1. Describe how well you thought the project allowed you
to understand the use of physical models and experiments to
simulate phenomena (particularly biological).

"It made me understand physical and experimental drug

"I never knew how to model physical phenomenon before I
enrolled in this course. It helped me improve my computer

"Oh great!!! Actually being new to grad studies, it was a
welcome breeze of air."

"The project was useful in .... . , ,,,, the concepts of
physical models and experiments."

"The project gave me some insight how to simulate a physi-
cal model and to develop mathematical model."

"I think this project is very helpful to understand how to
construct the physical models, carry out experiments, con-
struct and solve mathematical models."

"It helped me a lot in developing a good understanding
about the phenomena and in visualizing the problem."

"I believe now I understand the concept of having a physi-
cal model and fiii ,, mathematical model."

"Got a lot of knowledge about the mathematical model."

"Very well. P 1, .;,,, pure experimental work to comprehend
modeling was good."

"The project idea was good. Saw relationship."

Q2. Describe how well you thought the project allowed
you to understand the use of mathematical models to simulate
phenomena (particularly biological).

"It was important as we can now understand how important
it is to make simple assumptions. They make almost-real
models with less effort."

"The late equations could be easily coded in VBA, and the
experimental results could be studied easily."

"Yeah, it helped a lot, but biological systems are complex,
so we had to make several assumptions to mathematically
represent it."

"I think this project helped me to understand how to convert
the physical phenomena into mathematical models and
know the mechanism of the phenomena to make suitable
assumptions and negligibility."

"It assisted me tremendously. Now I relate the methods
learned to real-life situations."

"This seemed a bit vague .1-i ....,i ..t. . we had some dif-
ficulty j i.ii;,, physical to mathematical."

Q3. Discuss how the project helped you to improve your
numerical method skills.

"I am not sure it helps me much on my numerical method

"Not much, I was good at them."

"Gave me a better understanding of Runge Kutta and other

"The project gave opportunity to solve numerically the

"By solving project, I used a lot of numerical methods, like
how to solve algebraic equations, differential equations,
and regress the parameters."

"The project helped me in improving my numerical skills."

"I can now understand numerical analysis and relate to
real-life situations."

"Really helped us to model real-time problems.

"I felt like the project was very numerical-method light."

"I don't know how. But!!! It made my numerical method
skill to improve!"

Q4. Discuss how the project helped you to improve your
computational skills.
"It helps me use Excel and VBA."

"Had homework, projects, and a good discussion in the lab."

"May have helped more had I been serious i, J ii,, day 1."

"Yes it helped me to improve as we faced the real life prob-
lems. "

"This program needs to use Excel and VBA to solve the
program. It's helpful."

"Yes, it did improve my computational skills."

"I felt like the project required very little computational

Q5. Describe your ability to develop a cost-effective en-
capsulated drug.

"Cost effectiveness is dependent on how many people need
this and how I can optimize."

"At this point, on the basis of literature review, I can de-
termine the factors which can/might affect the cost of drug,
and they can be varied according to the requirement."

"I can develop a cost-effective encapsulated drug consider-
ing mathematical models, economic evaluation models."

"I couldn't actually develop an exact cost-effective project."

Vol. 42, No. 1, Winter 2008

"No knowledge whatsoever. I wish you had given back-
ground info about drug analysis or economics."

"I search the related information and get a rough knowl-
edge about the cost-effective encapsulated drug."

"The cost-effective part was difficult to evaluate."

Q6. Discuss the usefulness of the initial case study and Dr.
Maase's case study to the final project.

"The initial case study helped in understanding the aim of
the project."

"Dr. Maase's case study was very helpful in , ,111, the
code plus understanding the model."

"Very helpful. Great instructors and great interactive

"Very useful! Especially Dr. Maase's case study which gave
me a better perspective of the model.

"It was helpful as it gave insight how to simplify the real
system to develop mathematical models."

"They are all very useful to understand this project."

"Initial case study and Dr. Maase's case study both give me
an introduction, details, and hints to the final project."

"Helpful guide to know how to start."

Q7. Other comments.
"Give us all a good grade finally, please." 7

"I felt like the project was at times taken to be more dif-
ficult than it really was. Many of the memo requests seemed
identical." a

Chemical Engineering Education

Random Thoughts...



North Carolina State University
Education Designs, Inc.

Myths about student ratings of teaching abound on
every campus, usually accompanied by confident
assurances that ratings are just popularity contests
that reward entertainers and penalize the best teachers. (Re-
markably, the second group always seems to include the ones
doing the assuring.) Some years ago we surveyed the myths
and summarized the extensive research that showed most of
them to be wrong."1 Now it's 15 years later and a lot more
research has been done, with similar outcomes. Unfortunately
the myths are still alive and well, so here is the 2008 version
of what "everyone knows" about student ratings and how
much of that wisdom is supported by research.

Myth. Student ratings are not valid assessments of teaching
Fact. False. Thousands of research studies have shown that
student ratings correlate positively with every other measure
of teaching effectiveness, including alumni ratings, peer rat-
ings, administrator ratings, measures of learning (e.g., stan-
dardized tests, common exams in multi-section courses, and
ratings of student portfolios), and student motivation. 21 The
magnitude of the observed correlations varies considerably
across individual studies and a few studies report contradic-
tory results, but the weight of the evidence is clear. If students
consistently say someone's teaching is good or bad, they're
almost certainly right.
Myth. The highest i. i, ii % go to the easiest courses.
Fact. False. Up to a point, courses rated as more difficult on
average get higher ratings than easier courses, with ratings
only beginning to fall when courses reach levels of difficulty
beyond the backgrounds of most enrolled students. In a recent
study of 1,045 engineering, science, and humanities courses
at two universities, Dee 31 found that student perceptions of

� Copyright ChE Division of ASEE 2008
Vol. 42, No. 1, Winter 2008

course workload were not significantly different for courses
in the top and bottom quartiles of student ratings, with the
marginally higher workload rating going to the courses in
the top quartile.
Myth. Bad teachers who are easy graders get higher evalu-
ations than good teachers who are strict graders.
Fact. False. Individual instructors who give high grades rela-
tive to local averages may get higher ratings than they would
if their grades were lh i, .,'"1 but no studies have turned up
ineffective teachers who got high ratings just by giving high
grades.E5] However, the possibility that it could happen sup-
ports the common recommendation to use multiple sources
of assessment data.
Myth. They may not like me now because I'm 1 Ii ... -t n and
maintain high standards, but in afew years they'll appreciate
how good a teacher I was.
Fact. Generally false-it happens sometimes, but not often.1 2]
Alumni ratings correlate significantly with student ratings
given previously to the same instructor. If your students
think you're a great teacher now, most will still remember
you fondly in the future, and if they think you're lousy, don't

Richard M. Felder is Hoechst Celanese
Professor Emeritus of Chemical Engineering
at North Carolina State University. He is co-
author of Elementary Principles of Chemical
Processes (Wiley, 2005) and numerous
articles on chemical process engineering r
and engineering and science education,
and regularly presents workshops on ef-
fective college teaching at campuses and
conferences around the world. Many of his
publications can be seen at edu/felder-public>.
Rebecca Brent is an education consultant
specializing in faculty development for ef-
f fective university teaching, classroom and
computer-based simulations in teacher
education, and K-12 staff development in
language arts and classroom management.
She codirects the ASEE National Effective
Teaching Institute and has published articles
on a variety of topics including writing in un-
Sdergraduate courses, cooperative learning,
Public school reform, and effective university

expect to start getting holiday greetings from them in five or
10 years.

Not all common beliefs about evaluations are wrong, of
course. It's true that teachers who are enthusiastic and car-
ing tend to get better ratings than those who are reserved and
distant, but so " hal ' Enthusiasm and caring of instructors
also correlate with motivation and learning of their students, 21
suggesting that the higher ratings are probably deserved. It's
also true that other things being equal, elective courses tend to
get higher ratings than required courses, upper-level courses
tend to get higher ratings than lower-level courses, small and
moderately sized classes get higher ratings than very large
classes, student ratings in engineering and the sciences are
lower than ratings in other fields, and female instructors in
engineering and the physical sciences get lower ratings than
male instructors.[2] On average these effects are small, but
they exist and should be taken into account when ratings are
used to make decisions about such things as reappointment,
tenure, promotion, and merit raises.
In short, student evaluations have high levels of reliability
and validity and should always be part of the process used
to evaluate teaching. There are some aspects of a course that
students are in no position to evaluate, however, including
whether the course learning objectives are appropriate, the
content is current with the state of the field, and the course
adequately prepares the students for subsequent courses
in the curriculum. Those things can only be evaluated by
knowledgeable peers. Student ratings should therefore not
be the sole source of teaching assessment data but should be
supplemented with peer ratings and other measures of teach-
ing effectiveness.61] If different sources agree, as they usually
will, it's a good indication that the overall assessment is a fair
one; if they disagree, it's a red flag, and an effort should be
made to find out what's going on.
Since student ratings will undoubtedly remain central to
teaching evaluation (as they should), everything possible
should be done to make them as effective as possible. The
following recommendations -most of which are drawn from
the papers cited in Reference 2-address that goal.

* Use a rating form that has been developed with the as-
sistance of someone knowledgeable about educational
measurements. There is a science to survey construc-
tion in general and educational rating instrument con-
struction in particular. Either use a form that has been
developed and validated elsewhere, such as the IDEA
Student Ratings of Instruction system ( idea.ksu.edu/>) or TCETools (),
or work with an education specialist on your campus
or an external consultant.

* Don't trust i.,1,, collected from fewer than 10
students or less than 2/3 of a class, and don't make
personnel decisions based on 1.,11, from a single

* Use a few global or summary items with Likert-scale
(1-5) i ,r1,,, for summative evaluation (evaluation used
to help inform personnel decisions), and a longer list of
more specific items for formative evaluation (diagnostic
evaluation used to help instructors improve their teach-
ing). Global items correlate more strongly with student
learning than more specific items do, and you'll get a
higher rate of return if there are fewer questions.

* When i. i, i1,, , .i;,u, i, . remember that they may be
i, 7/ti affected byfactors other than the quality of
instruction, such as the nature and level of the course,
the class size, and the gender of the instructor. The
IDEA system and TCETools include provisions for
taking these factors into account.

* Try to persuade students that their j.,1,,, will be
considered carefully and may have an impact on fac-
ulty personnel decisions and decisions about teaching
assignments. If you can make this case convincingly,
most students will take the ratings seriously and you
should get a good rate of return. If you can't make the
case, there is no reason the students should take the
ratings seriously and you should not be surprised if
they don't.

1. Felder, R.M., "WhatDo They Know, Anyway?"( ,.... 3 /., 2( ,),
134-135 (1992), html>
2. (a) Cashin, WE., "Student Ratings of Teaching: Recommendations
for Use," IDEA Paper No. 22, Center for Faculty Evaluation and De-
velopment, Kansas State University, January 1990, ksu.edu/papers/Idea_Paper_22.pdf>; (b)Ory, J.C., i i..II. Thoughts
and Concerns About Student Ratings," New Directionsfor Teaching
and Learning, 87, 3-15 (Fall 2001); (c) McKeachie, W.J., "Student
Ratings: Their Validity of Use," American Psychologist, 52(11), 1218-
1225 (1997); (d) McKeachie, WJ., Teaching Tips, llth Ed., Boston,
Houghton Mifflin Co., 326-328 (2002)
3. Dee, K.C., "Student Perceptions of High Course Workloads Are Not
Associated With Poor Student Evaluations of Instructor Performance,"
J. Eng. Ed., 96(1), 69-78 (2007)
4. Greenwald, A.G., and G.M. Gillmore, "Grading Leniency is a Remov-
able Contaminant of Student Ratings," American Psychologist, 52(11),
1209-1217 (1997)
5. Abrami, PC., WJ. Dickens, R.P Perry, and L. Leventhal, "Do Teacher
Standards forAssigning Grades Affect Student Evaluations of Instruc-
tion?" J. Educational Psychology, 72, 107-118 (1980)
6. Felder, R.M., and R. Brent, "How to Evaluate Teaching," Chem. Eng.
Ed., 38(3), 200-202 (2004), Columns/Teacheval.pdf>. See also "A Protocol for Peer Review of
Teaching," ( Review).pdf>) and "If You've Got It, Flaunt It: Uses and Abuses of
Teaching Portfolios" ( Portfolios.html>) O

All of the Random Thoughts columns are now available on the World Wide Web at
http://www.ncsu.edu/effective_teaching and at http://che.ufl.edu/~cee/

Chemical Engineering Education


Introducing Stochastic Simulation of Chemical Reactions

Using the Gillespie Algorithm and MATLAB:


Kansas State University * Manhattan, KS 66506

I n their contribution, Martinez-Urreaga and his collabora- Andres Argoti is a research associate in the Department of Chemical
tors'11 have presented the stochastic simulation of a chemi- Engineering at Kansas State University. He received his B.S. in Chemical
cal reaction-specifically, a simple reversible chemical Engineering from Universidad Nacional de Colombia, Bogota, and his
M.S. and Ph.D. from Kansas State University, both in Chemical Engineer-
reaction whose rate law is of the first-order, or linear. Their ing. His major research interest is in the stochastic analysis, modeling,
work is indeed useful and informative: Randomly or stochas- and simulation of chemical and biochemical processes.
tically behaving processes or systems are ubiquitous in the L.T. Fan, a University Distinguished Professor, holds the Mark H. and
chemical or allied industries.[29] Margaret H. Hulings Chair in Engineering and is director of the Institute
of Systems Design and Optimization at Kansas State University. He
The current contribution involves two aspects and is served as department head of the Department of Chemical Engineering
intended to complement and augment the work of Martinez- between 1968 and 1998. Currently, his major research interests are in
stochastic analysis and modeling and in process-network synthesis
Urreaga and his collaborators[11 to enhance its usefulness, based on process graphs (P-graphs). Fan received his B.S. in chemical
First, the reacting system (i.e., the simple reversible chemi- engineering from National Taiwan University, M.S. in chemical engineer-
ing from Kansas State University, and M.S. in Mathematics and Ph.D. in
cal reaction) illustrated by them is modeled as a stochastic chemical engineering from West Virginia University.
process, specifically a birth-death process -the most general
prbcss, ifical a ir process he e stnt del ien Juan Cruz is a graduate research assistant in the Department of Chemical
subclass of Markovian processes.10 "] The resultant model, in Engineering at Kansas State University where he is currently carrying out
turn, gives rise to the master, or governing, equation11, 12] of his doctoral studies. Cruz received his B.S. in chemical engineering from
the birth-death process whose solution renders it possible to Universidad Nacionalde Colombia, Bogota; subsequently, he had a two-
year stint performing research at University of Puerto Rico, Mayagaez. His
obtain the analytical expressions for the process' means and major research interests include nonaqueous enzymology, biomolecular
higher moments about the means, e.g., the variances, which engineering, and molecular modeling and simulation.
are collectively a manifestation of the process' inherent fluc- Song-tien Chou obtained his B.S. in chemical engineering from National
tuations. Second, the master equation of the birth-death process Taiwan University, and earned M.S. degrees in both chemical engineering
is stochastically simulated via the Monte Carlo method. and statistics, and a Ph.D. in statistics, all from Kansas State University.
is stochastically simulated ia the Monte Carlo method1He is teaching in the Department of Finance and Banking at Kun Shan
The method is implemented by the time-driven approach[14 University in Taiwan; he was the department chair between 2000 and
5] in addition to the event-driven approach[10, 1523] adopted by 2002. His research interests include the application of stochastic pro-
cesses, risk analysis, and environmental engineering.
* Department of Finance and Banking, Kun Shan University,
Yung-Kang City, Tainan Hsien, 71003 Taiwan � Copyright ChE Division of ASEE 2008
Vol. 42, No. 1, Winter 2008 3.

Martinez-Urreaga and his collaborators11l on the basis of the
Gillespie algorithm.l0, 17 18] The event-driven approach entails
the determination of the probability distribution of a random
waiting time or period of quiescence19, 20]; the simulation clock
is advanced according to this random waiting time. Unlike
the event-driven approach, the time-driven approach advances
the simulation clock by a fixed time increment,114 15] which can
be estimated as a function of the intensities of transition
embedded in the master equation of the process. The time-
driven approach is especially useful when the probability
distribution for the random waiting time in the event-driven
approach is exceedingly complex to determine. The means
and variances of the random variables characterizing the
birth-death process have been computed in implementing the
event-driven and time-driven approaches. For validation, the
analytical solutions from the stochastic model as well as the
numerical results from the Monte Carlo simulation are com-
pared with each other. In addition, these results are validated
by comparing the values of the means obtained from the
solution of the deterministic model for the simple reversible
chemical reaction as presented by Martinez-Urreaga and his
The two approaches for stochastic simulation are further
illustrated with the photoelectrochemical disinfection of
bacteria whose rate of decaying is linear. The results from
simulation are validated by comparing them with the available
experimental data[241 as well as with the analytical solutions
derived from the corresponding stochastic model.

The system under consideration comprises the populations
of molecules of species A and B per unit volume that are be-
ing chemically transformed into one another according to the
simple reversible reaction,

A B (1)

Clearly, the population of A decreases and that of B increases
when one molecule of A transforms into one molecule of
B with forward reaction-rate constant k1. Similarly, the
population of B decreases and that of A increases when one
molecule of B transforms into one molecule of A with reverse
reaction-rate constant k2. The system constitutes, therefore, a
special class of Markov processes, the birth-death process.11"
12] Thus, the numbers of molecules of A and B at time t, NA(t)
and NB(t), respectively, can be taken as the random variables
of the birth-death process of interest. A realization of NA(t) is
denoted by nA, and that of NB(t) by ng. At any time t, the total
number of molecules in the system remains constant; hence,
the following relationship holds;
NA (0)+N(0)= NA (t)+N,(t)

nAo +nBo = NA (t)+ N (t)


where NA(O) or nA0 is the number of molecules of A at the
outset of the reaction, i.e., at t = 0. Similarly, NB(0) or ni0 is
the number of molecules of B at t = 0. These two quantities
are constant, thereby indicating that only one of the random
variables in Eq. (2) is independent. By denoting the total
number of molecules at t = 0, i.e., (nA0+nBO), as nT0, Eq. (2)
can be rewritten as
nT = NA t)+ NBt) (3)

If NA(t) is selected as the independent random variable, NB(t) can
be expressed in terms of NA(t), from the above equation, as
NB (t) = nTO -NA (t) (4)

whose realization is naturally
nB nTO -nA

Consequently, the state space for the current system is all
the possible numbers of molecules of A, i.e., {nTO, (nT - 1),
..., 2, 1,0}.

The derivation of the master equation of a stochastic process
has been extensively elaborated in our previous contribution
to this journal.12s5 For the birth-death process of interest, the
probability balance around any arbitrary state nA of the inde-
pendent random variable, NA(t), gives rise to1" 12]
-P(nA;t)= W+ (nA -1;t)p(n -1;t)
+W (n +l;t)p(nA +l;t)
-[W (nA;t)+W (nA;t) p(nA;t) (6)

This derivation is also detailed inAppendix A. The term, p(nA;
t), in the above equation denotes the probability of nA mol-
ecules being present at time t. On the basis of the first-order
rate law, the intensity of birth, W+(nA; t), and the intensity of
death, W (nA; t), have been defined as"'
W+ (nA;t) = = k2nB = k2 (nTO -nA) (7)

W (nA;t) dnA kn, ,A

respectively. It is worth noting that the intensities of birth and
death can be defined on the basis of not only linear but also
nonlinear rate laws, thereby giving rise to a variety of linear
and non-linear stochastic models. In fact, the formulation of
stochastic models for processes obeying nonlinear rate laws
has been presented in some of our earlier contributions.J26 28]

The solution of the master equation, Eq. (6), renders it pos-
sible to obtain the mean and higher moments about the mean
of the random variables, N (t) and NB(t); see Appendix B.
Chemical Engineering Education

Among these higher moments, the second moment about the mean, i.e., the variance, is
of utmost importance: It signifies the fluctuations, or scatterings, of the random variable
about its mean.[29]
Eq. (6) in conjunction with Eqs. (7) and (8) yields the expression for the mean of
N,(t), E[NA(t)] or (NA(t)), as

(NA (t))= n AO (k +k2)-k2 exp[-(k +k2)t +k2 (9)
(k +k2 nTO
Note that in the particular case where nA0 = nT0, i.e., where nBO = 0, the above expression
reduces to30, 31]
(NA (t))= - TO {k exp[-(k, +k2)t +k2, (10)
(kl +k2)
In addition, the mean of NB(t), E[NB(t)] or (NB(t)), is, from Eqs. (4) and (9),

(NB (t))=nTO 1- I AO (k +k2)-k2 exp[-(k +k2)t +k2 (11)
(k, +-k2) nTO J
Similarly, the expression for the variance of NA(t), Var[NA(t)] or Ao2(t), is evaluated as
( (t)=k nTO 1-exp[-(k +k2)t]}
2 (k +k2)
[nA k + k2) + k2 exp[-(k+ k2) t + klk (12)
For the case, nA0 = n, this equation becomes[301

S(kt)= k 1 {1-exp-(k1 +k2)t]}NA(t)) (13)
(7, (t) -*kj +k2

where (NA(t)) is given by Eq. (10). The standard deviation, GA(t), is the square root of
the variance, A2(t); thus, from Eq. (12),
A jt nw = TO (1- exp[-(k,1 +k2) t]1/2
t (k + k2)
nAO k1 -k )+k2 exp[-(k1+k2)t +klk2 (14)

From Eqs. (4) and (12), the variance of NB(t), Var[NB(t)] or (B2(t), is

St) nTO -exp[-(k +k2)t}
S(k , +k2)

SnAO j(ki-k) + k exp[-(k +k2)t +kk2 = (t) (15)

From this expression, oB(t) is given by
TB (t)= nTO {1-exp[-(k k2)t}/2
(k, + k2)

2 22
[T ( k2) +k2 exp[-(k +k2)t+ klk2

Vol. 42, No. 1, Winter 2008


Randomly or



processes or

systems are

ubiquitous in

the chemical

or allied


In general,
the stochastic
simulation of
kinetics can
focus on the
evolution of
the number
of entities
comprising a
reacting system
or the spatial
among these

The coefficient of variation, CV(t), signifies the relative fluctuations of a random variable
about its mean[32]; it is computed as the ratio between the standard deviation, o(t), and the
mean, (N(t)). From Eqs. (9) and (14), the coefficient of variation of NA(t), i.e., CVA(t), is
CVA (t)
A (t)
(NA (t))

nTk0 1- exp
(k +k2)

(k + k2)t 1/2 nA k k2 -k +k
SnT 1 2 *+k

TOk, + k) n k -k2) k2 exp -(k,
(k l -k 2) n TO 2

CVA (t)= 1 -exp k

exp -(k, + k2)t]

kk 1/2
kL k2


k2')t /2

no (nAk k -k) +k exp[- (k1 +k2)t +kk2
nAO (k -k2)-k2 exp[-(k +k2)t +k2

By substituting the expression for (NA(t)) as given by Eq. (9), the above equation can be
rewritten as

CVA (t)= A ) [ k2 jexp
(NA (t)) k, +k,2

(k,+ k2)t] -[k 1 (NA t))
k, +k2

Sk2 J nHAoexp[-(k +k2)t]

Sk nA exp-2(k +k2)t}1/2
k, +_k2 j

When n0 = nA0, this expression becomes

CVA( t)= 1exp- k
k, +-k2)

2) (NA (t))1/2
1 (NA (t))

where (NA(t)) is given by Eq. (10). In their contribution, Martinez-Urreaga and his col-
laborators'11 indicate that the relative fluctuations of NA(t) about its mean (NA(t)), i.e., the
coefficient of variation, CVA(t), are approximately of the order of ((NA(t)))1/2N (t), or(
NA(t) 1/2 as revealed through Eqs. (18) and (19).
The master equation of the birth-death process, Eq. (3), is stochastically simulated by
means of the Monte Carlo method. The stochastic simulation of chemical kinetics by this
method has been known for many years.17, 18 21 33-36] A number of the works devoted to it,
however, fail to deal with the crux of any stochastic simulation, which is its capability to
estimate not only the mean but also higher moments about the mean, especially the vari-

Chemical Engineering Education

ance.[25] According to Haseltine and Rawlings,[37 "... stochastic
simulation of chemical kinetics has received an increased amount
of attention from the modeling community. .. ." In general,
the stochastic simulation of chemical kinetics can focus on the
temporal evolution of the number concentrations of entities
comprising a reacting system or the spatial interactions among
these entities. The former can be accomplished by resorting to
the classical formulation of the Monte Carlo method,131 and the
latter can be carried out by means of Kinetic Monte Carlo (KMC)
methods, which take into account the spatial distribution of the
entities at the atomic or molecular scale.[38 39] In this work, the
Monte Carlo method is deployed for stochastically simulating
the temporal evolution of the numbers of molecules of A and B
for the simple reversible reaction; thus, the simulation of spatial
interactions among these molecules is not considered. The two
basic procedures to implement the method, one resorting to
the event-driven approach10 15-23] and the other resorting to the
time-driven approach,[14 15] are described in detail; these two ap-
proaches differ in the manner of updating the simulation clock
of the process's temporal evolution.
Event-Driven Approach
The event-driven approach advances the simulation clock by a
random waiting time, T, which has an exponential distribution.10, 19]
No event takes place during the time interval, (t, t + T), and only
one event occurs at the end of this time interval at which the
state of the system is specified by the probability of transition
corresponding to each event. The Gillispie algorithm[10 17 18] com-
prises a series of steps to perform the Monte Carlo simulation of
Markov processes by the event-driven approach. This algorithm
has been recently extended[351 to accelerate the computing speed
of the original formulation.17,18] Herein, Gillespie's original algo-
rithm, as presented by Martinez-Urreaga and his collaborators, '1
has been revised and augmented to include the evaluation of
the means, variances, and standard deviations of the birth-death
process of interest as described below.
Step 1. Define the total number of molecules per unit
volume, nTO, the total number of simulations, or
trajectories, Zf, and the length of each simulation,
tf. Initialize the simulation counter as Z <- 1.
Step 2. Initialize clock time t, data-recording time 6,[31]
the realizations of NA(t) and NB(t) at time t and
simulation Z, i.e., nAz(t) and nBz(t), respectively,
and the realizations of NA(6) and NB(6) at time 6
and simulation Z, i.e., nA z(6) and nB z(6), respec-
tively, as follows:
00 to
nAZ (to) nA0
nB.Z (to) TO- nAZ (t)
nA.Z (0) nAZ(t0)
nBZ ) nBZ to)

Step 3. Sample a realization u from the uniform random
variable, U, on interval (0, 1). Estimate T according
to the following expressioni10;


where W [nAZ(t);t] = k2[n T-nA(t)] and
W [nA,(t);t] = knAz(t) are the intensities of birth
and death, Eqs. (7) and (8), respectively; see Ap-
pendices C and D.

Step 4. Advance clock time as t <- (t + T).

Step 5. If (6 < t), then continue to the next step; otherwise,
continue to Step 8.

Step 6. Compute the sample means, variances, and standard
deviations at time 6 as follows:
a. Record the value of realizations at 6:

nAZ () nA,Z t- T)
nB.z () nB.Z (t- T)

b. Store the sum of realizations at 6:
A.z A nAZ (0)
B,z () - nBZ ()

c. Store the sum of squares of realizations at 6:
4A,Z Zn) n2
Zn,z 0) nB ,z (0)

d. Store the square of sum of realizations at 6:


4'BZK0) ZnB

[ A,Z (0)2

B,Z(0) 2

e. Compute the sample means at 6[13 19 20]:
1 z 1
Z z-1 Z
1 z 1
mBZ )z () -
MBZ()Z-EB,Z() =-BZ()
Z Zl Z



f. If 1 < Z < Zf, then compute the sample variances
and standard deviations at 0[13 19, 20]:

Vol. 42, No. 1, Winter 2008

-W+ n (t);t +W n t);t n(l-u)
- W [nA, (t); t+W_[n,(t);t]

sz (e)

- -� nAZ AZZ) - Z AZe)-1 )
z- 1) z znz O (Z- 1) Z

sB,z (0)

Z1 Z-
(Z- n ))

- nz(O)
z Z1

(Z 1) z z (


SAz(0) _ s zO (0)1/2

SBZ 0) _ 5s2z( ) ( 2

Step 7. Advance 6 by a conveniently small AO as 6 <- (6+A6). If (6
< t), then return to Step 5; otherwise, continue to Step 10.
Step 8. Determine the state of the system at the end of time interval
(t, t +T). To accomplish this, sample a realization u' from the
uniform random variable, U, on interval (0, 1) and compute
the probability of transition corresponding to the birth event,
w [nA,(t);t], as follows:
W+ [nA, (t);t (2
w [nAZ (t); tl= 27)
W+ nA, (t);t+ W [nA,Z (t); tK
If {0 population of molecules of A increases by one; thus,
nA, (t) [nA, (t-- T)
nBZ (t) [�nT- nAZ (t)
nAZ() nAz (t)
nB () nZ (t)
Otherwise, a death event occurs, i.e., the population of mol-
ecules of A decreases by one; thus,
nA, (t) nA,Z (t-- )--
nB.Z (t) [nT - nA, (t)
nAZ() nA,Z (t)
nBZ (0) nB (t)

Step 9. Repeat Steps 3 through 8 until tf is reached.
Step 10. Update simulation counter as Z - (Z + 1).
Step 11. Repeat Steps 2 through 10 until Zf is reached.
Given in Appendix E is the computer code in MATLAB for performing
Monte Carlo simulation on the basis of the Gillespie algorithm via the
event-driven approach as presented above.
Time-Driven Approach
As briefly indicated at the outset of this contribution, the time-driven

approach[14 15] differs from the event-driven ap-
proach: It advances the simulation clock by a
fixed time increment of At, which is sufficiently
small so that at most one or no event occurs
during time interval (t,t + At). At the end of this
interval, the state of the process is determined
by the probability of transition corresponding to
each event. The series of steps for implementing
Monte Carlo simulation of Markov processes
in general and the birth-death process of inter-
est in particular by the time-driven approach is
given below.
Step 1. Define the total number of molecules
per unit volume, nTO, the total num-
ber of simulations, or trajectories, Z,,
and the length of each simulation, t,.
Initialize the simulation counter as Z
<- 1. Compute time increment At
as follows:
At=1 (28)
WMt nA (hA;t) + wM (nA; t) I

where WM(nA;t) and WM(nA;t) are the
maximum possible values of the inten-
sities of birth and death, respectively.
Step 2. Initialize time t and the realizations
of NA(t) and NB(t) at time t and
simulation Z, i.e., nA,(t) and nz(t),
respectively, as follows:
t- to
nAZ (t) nAO
nB.Z (to) nTO - A.Z (t)

Step 3. Compute the sample means, vari-
ances, and standard deviations at
time t as follows:
a. Record the values of realizations at t:
nAZ (t) nA.Z
nB. (t) - nB.Z
b. Store the sum of realizations at t:

A.Z (t) - nAZ (t)

c. Store the sum of squares of realiza-
tions at t:
QA.Z (t) Z- T Z (t)
^B.Z (t) n- n.z (t)

Chemical Engineering Education

d. Store the square of sum of realizations at t:


Bz (t) n BZ(t)

[ Kz(t)2

. Z(t)12

e. Compute the sample means at t:

mA.z (t) nA.z () = -7.z (t
7=1 Z

z 1
mBz(t) ZnBz(t)= 7, (t)
7-1 Z

Step 7. Repeat Steps 3 through 6 until
tf is reached.
Step 8. Update the simulation counter
as Z <- (Z + 1).
Step 9. Repeat Steps 2 through 8 until
Zf is reached.
Given in Appendix F is the computer code
in MATLAB for performing Monte Carlo
(29) simulation via the time-driven approach as
presented above.


f. If 1 < Z < Z,, then compute sample the variances and standard
deviations at t:

z (tt)

Z-1) nZzt(t)

sB.z (t)
n zt)
S(z-i { B
(Z - 1) z t1 )

1 z 2 I
Z 1aTz1 = (Z-1)L

- i nB (t)
Z z-1

1 )
_ 1) {^B.z(t)-
(z-i~l -

-Z Az(t) (31)

-I Z (t) (32)

SAz (t) sz ( t)l2 (33)

SBz(t) S2z(t)2 (34)

Step 4. Advance time as t <- (t + At).
Step 5. Estimate the probabilities of transition corresponding to the birth
and death events as {W+[na.z(t);t]}At and {W [naz(t);t]+W
[naz(t);t]}At, respectively.
Step 6. Determine the state of the system at the end of time interval (t,t + At).
To accomplish this, sample a realization u from the uniform random
variable U on interval (0, 1) and compare it with the probabilities of
transition for the birth and death events. If u < [W (nA;t)]At, then a
birth event occurs; thus,
nA, (t) nAz (t- At) + 11

nB (t) [nTO - nAZ (t)

If u < [W (nA;t) + W (nA;t)]At, then a death event occurs; thus,
nAz (t) [nAz (t - At) - 1
nB (t) [nT - n-A. (t)

Otherwise, no event occurs; thus,
nAz (t) nz (t- At)
nB.Z (t) nT - nAz (t)

Vol. 42, No. 1, Winter 2008

The means, variances, and standard de-
viations of the random variables, NA(t) and
NB(t), characterizing the simple reversible
reaction have been computed with the values
for k1 and k2, and the two sets of values for
nA0 and nB0, given by Martinez-Urreaga and
his collaborators."1 Thus, k1 is 4 - 103 s 1
and k2 is 1 - 103 s-1; the first set of values
for nA0 and n0O comprises 175 molecules of
A and 25 molecules of B, and the second
set comprises 3,500 molecules of A and 500
molecules of B.
The means, (NA(t)) and (NB(t)), have
been computed from the stochastic model
according to their corresponding analytical
expressions, Eqs. (9) and (11); the results
are illustrated for the first set of nA0 and nB0
in Figure 1 (next page) and for the second
set of nA0 and nB0 in Figure G-1 in Appendix
G. In these figures, the numbers of molecules
of A and B obtained from the solution of the
deterministic model as presented by Mar-
tinez-Urreaga and his collaborators1lI are
superimposed for comparison. Clearly, the
values of (NA(t)) and (NB(t) )obtained from
the stochastic model are in excellent accord
with the corresponding deterministic values,
thereby verifying that the mean component
of the stochastic model for the birth-death
process is equivalent to the deterministic
The standard deviations, oa(t) and o,(t),
signify the variability attributable to the in-
ternal or characteristic noises of the process
as predicted by the stochastic model. The
values of oa(t) and o (t) have been computed
from Eqs. (14) and (16), respectively; they
are plotted in Figures 1 and G-1 as the stan-
dard deviation envelopes, [(N (t)) � o(t)]
for N (t) and [(NB(t)) � o,(t)] for NB(t),

0 200 400 600 800 1000
Figure 1. Temporal evolution of the numbers of molecules
of A and B per unit volume with noA = 175 and nBo = 25:
The solutions from the stochastic model are compared
with those from the deterministic model.

respectively. Note that the variability or dispersion of the
process is more pronounced in Figure 1 than in Figure G-1:
As indicated by Martinez-Urreaga and his collaborators"11
and earlier in the current contribution, the relative extent of
fluctuations of a random variable N(t) about its mean (N(t)) is
approximately of the order ((N(t))' - (Nio), or (N(t) "12, i.e.,
it is inversely proportional to the square root of the popula-
tion size. In other words, the smaller the population size, the
greater the extent of variability or dispersion of the process
about its mean and vice versa.
The sample means, mA z(6) and m. z(6), have been comput-
ed through the Monte Carlo method on the basis of the event-
driven approach by averaging the realizations of NA(t) and
NB(t) from Z successive simulations. The results are illustrated
for Z = 2 in Figure 2 and Z = 50 in Figure 3 for the first set of
nAo and n B; the values of (NA(t)) and (NB(t)) computed from
the stochastic model as well as the numbers of molecules of A
and B obtained from the solution of the deterministic model are
superimposed in both figures for comparison. Clearly, mA z(6)
and m z(6) approach (NA(t)) and (NB(t)), respectively, as the
number of simulations, Z, varies from 2 to 50. The results
are in line with the Weak Law of Large Numbers (WLLN),
stating that the sample mean approaches the population mean
as the size of the sample becomes large.[29] In stochastically
simulating the birth-death process, the size of the samples,
i.e., the numbers of realizations of NA(t) and NB(t), increases
as the number of simulations, Z, increases. Similarly, m,(6))
and mn z(6) are essentially identical to (NA(t)) and (NB(t)),
respectively, when they are computed by averaging 50 simula-
tions, i.e., Z = 50, for the second set of nA0 and n0O; see Figure
G-2 in Appendix G. The sample standard deviation envelopes,
[mA, (6) � sAz(6)] for NA(t) and [mBz(6) � SBz(6)] for NB(t),
are also plotted in Figures 2 and 3 as well as in Figure G-2 in
Appendix G; they are essentially identical to [(NA(t)) o,(t)]

0 200 400 600 800 1000

Figure 2. Temporal evolution of the means and standard
deviation envelopes of the numbers of molecules of A
and B per unit volume with n o = 175 and nBo = 25: The
average of two Monte Carlo simulations via the event-
driven approach is compared with the analytical solu-
tions resulting from the stochastic model; the solution of
the deterministic model is represented by the solid lines.

0 200 400 600 800 1001
_t or e, s

Figure 3. Temporal evolution of the means and standard
deviation envelopes of the numbers of molecules of A
and B per unit volume with n o = 175 and nBo = 25: The
average of 50 Monte Carlo simulations via the event-driv-
en approach is compared with the analytical solutions
resulting from the stochastic model; the solution of the
deterministic model is represented by the solid lines.

and [(NB(t)) � OB(t)], respectively, for Z = 50.
Figures G-3 and G-4 in Appendix G exhibit mAz(t) and
mn.z(t) as well as [mA.z(t) � s .z(t)] and [mn.z(t) �_ S.z(t)],
which have been computed through the Monte Carlo method
by resorting to the time-driven approach. For the first set of
nAO and n B, the results with Z = 2 are presented in Figure
G-3, and those with Z = 50 are presented in Figure G-4; they
are compared with the solutions of the deterministic and
stochastic models. Note that mA.z(t) and mBz (t) in Figure G-4
Chemical Engineering Education

are equivalent to (NA(t)) and (NB(t)), respectively,
for Z = 50; this is also the case for the second set
of nA0 and nB0 as illustrated in Figure G-5 in Ap-
pendix G.
Table 1 lists the computational times, required by
the Monte Carlo method via the event-driven and
time-driven approaches, for estimating the sample
means, variances, and standard deviations of the
random variables, NA(t) and NB(t), characterizing
the simple reversible reaction. To perform the Monte
Carlo simulations by resorting to these approaches,
the MATLAB codes given in Appendices E and F
were implemented on a PC (Pentium IV 3.0 GHz
with 512 MB RAM;Windows XP). Clearly, the time
required by the event-driven approach was shorter
than that required by the time-driven approach.
Figures 2 and 3 in conjunction with Figures G-2
through G-5 in Appendix G reveal a fundamental
distinction between the variations attributable to
the inherent or internal noises of a random process
and those attributable to the numerical method to
stochastically simulate it. The former can only be
quantified in terms of higher moments about the
means of the random variables characterizing the
process, especially, the variances or standard devia-
tions. The latter arise from the generation of samples
on the basis of random numbers; these variations
tend to vanish as the number of random samples
increases, i.e., as the number of simulation runs
becomes large.
The efficacy of the two approaches for stochas-
tic simulation via the Monte Carlo method, i.e.,
the event-driven and time-driven approaches, are
further illustrated with an example of photoelec-
trochemical disinfection of bacteria whose rate
of decaying is assumed to be of the first-order, or
linear. In contrast to a chemical reaction involving
molecules that are discrete and microscopic in size;
the disinfection of bacteria deals with microorgan-
isms which are also discrete but mesoscopic in size.
In general, the inherent fluctuations of a system
comprising mesoscopic entities are appreciably
more pronounced than those of a system comprising
microscopic entities, especially, when the number
concentrations of such entities are minute.
Naturally, the system of interest is the population
of bacteria per unit volume. Catalyzed photoelec-
trochemically, these bacteria die one at a time;
moreover, the bacteria have ceased to grow and do
not reproduce throughout the deactivation. This sys-
tem constitutes a special instance of the birth-death
processes in which the intensity of birth is absent;
hence, it is termed the pure-death process.[10 11 25]
Vol. 42, No. 1, Winter 2008

In modeling the photoelectrochemical disinfection of bacteria as
a pure-death process, the number of live bacteria at time t, Nc(t), is
taken as the random variable of the process; a realization of Nc(t) is
denoted by nc. All the possible numbers of bacteria are the states of
the process, and the collection of these numbers, {nco, (n0 - 1), ...,
2, 1, 0}, constitutes the state space, where nco is the initial number
of bacteria, i.e., N,(0) or nc at t = 0. For the pure-death process, the
probability balance around any arbitrary state nc of the independent
random variable, Nc(t), leads to"11

dp(nc;t)=W (nc +1;t)p(nc +1;t)-W (nc;t)p(nc;t)


which is the master equation of the process11, 25]; see Appendix A. On
the basis of the first-order rate law, the intensity of death, W (nc;t) ,
in the above expression is

W (nc;t) dnc - k3nc

where k3 is the first-order rate constant. For the pure-death process,
the expressions of the mean, E[Nc(t)] or (Nc(t)), and the variance,
Var[Nc(t)] or oc2(t), are obtained from Eqs. (35) and (36) as[25 40]
(Nc (t))= nco exp(-k3t) (37)


7c (t)= nco exp(-k3t)[1- exp(-k3t)],


respectively; see Appendix B. Obviously, oc2(t) can be related to
(Nc(t)) as

C( (t) = [- exp(-k3t)l(Nc (t))

The standard deviation, oc(t), is, from Eq. (38),

(Tc (t)= n0 {exp(-k3t)1- exp(-k3t) }12



Computational Times for Estimating the Sample Means,
Variances, and Standard Deviations of the Random Variables,
N (t) and NB(t), Characterizing the Simple Reversible Reaction
by the Monte Carlo Method via the Event-Driven
and Time-Driven Approaches
Initial values Event-driven approach* Time-driven approachW
Z=2 Z=50 Z=2 Z=50
n = 175 < s < s Is Is
nBO= 25
nA = 3,500 1 s 9s 5s 17s
nB0= 500
* In the event-driven approach, the sample means, variances, and standard
deviations have been recorded every two seconds, i.e., AO = 2 s.
� In the time-driven approach, the fixed time increment, At, has been com-
puted according to Eq. 28 in the text with WM(nA;t)=k2nT = k2(nAO+nBO)
and W M(nA;t)=k2nAO. Thus, At = 1.111 s with nA = 175 and nB = 25,
and At = 0.056 s with n = 3,500 and nB = 500.


or, from Eq. (39),

(c (t)= l-exp (-k 3t)12 (Nc (t))1/2 (41)

From Eqs. (37) and (40), the coefficient of variation, CVc(t),
is computed as

1 1-exp(-k3t) 242)
C exp(-k3t)

which can be rewritten as

CVc (t)= [1- exp(-k3t) 2 (Nc (t))2 (43)
(Nc (t))

Note that for the pure-death process the relative fluctuations
of Nc(t) about its mean (Nc(t)), as given by CVc(t), are ap-
proximately of the order of ((N(t) ))1 2/ (Nc(t)), or (Nc(t)) 12,
as indicated by Martinez-Urreaga and his collaborators. 1]
The temporal mean of the model, Eq. (37), has been re-
gressed on the experimental data for the photoelectrochemical
disinfection of E.coli[24] near the termination period of the
process. By identifying 140 cells per milliliter as the initial
population of bacteria, nco, the nonlinear regression by means
of the adaptive random search procedure41 42] has resulted in 0.091
min1 for the value of k3. With these values of no and k3, (Nc(t)
) and oc(t) have been computed from Eqs. (37) and (40), respec-
tively. Figure 4 illustrates the resultant values of (Nc(t))
and [(Nc(t)) � oc(t)] in conjunction with the experimental
data. 24] In addition, the sample means, mcz(6) and mc.z(t), and
standard deviation envelopes, [mcz(0) � sc.z(6)] and [mc.z(t) �
scz(t)], have been computed via the Monte Carlo method by
resorting to both the event-driven and time-driven approaches,
respectively. The results from averaging 50 simulations, i.e.,
Z = 50, are presented in Figure 4 for comparison. Note that
the results from Monte Carlo simulation closely approximate
the analytical results from the stochastic model. As expected,
the deviations or fluctuations of the experimental data[24] are
substantially more pronounced than those predicted by the
stochastic model: The overall deviations of the experimental
data include those attributable not only to the internal noises
of the process but also to the external noises arising from
inevitable manipulation errors and instrumental limitations
that can never be totally suppressed.

The contribution of Martinez-Urreaga and his collabora-
torsM11 on the stochastic simulation of chemical reactions,
specifically, a simple reversible chemical reaction, has been
complemented and augmented in two aspects. First, a sto-
chastic model for a simple reversible chemical reaction has
been derived based on the first-order, i.e., linear, rate law as
a birth-death process. The master equation arising from the
model has rendered it possible to analytically obtain the ex-

pectedmeans and variances of the process. Second, the master
equation of the process has been stochastically simulated
through the Monte Carlo method via the time-driven approach
in addition to the event-driven approach on the basis of the
Gillespie algorithm. The process's means and variances from
numerical simulation are in accord with the solutions of the
deterministic and stochastic models. To further reveal their
efficacy, the two approaches for stochastic simulation by the
Monte Carlo method are illustrated with the photoelectro-
chemical disinfection of bacteria obeying the first-order rate
law. The mean computed from stochastic simulation are in
line with the available experimental data as well as with the
analytical results derived from the corresponding stochastic
model. As expected, the fluctuations of the data about their
mean are large when compared with those computed from
the stochastic simulation.

This is contribution No. 07-123-5, Department of Chemical
Engineering, Kansas Agricultural Experiment Station, Kansas
State University, Manhattan, KS 66506.

E[NA(t)] mean, expected value, or the first moment of the
random variable, Na(t)
E[NB(t)] mean of the random variable, NB(t)
E[Nc(t)] mean of the random variable, Nc(t)

t or 0, min

Figure 4. Temporal evolution of the mean
and standard deviation envelope near the tail-end
of photoelectrochemical disinfection of a population
ofE. coli per unit volume: The average of 50
Monte Carlo simulations by resorting to the
event-driven and time-driven approaches
is compared with the analytical solutions resulting
from the stochastic model;
the initial population of bacteria,
nco, is 140 cells per milliliter.

Chemical Engineering Education






p(nA; t)

forward reaction-rate constant, s 1
reverse reaction-rate constant, s 1
first-order rate constant, min
random variable representing the number of mol-
ecules of A at time t
random variable representing the number of mol-
ecules of B at time t
random variable representing the number of live
bacteria at time t
realization of the random variable, NA(t)
realization of NB(t)
realization of Nc(t)
number of molecules of A at t = 0, NA(0)
number of molecules of B at t = 0, NB(0)
total number of molecules at any arbitrary t, (n A
+ nBO)
number of live cells at t = 0, Nc(0)
probability of nA molecules being present at time

p(nc; t) probability of nc live bacteria being present at
time t

W (nA; t)

uniform random variable on interval (0,1)
realization of U
variance of the random variable, NA(t)
variance of NB(t)
variance of Nc(t)
intensity of birth for the birth-death process in
state nA at time t

W (nA; t) intensity of death for the birth-death process in
state n at time t
W (nc; t) intensity of death for the pure-death process in
state nc at time t
Greek Letters
0 data-recording time in the event-driven approach
oA2(t) Var[NA(t)]
oB2(t) Var[NB(t)]
oc2(t) Var[Nc(t)]
GA(t) standard deviation of the random variable, NA(t)
OB(t) standard deviation of NB(t)
oc(t) standard deviation of Nc(t)
T random waiting time in the event-driven approach


1. Martinez-Urreaga, J., J. Mira, and C. Gonzalez-Fernandez, "Introducing
the Stochastic Simulation of Chemical Reactions Using the Gillespie
Algorithm and MATLAB," Chem. Eng. Educ., 37, 14 (2003)
2. Fan, L.T., S.H. Hwang, S.T. Chou, and R. Nassar, "Birth-Death Model-
ing of Deep Bed Filtration: Sectional Analysis," Chem. Eng. Comm.,
35, 101 (1985)
3. Fan, L.T., B.C. Shen, and S.T. Chou, "Surface-Renewal Theory of

Vol. 42, No. 1, Winter 2008

Interphase Transport: A Stochastic Treatment," Chem. Eng. Sci., 48,
3971 (1993)
4. Chou, S.T., L.T. Fan, and R. Nassar, "Modeling of Complex Chemical
Reactions in a Continuous-Flow Reactor: A Markov Chain Approach,"
Chem. Eng. Sci., 43, 2807 (1988)
5. Fox, R.O., and L.T. Fan, "Stochastic Modeling of Chemical Process
Systems. Part I: Introduction." Chem. Eng. Educ., 24, 56 (1990)
6. Woodle, G.R., and J.M. Munro, "Particle Motion and Mixing in a
Rotary Kiln," Powder Tech., 76, 241 (1993)
7. Diaz, E., J. Szepvolgyi, and J. Gyenis, "Modeling and Dynamic Simula-
tion of the Stochastic Behavior of Bulk Solid Mixing," Hungar. J. Ind.
Chem., 25, 115 (1997)
8. Alonso, M., and FJ. Alguacil, "Stochastic Modeling of Particle Coat-
ing," AIChE J., 47, 1303 (2001)
9. Berthiaux, H., and V. Mizonov, "Applications of Markov Chains in
Particulate Process Engineering: A Review," Can. J. Chem. Eng., 82,
1143 (2004)
10. Gillespie, D.T., Markov Processes-an Introduction for Physical
Scientists, Academic Press, San Diego, pp. 48, 60, 226, 328, 330, 375,
11. van Kampen, N.G., Stochastic Processes in Physics and ( I ...
North-Holland, Amsterdam, pp. 55-58, 96-97, 134-136, 139, 163
12. Oppenheim, I., K. E. Shuler, and G. H. Weiss, Stochastic Processes in
Chemical Physics: The Master Equation, The MIT Press, Cambridge,
MA, pp. 53-61 (1977)
13. Sobol', I.M., A Primerfor the Monte Carlo Method, CRC Press, Boca
Raton, FL, pp. i-1, ix, 15-16, 42 (1994)
14. Rod, V., and T. Misek, "Stochastic Modeling of Dispersion Formation
in Agitated Liquid-Liquid Systems," Transactions of the Institution of
Chemical Engineers, 60, 48 (1982)
15. Rajamani, K., WT. Pate, and D.J. Kinneberg, "Time-Driven and
Event-Driven Monte Carlo Simulations of Liquid-Liquid Dispersions:
A Comparison," Ind. Eng. Chem. Fund., 25, 746 (1986).
16. Kendall, D.G., "AnArtificial Realization of a Simple 'Birth-and-Death'
Process," J. Roy. Stat. Soc. B, 12, 116 (1950)
17. Gillespie, D.T., "A General Method for Numerically Simulating the
Stochastic Time Evolution of Coupled Chemical Reactions." J. Comput.
Phys., 22, 403 (1976)
18. Gillespie, D.T., "Exact Stochastic Simulation of Coupled Chemical
Reactions." J. Phys. Chem., 81, 2340 (1977)
19. Shah, B.H., J.D. Borwanker, and D. Ramkrishna, "Monte Carlo Simula-
tion of Microbial Population Growth," Math. Biosci., 31, 1 (1976)
20. Shah, B.H., D. Ramkrishna, and J.D. Borwanker, "Simulation of
Particulate Systems Using the Concept of the Interval of Quiescence,"
AIChE J., 23, 897 (1977)
21. Rao, C.V., andA.P Arkin, "Stochastic Chemical Kinetics and the Quasi-
Steady-State Assumption: Application to the Gillespie Algorithm," J.
Chem. Phys., 118, 4999 (2003)
22. Cao, Y., D.T. Gillespie, and L.R. Petzold, "Accelerated Stochastic
Simulation of the Stiff Enzyme-Substrate Reaction," J. Chem. Phys.,
123, 144917(2005)
23. Ullah, M., H. Schmidt, K.-H. Cho, and 0. Wolkenhauer, "Deterministic
Modelling and Stochastic Simulation of Biochemical Pathways using
MATLAB," IEE Proc. Syst. Biol., 153, 53 (2006)
24. Harper, J.C., PA. Christensen, T.A. Egerton, T.P Curtis, and J. Gun-
lazuardi, "Effect of Catalyst Type on the Kinetics of the Photoelec-
trochemical Disinfection of Water Inoculated with E. Coli," J. Appl.
Electrochem., 31, 623 (2001)
25. Fan, L.T., A. Argoti Caicedo, S.T. Chou, and WY. Chen, "Stochastic
Modeling of Thermal Death Kinetics of a Cell Population: Revisited,"
Chem. Eng. Educ., 37, 228 (2003)
26. Duggirala, S.K., and L.T. Fan, "Stochastic Modeling of Non-Linear
Sieving Kinetics," Chem. Eng. Comm., 61, 59 (1987)
27. Fox, R.O., and L.T. Fan, "Application of the Master Equation to


Coalescence and Dispersion Phenomena," Chem Eng Sci., 43, 655
28. Chou, S., L.T. Fan, A. Argoti, R. Vidal-Michel, andA. More, "Stochastic
Modeling of Thermal Disinfection of Bacteria According to the Logistic
Law,"AIChE J., 51, 2615 (2005)
29. Casella, G., and R.L. Berger, Statistical Inference, Duxbury, Pacific
Grove, CA, pp. 54, 55, 59, 89-91, 233 (2002)
30. McQuarrie, D.A., "Kinetics of Small Systems. I," J. Chem. Phys., 38,
433 (1963)
31. Steinfeld, J.I., J.S. Francisco, and W.L. Hase, Chemical Kinetics and
Dynamics, Prentice Hall, Englewood Cliffs, NJ, pp. 97-102 (1989)
32. Pang, W-K., P-K. Leung, W.-K. Huang, and W Liu, "On Interval
Estimation of the Coefficient of Variation for the Three-Parameter
Weibull, Lognormal, and Gamma Distribution: A Simulation-Based
Approach," Eur. J. Oper. Res., 164, 367 (2005)
33. Dixon, D.A., and R.H. Shafer, "Computer Simulation of Kinetics by
the Monte Carlo Technique," J. Chem. Educ., 50, 648 (1973)
34. Moebs, W.D.C., "A Monte Carlo Simulation of Chemical Reactions,"
Math. Biosci., 22, 113 (1974)
35. Gillespie, D.T., "Approximate Accelerated Stochastic Simulation of
Chemically Reacting Systems," J. Chem. Phys., 115, 1716 (2001)
36. Rawlings, J.B., and J.G. Ekerdt, Chemical ReactorAnalysis andDesign
Fundamentals, Nob Hill Publishing, Madison, WI (2002)
37. Haseltine, E.L., and J.B. Rawlings, "Approximate Simulation of
Coupled Fast and Slow Reactions for Stochastic Chemical Kinetics,"
J. Chem. Phys., 117, 6959 (2002)
38. Drews, T.O., A. Radisic, J. Erlebacher, R.D. Braatz, PC. Searson, and
R.C. Alkire, "Stochastic Simulation of the Early Stages of Kinetically

Limited Electrodeposition," J. Electrochem. Soc., 153, 434 (2006)
39. Mastny, E.A., E.L. Haseltine, and J.B. Rawlings, "Stochastic Simulation
of Catalytic Surface Reactions in the Fast Diffusion Limit," J. Chem.
Phys., 125, 194715 (2006)
40. Spilimbergo, S., and D. Mantoan, "Stochastic Modeling of S. Cerevi-
siae Inactivation by Supercritical CO2," Biotechnol. Prog., 21, 1461
41. Fan, L.T., H.T. Chen, and D. Aldis, "An Adaptive Random Search
Procedure for Large-Scale Industrial and Process Systems Synthesis,"
Proceedings of the Symposium on Computers in Design and Erection
of Chemical Plants, Aug. 31-Sep. 4, Karlovy Vary, Czechoslovakia,
pp. 279-291 (1975)
42. Chen, H.T., and L.T. Fan, "Multiple Minima in a Fluidized Reactor-
Heater System," AIChE J., 22, 680 (1976).
43. Taylor, H.M., and S. Karlin, Introduction to Stochastic Modeling,
Academic Press, San Diego, pp. 25, 357 (1998)
44. Jeffreys, H., and B. Jeffreys, Methods of Mathematical Physics, 3rd
Ed., Cambridge University Press, Cambridge, pp. 22-23 (1999)
45. McQuarrie, D.A., Mathematical Methods for Scientists and Engineers,
University Science Books, Sausalito, CA, pp. 4, 1053 (2003)


Seven appendices (Appendices A through G) can be viewed
or downloaded from the corresponding author's Web site at
, or can
be obtained as hard copies by writing to f.Lai kiu \.diu

Chemical Engineering Education

IM!1- class and home problems

Applications of



Using Mathematica

National Institute of Applied Sciences and Technology
Consider a single component characterized by its critical
pressure, temperature, acentric factor, Po T , and o. The
Peng-Robinson equation of state[2 4] is given by:
P = (1)
(V-b) [V(V+b)+b(V-b)]

b = 0.07780 R
(RTo )2
a= 0.45724 L 1+m - )

m= 0.37464 + 1.54226 - 0.26992 2. (5)
The Peng-Robinson EOS* is part of a family of equations
called cubic because the compressibility factor, Z = PV/RT,
is a solution of the following cubic equation derived from Eq.
(1) and written for a multi-component mixture:

3 + (- B)Z2 +(A- 3B2 -2B)Z + (-AB +B2

* Tunis, Tunisia

c c c c
A=EEy.yA, or EExxA,
I-1 ji1 I1 j1i
A = AAA 0 (1-k )
B= yB, or x,B,
A 1 = 1
P- Pr
A =0.45724 a- and B =0.07780-
ST2 ' T

S For each component, we define the reduced pressure and
temperature using P,, =P PP, and T = T/ , and a is
(4) identified by an equation similar to Eq. (3) - given for the
pure component case. The binary interaction parameter, k ,

Housam Binous is a full-time faculty
member at the National Institute of Applied
Sciences and Technology in Tunis. He
earned a Dipldme d'ingenieur in biotech-
nology from the Ecole des Mines de Paris
and a Ph.D. in chemical engineering from
the University of California at Davis. His
research interests include the applications
Sof computers in chemical engineering.

� Copyright ChE Division of ASEE 2008

B3= 0

*All these computations are available in the form of notebooks upon
requestfrom the author or the Wolfram Research Information Center.11
Vol. 42, No. 1, Winter 2008

The object of this column is to enhance our readers' collections of interesting and novel prob-
lems in chemical engineering. Problems of the type that can be used to motivate the student by
presenting a particular principle in class, or in a new light, or that can be assigned as a novel home
problem, are requested, as well as those that are more traditional in nature and that elucidate dif-
ficult concepts. Manuscripts should not exceed 14 double-spaced pages and should be accompanied
by the originals of any figures or photographs. Please submit them to Professor James O. Wilkes
(e-mail: wilkes@umich.edu), Chemical Engineering Department, University of Michigan, Ann
Arbor, MI 48109-2136.

is obtained from HYSYS 3.2, a major process simulator by
Aspen Technology, Inc.,[1] or assumed to be equal to zero if
the simulator is not available.
The liquid and vapor mole fractions are related by:
y, =K,x, with i= 1,2,.., C (11)
where K is the equilibrium constant.
This equilibrium constant is obtained using the t-
method as follows:

K,=- fori= ,2,...,C (12)
where the fugacity coefficients are:

(Zv 1) - In(Zv-B) A
B 241B
v, = exp 2 yA B 1 B1 (13)
-1 --- -- In -+(+-52)-
A B Z - B

The vapor-phase and liquid-phase fugacity coefficients,
which measure deviation from ideality, are defined from the
liquid-phase and gas-phase fugacities, fv andfl , as follows:
v f f l
S and , (14)
y,P x,P
The - method is a direct result of expressing vapor-
liquid equilibrium conditions (i.e., temperatures, pressures,
and fugacities are equal in both phases) and the definition of
fugacity coefficients.
A similar expression to Eq. (13) is obtained for the liquid-
phase fugacity coefficient, i,, by replacing the gas-phase
compressibility factor, Zv, with its liquid-phase counterpart,
ZI. These two compressibility factors are the largest and
smallest roots of Eq. (6), respectively. For specified values
of P and T, the smallest root, Zl , allows the determination of
the liquid molar volume while the highest root, Z , permits
the computation of the vapor molar volume. The intermediate
root has no physical significance.
Finally, one needs the following expression for the en-
thalpy departure function from ideality in order to compute

HD = RT(Z - 1) 1+
24B RT
L +(I + 2)B T d A(RT)2 (RT)2' (
LoZ+ 1- )B dT P PA (15)

Problem statement
A quaternary mixture, at 485 psia and 100 �F, is composed
of 0.41% hydrogen, 5.71% methane, 70.97% benzene, and

22.91% toluene. This mixture is sent to a stabilizer to remove
hydrogen and methane. The feed pressure is decreased adia-
batically from 485 psia to 165 psia by valve and pipe line
pressure drop. Find the vapor-phase fraction, the temperature
and other relevant variables.

This problem has been solved using a tedious iterative tech-
nique by Henley and Seader.[6] The unknowns in this problem
are the mole fractions in the two phases, the temperature,
the vapor phase fraction, and the compressibility factors. We
have 12 nonlinear algebraic equations to solve simultane-
ously. These equations are three equilibrium relations, three
component mass balances, two summation rules, two cubic
equations of the compressibility factors, the enthalpy balance,
and the famous Rachford and Rice equation given by:

zK, ) = 0 (16)
Ii+i vap(K -1)
where xvap is the fraction vaporized, V/F The Mathematica
commands, which allow the determination of the 12 un-
knowns, are based on the built-in function FindRoot and are
given by:

Eq[i_]:= y[i]
EQ[i_]:= x[i] == z[i]


== K[i] x[i]
/(1+ vap (K[i] - 1));

sEq[1],Eq[2],Eq[3], y[i]
sol = FindRootEQI,EQ[2,EQ[3

Sz[i](K[i] + (K[i]1))== 0,
1,> z[i] (K[i] - 1) /(1 (a+ vp (K[i] - 1)) == 0,

Z ^3 + (-1+ B,)Z A2 Z (A -3B A2- 2B,)
-Al B+B A 2 + B A3 == 0,
Zv 3 +(-1+ B) Zv 2 + Zv (A - 3Bv 2 - 2Bv)
-Av Bv + B 2 + B A3 == 0,
HF == vap HV + (1- vap) HL},{x[1],0.0004},
,vap0.01o , {y[1],0.1}, {y[2],0.85},{y[3],0.02},
{y[4], 0.001}, {Z,,0.05}, {Z,,0.95}, {T,540},
Maxlterations- > 1000] / /Chop
We find the following output,

x[l] - 0.00020583, x[2] - 0.0259441,
x[3] - 0.736087, x[4] - 0.237763
y[l] - 0.106105, y[2] 0.873206,
y[3] - 0.0185161,y[4] - 0.00217289,
tvap - 0.0367724,Z, - 0.0405552,1
Z, - 0.981151, T - 560.297

Chemical Engineering Education

The vapor-phase fraction vap and temperature are 0.0367
and 560.297 �R, respectively.


Problem statement
Methanol is used as a solvent in the sweetening of natural
gas (i.e., removal of carbon dioxide and hydrogen sulfide).
Loss of methanol by transfer from liquid phase to gas phase
can become an important economic consideration. Compute
the gas phase mole fraction of methanolP31 in treated natural
gas leaving the sour gas absorption unit vs. temperature at
a pressure of 145.038 psia. Treated natural gas is composed
of 70% methane, 10% nitrogen, 10% carbon monoxide and
10% hydrogen.

We perform several isothermal flash computations at vari-
ous temperatures. The problem's unknowns are the ten mole
fractions in the gas and liquid phases, the vapor phase fraction,
and the two compressibility factors. The available equations
are the five equilibrium relations, the two summation rules,
the five mass balances, the two cubic equation of state of the
compressibility factors, and the Rachford and Rice equation.
The calculation uses the built-in function FindRoot to solve
13 nonlinear algebraic equations simultaneously. The Math-
ematica commands, which allow the determination of the 13
unknowns, are the following:



D 0.001


350 400 450 500 550
Temperature in �R

Figure 1. Solubility of methanol in treated natural gas
vs. temperature at 145.038 psia.

Eq[i_]:= y[i] == K[i] x[i]
EQ[i_]:= x[i] == z[i] /(1+ vap (K[i] -1));
For[i = ,i < 15,
T= 1329.67,347.67,365.67,383.67,401.67,
491.67,509.67,527.67,545.67,563.67} [[i]l,

sol[i]= FindRoot[ Eq[l],Eq[2],Eq[3],Eq[4], y[i]== 1,

EQ[1],EQ[2],EQ[3],EQ[4], x[i]== 1,

Lz[i] (K[i]- 1)/(1 + vap (K[i] - 1))== 0,
Z1A 3 + (-I+ 1B) Zl A2 + Z, (A, -3B132 2B1)
-A, B, +B,^2 +B 3== 0,
ZA3+(-1+ B)ZA2+ Z (A -3B, A2 -2B,)
-Av B +BA 2 +B A3== 0
{x[1],x[1] /.sol[i -1]},{x[2],x[2] .sol[i-1]},
{x[3],x[3] .sol[i - 1]},
{x[4],x[4] .sol[i -1]},{x[5],x[5] /.sol[i -1]},
{y[1],y[1] /.sol[i -1]}, {y[2],y[2] .sol[i- 1]},
{y[3],y[3]/ .sol[i - 1]},
{y[4],y[4] .sol[i -1]}, {y[5],y[5] /.sol[i -1]},
{ vap vap/ .sol[i -1]}, {Z1, Z / .sol[i -1]},
{Zv,,v /.sol[i- 1]},
Maxlterations- > 1000] / Chop,Print[i],i + +}]
Figure 1 depicts the computed solubility of methanol
in treated natural gas vs. temperature at 145.038 psia. We
observe an increase in solubility at higher temperatures. We
see that these solubility values can become large enough to
present serious economic drawback due to excessive loss of
Problem Statement
Nitrogen and hydrogen react to form ammonia. This reac-
tion is favored by low temperatures and high pressure. Kinetic
considerations, however, oblige us to use high temperature.
Thus, reactors are operated at very high pressures to get a
reasonably high conversion of reactants. High gas-phase
pressures imply significant deviation from ideality and the
need to take into account the fugacity coefficients.[4] In fact,
the equilibrium constant depends on K as follows:
a a YNH, 1 X(2-X) 1
Ka= 05 15 K = K (17)
=1 YN2 YH2 p 1 - X 3(1 X) P
2 2

Vol. 42, No. 1, Winter 2008

where Kv is given by:

K=v , (18)
1.05 15

a, = 4,y, (19)
and X is the extent of reaction. Plot K and X vs. pressure at
temperatures ranging from 900 OR to 1440 OR.
The unknowns in this type of problems are five: the mole
fraction in the gas phase, the extent of reaction, and the gas
phase compressibility factor. Once again, the calculation
uses the built-in function FindRoot to solve five nonlinear
algebraic equations simultaneously. The Mathematica com-
mands, which allow the determination of the five unknowns,
are the following:
P = 0.001;
For[i = 1,i < 17,{so11260=
Findo Ka[T] (0.5(1- X)) 0.5(1.5(1- X))
A1F.5 4,[1] 1.5 4,[2] A 0.5==
X(2- X)/P 4,[3],y[3](2- X)
== X,y[2] (2- X) == 0.5 (1 - X),
y[l](2- X) == 1.5(1- X),
Z, 3 + (-1+ B,) Z A 2+ Z, (A - 3BA 2 - 2B,)
-A, Bv+ B, 2+BA 3== 0}
{Zv, 0.95},MaxIterations - 10000],
Pp[i] = P,Kv[i] = v,[3] / ([1] A 1.5 v,[2] A 0.5)/
.so11260,X[i]= X /.so11260,
Print[i," ",P],P = P + 1000,i= i +1}]

In Figure 2, we plot Kv vs. pressure at various temperatures
(900, 1080, 1260 and 1440 �R). Values ofKK are significantly
different from unity, indicating that this factor must be taken
into account at high pressures. The extent of reaction at equi-
librium vs. pressure, for the same temperatures, is represented
in Figure 3. The extent of reaction approaches unity at high

We have shown through simple examples how one can
take advantage of the numerical and graphical capabilities
of Mathematica to perform flash calculations for multi-com-
ponent mixtures, solubility determination, and high-pressure
chemical equilibrium computations. Similar computations can
be performed very easily using Matlab. 7

x: liquid mole fraction
y: vapor mole fraction

Po,: critical pressure (psia)
To.: critical temperature (�R)
Pr,: reduced pressure
Tr,: reduced temperature
m: acentric factor
Z: compressibility factor

4, 4,: fugacity coefficients

K: equilibrium constant

Pressure in psia
Figure 2. K vs. pressure at various temperatures.





900� R
0.2 \ 1080� R
1260� R
1440� R
0 2000 4000 6000 8000 10000 12000 14000

Pressure in psia
Figure 3. Extent of reaction vs. pressure
at various temperatures.
Chemical Engineering Education

0 2000 4000 6000 8000 10000 12000 14000









0vap :

departure from ideality enthalpy (cal/mol)

binary interaction parameter

number of components

feed mole fraction

universal gas constant (cal/mol.�R)

activity of species i (psia)

stoichiometric coefficient

liquid phase fugacity coefficient

vapor phase fugacity coefficient

fraction vaporized = V/F

1. results= 1; search_person_id= 1536>
2. Tester, J.W., and M. Modell, Thermodynamics and Its Applications,
3rd Ed., Prentice Hall, Upper Saddle River, NJ (1996)
3. Prausnitz, J.M., R.N. Lichtenthaler, and E.G. deAzevedo, Molecular
Thermodynamics of Fluid-Phase Equilibria, 3rd Ed., Prentice-Hall,
Englewood Cliffs, NJ (1998)
4. Sandler, S.I., ( /.... .i....t i ........ ring Thermodynamics, Wiley, NY,
3rd Ed. (1999)

6. Henley, E.L., and J.D. Seader, Equilibrium-Stage Separation Opera-
tions in Chemical Engineering, Wiley, NY (1981)
7. do?objectType= author&objectId=1093893> 7

Vol. 42, No. 1, Winter 2008

This one-page column will present practical teaching tips in sufficient detail that ChE educators can
adopt the tip. The focus should be on the teaching method, not content. With no tables or figures
the column should be approximately 450 words. If graphics are included, the length needs to be
reduced. Tips that are too long will be edited to fit on one page. Please submit a Word file to Phil
Wankat , subject: CEE Teaching Tip.


The University ofAlabama

Every semester I survey my students to learn why they
chose to attend UA, why they chose to major in chemi-
cal engineering, and what kind of career they aspire
to. I also ask about their major interests and hobbies. I found
that chemical engineering students organized COE Does
ART, an official campus organization in which engineering
students put on a musical each semester. Several chemical
engineering students competed in a recent university "battle
of the bands" as members of the rock group Liquid Metal.
Others are members of the university's marching band, sing
in a musical, write short stories, create art, or make crafts.
Engineering is an incredibly creative profession so it should
not surprise anyone that their activities outside the classroom
reflect very creative personalities. Can we harness this creativ-
ity for the benefit of chemical engineering education?
Many years ago, CEE published the poem "The Reynolds'
Number Song"'11 that humorously celebrates the Reynolds
Number. Since most current students enjoy rap music, I
found it useful to modify the
poem slightly and perform it
in class as the "Reynolds Rap"
(pun intended). Imagine the
class pounding out a rhythm on
their desks and joining me in the
chorus, "Take a 'd' times a 'v
and a rho by a mu, put them all
together with a little bit of glue,
now you got a number gonna see
you through, and tell you what
the [insert colorful language of
your choice] fluid's gonna do."
It is certainly hard to forget the
equation after that exercise!

� Copyright ChE Division ofASEE 2008

In the fall of 2006, I recorded the "Reynolds Rap" in my
home studio and played it for my 90 freshmen in the Intro-
duction to Chemical Engineering seminar. During the 2007
summer lab several students asked me to write a song about
distillation columns, suggesting the title "Steppin' Off the
Stages." Since our end-of-summer-lab picnic was themed The
Fiesta before the Siesta, I set it to a Latin rhythm and played
it at the picnic. The students very much enjoyed both songs,
and many downloaded them.
These songs are available free as MP3 downloads at
for your listen-
ing pleasure and-if you dare-to use in your classroom. I
believe there are many would-be rock stars in the profession
that have also written or will write a chemical engineering-
themed song. Make a good recording and send me the MP3
file and lyrics by e-mail (alane@eng.ua.edu), and you just
might find it posted on this Web site and played on chemical
engineers' iPods around the globe!

1. Harriott, P., "The Reynolds' Num-
berSong,"( i ..... in . 8(1),12
(1979) 1

Figure 1. Doobie 'Doghouse' Wilson (AKA Alan Lane)
with his dobro.
Chemical Engineering Education

Alan M. Lane is a professor of
chemical and biological engineer-
ing at The University of Alabama.
He conducts research in hetero-
geneous catalysis and teaches
chemical reaction engineering.
At night he performs as Doobie
"Doghouse" Wilson in local bars
and coffee houses housewilson.com>.

Me,= Rteaching tips )

M]n1 class and home problems

The object of this column is to enhance our readers' collections of interesting and novel prob-
lems in chemical engineering. Problems of the type that can be used to motivate the student by
presenting a particular principle in class, or in a new light, or that can be assigned as a novel home
problem, are requested, as well as those that are more traditional in nature and that elucidate dif-
ficult concepts. Manuscripts should not exceed 14 double-spaced pages and should be accompanied
by the originals of any figures or photographs. Please submit them to Professor James O. Wilkes
(e-mail: wilkes@umich.edu), Chemical Engineering Department, University of Michigan, Ann
Arbor, MI 48109-2136.


An Exercise in Validation

of Computational Results

Ben-Gurion University of the Negev * Beer- I,. .., 84105, Israel
Tel-Aviv University * Tel-Aviv 69978, Israel
Auburn University * Auburn, Alabama 36849
University of Connecticut * Storrs, CT 06269
M mathematical software packages such as Polymath, Mordechai Shacham is the Benjamin H. Swig professor of chemical
MATLAB and Mathcad are widely used for engi- processes and former chair at Ben-Gurion University of the Negev in
neering problem solving. The application of several Israel. He received his B.Sc. and D.Sc. degrees from the Technion, Israel
Institute of Technology. His research interests include analysis, modeling
packages to typical chemical engineering problems has been and regression of data, applied numerical methods and prediction, and
demonstrated by Cutlip, et al.[' The main characteristic of consistency analysis of physical properties. He is past president and an
honorary fellow of the Israel Institute of Chemical Engineers, and the recipi-
these packages is that they provide a "problem solving envi- ent of the 2000 CACHE Award for excellence in computing in chemical
ronment (PSE)"[2] rather than just scientific subroutine librar- engineering education.
ies callable from programming languages that were popular Neima Brauner is a professor at Tel-Aviv University in TelAviv, Israel. She
a few years ago. For routine use of the software packages in received her B.Sc. and M.Sc. in chemical engineering from the Technion
problem solving, there is no need to be an expert in either Institute of Technology in Haifa, Israel, and her Ph.D. in mechanical
problem solving, there no need to be an expert in either engineering from the Tel-Aviv University. Her research interests include
numerical methods or programming. In most cases, it is suf- hydrodynamics and transport phenomena in two-phase flow systems
ficient to specify the mathematical model of the problem using and development of interactive statistical and numericalmethods for data
analysis in process analysis and design. She is an editor for Reviews in
the syntax required by the package. The technical details of Chemical Engineering, and associate editor of Heat Transfer Engineering.
the solution are carried out by the package. There are cases,
SW. Robert Ashurst is currently an assistant professor of chemical
however, where the solution is obviously incorrect and/or the engineering at Auburn University. His research focuses on design of
program stops with an error message. molecular precursors for advanced monolayer films, tribology at the
micro- and nano-scale, and the influence of surface chemical treatments
The failure of a software package to reach a correct solu- on micro- and nano-scale devices. He completed his Ph.D. in chemical
tion is commonly due to the presence of errors in the math- engineering at the University of California at Berkeley in 2003 under the
support of a National Science Foundation Graduate Fellowship. He is an
ematical model or the selection of an inappropriate solution active member of AIChE.
algorithm and/or inappropriate algorithm parameters (such Michael B. Cutlip has spent his entire academic career at the University of
as initial estimates, error tolerances, etc). Validation of the Connecticut, where he served as department head of chemical engineer-
mathematical models has been discussed in detail by Brauner, ing and director of the university's honors program. He is a former chair
and national program chair for the Chemical Engineering Division, and he
et al.[3] In this paper, we present an exercise that can be used co-chaired the ASEE Summer School for chemical engineering faculty in
Copyright ChEDivi2002. He is also co-author of the POLYMATH software package.
SCopyright ChE Diision of ASEE 2008
Vol. 42, No. 1, Winter 2008 5j

in educating students to identify difficulties associated with
an algorithm for the solution of a problem and help them to
select appropriate parameters for the algorithm.
This exercise involves numerical solution for a system of
ordinary differential equations (ODE's). It can be used at
two levels.
1. The first is in courses where students routinely use
numerical software for solving ODE's (i.e., Reaction
Engineering, Process Control and Process Simula-
tion). In such courses, this part of the exercise can be
used to demonstrate that computer solutions can be
inaccurate or even completely incorrect, and therefore
it is always necessary to validate a numerical solu-
2. At the more advanced level, the second part of the
exercise can be used in courses related to mathemati-
cal modeling and numerical analysis. In such courses,
this part of the exercise provides training and demon-
stration of advanced topics such as error estimation
and step size control, error ;1,. *. ,, i.. . and stiffness
of ODE's.
In this work, the mathematical software packages Poly-
math 6.11 and MATLAB2 are used for the demonstration,
but other packages may be used in a similar manner.

The problem considered in this exercise involves the
dynamic simulation of a system of three interacting tanks,
as depicted in Figure 1. The model equations (as were
provided Ashurst and Edwards'41) that can be derived from
mass balances on each tank are
-- -h D -D,2 dhl
q, -q, = 2 +h, (1)
q q 4 H2 dt (1)

q2 + q4 - q5 = Te2h2 (2)
q5 + q3 - q6 = h 3 (3)

where q represents the in/out volumetric flow rate (m3/s),
H the tank height (m), h the liquid height (m), D1,2 the
upper/lower diameter ratio of tank 1 (m), and t is the time
(s). The volumetric flow rates out of the tanks, through a
drainage area, Ao (m2), are given by

q4 = Ao 2g(h- h2) ifh, > h2

q4 = -Ao2g(h2 - hi) otherwise (4)

q5 = Ao 2g(h -h3) ifh2 > h3

q, = -Ao2g(h3 - h) otherwise (5)

q6 = Ao2gh (6)

The response of the system to a step change in the inlet flow
rates is to be investigated. The model equations and the nu-
merical data pertinent for a particular case, as they have been
entered into the Polymath 6.1 software package, are shown in
Table 1. Note that the equations and comments in the Polymath
input provide a complete documentation of the problem (in
Polymath the comments are indicated by the sign: #).

Polymath Model for the "Train of Interactive Tanks" Problem
No. Equations and # Comments
1 # Mass balances on the tanks
2 d(hl)/d(t) = If hl >= H1 Then (0) Else (ql - q4) / ((pi / 4) *
((D2 + ((Dl - D2) / HI) * hi) A 2)) # Liquid level in Tank 1 (m)
3 d(h2)/d(t) = If h2 >= H2 Then (0) Else (q2 + q4 - q5) (pi *
exp(2 * h2) ) # Liquid level in Tank 2 (m)
4 d(h3)/d(t) = If h3 >= H3 Then (0) Else (q5 + q3 - q6) (h3 A 2)
# Liquid level in Tank 3
5 # Outlet flow rates are functions of potential flow theory.
6 q4 = If hl < h2 Then - Ao * sqrt(2 * g * (h2 - hl)) Else Ao *
sqrt(2* g * (hl -h2))#
Volumetric flow rate out of Tank 1 (m^3/s)
7 q5 = If h2 < h3 Then - Ao * sqrt(2 * g * (h3 - h2)) Else Ao *
sqrt(2 * g * (h2 - h3)) #
Volumetric flow rate out of Tank 2 (m^3/s)
8 q6 =Ao * sqrt((2 * g * h3))#
Volumetric flow rate out of Tank 3 (m^3/s)
9 # Input variables
10 ql = .05 # Volumetric flow rate into Tank 1(m^3/s)
11 q2 = .05 # Volumetric flow rate into Tank 2(mA3/s)
12 q3 = If t <= 50 Then .05 Else.05 +.05* .1 #
Volumetric flow rate into Tank 3(mA3/s)
13 # Constants
14 H1 = 4 # Height of Tank l(m)
15 H2 = 2.5 # Height of Tank 2 (m)
16 H3 = 2.5 # Height of Tank 3 (m)
17 D1 = 4 # Upper diameter of Tank l(m)
18 D2 = 1 # Lower diameter of Tank 2 (m)
19 Ao = pi / 100 # Area of drainage orifice at bottom of tanks
20 pi =3.141592654
21 g=9.81
22 # Initial Conditions
23 hl(0) = 1.809 # Liquid level in Tank 1 at steady state(m)
24 h2(0) = 1.68 # Liquid level in Tank 2 at steady state(m)
25 h3(0) = 1.163 # Liquid level in Tank 3 at steady state(m)
26 t(0) = 0
27 t(f)= 2400.

'POLYMATH is a product of Polymath Software (< , .. ........ t, .....
2MATLAB is product of The Math Works, Inc. ()

Chemical Engineering Education

The system is started at steady-state values of h, h2 and h3
(see lines 23 through 25), where the inlet flow rates to the
different tanks are ql = q2 = q3 = 0.05 m3/s. At t = 50 s, a step
change of +10% is introduced into q3. The system is simulated
up to the final time of tf = 3600 s.
The integration of the model equations using the RKF45
algorithm of Polymath with the default parameters yields the

Figure 1. A train of three interactive tanks.!4'

Change of the Variable Values after Introduction of a 10%
Step Change in Inlet Flow Rate to Tank 3
Variable Initial value Minimal Maximal Final
value value value
t 0 0 2400 2400
h, 1.807 1.807 1.854328 1.854328
h, 1.68 1.679873 1.725494 1.725494
, 1.162 1.162 1.221394 1.221394
q 0.0495908 0.0495908 0.0499622 0.0499475
q, 0.1001531 0.0971087 0.1001531 0.0988002
q 0.1500039 0.1500039 0.1537897 0.1537897
dhl/dt 9.39E-05 8.66E-06 9.39E-05 1.17E-05
dlh/dt -6.22E-06 -6.22E-06 3.06E-05 1.16E-05
dh/dt 1.11E-04 -5.39E-07 2.71E-04 7.08E-06

0 500 1000 1500
time (s)

2000 2500 3000

error message L i i. . sqrt ofn ,... ,a,.' number is not allowed,
,.'p q 2, / 3)," which indicates that the height of the liquid level
in the 3rd tank (h3) has somehow reached a negative value dur-
ing the numerical integration. A quick check can be made on
the dynamics by reducing the integration final time to 2400 s.
As there is no abnormal behavior or unexpected trends during
this numerical integration, as shown in Table 2 and Figure 2,
the assignment is to check
and verify the error in the
q3 integration up to the time
H H3 of 3600s. The options that
need to be investigated for
This problem are:
1. There are possible er-
H3 rors in the implementa-
\ h3 q tion of the Runge-Kutta
q5 q6 1 . ,, ...

2. The default truncation
error tolerance of the program is too large
for this problem, and the negative value of
h3 is a result of large truncation error and
propagated error in some of the variables.
3. The system of equations is stiff. The
explicit RKF45 may exhibit oscillatory be-
havior if the ;i,. . inr. -t, step size is ... ,..a
than that dictated by the stability limit.

1. Most software packages include some
variation of the Runge-Kutta algorithm
- with truncation error estimation and step
size control - as one of the ODE solver tools.
Using the same or a similar algorithm in a
different software package can provide ad-

ditional information regarding the source of the error.
2. The appropriateness of the truncation error tolerance
can be tested by rerunning the problem using various error
tolerances and comparing the variable values obtained at a
specified time (say tf). It is preferable to check the values at
the final time, as at this point the error propagated throughout
the integration interval is considered.
3. The stiffness of the problem can be determined by defin-
ing the f vector as

f= dh (7)

Figure 2. Change of the liquid level in the third tank
after introduction of a 10% step change
in its inlet flow rate (Polymath Solution).
Vol. 42, No. 1, Winter 2008

where dhl/dt, dh2/dt, and dh3/dt are given by Eqs. (1-3),
respectively. The stiffness of the system can be assessed by
calculating the stiffness ratio of the df/ah (Jacobian) matrix.
The stiffness ratio (SR) is defined by (see, for example, p.
245 in Rice and Do[5])
max IX,
SR -- (8)
min X,

where max X, is a negative eigenvalue of the 8f/ah matrix
with maximal absolute value, and min JI, is a negative ei-
genvalue of the 8f/ah matrix with minimal absolute value. If
SR>1000, the system of ODE's is considered a stiff system. If
the system is stiff, using the appropriate algorithms (i.e., stiff
or stiffbs in Polymath) should prevent obtaining the negative
h, value during integration.


Polymath 6.1 can be used to convert the model of the
problem into a MATLAB function. The Polymath generated
function is shown in lines 21 through 42 of Table 3. Note that
Polymath reorders the equations and changes the syntax of
the model according to the requirements of MATLAB. Some
manual editing of the Polymath generated function was done
in order to enable more concise presentation (the comments
were put in the same line with the commands and the if-else
statements were put in a single line).
The "main program" that runs the Polymath generated
function (shown in lines 1 through 19 of Table 3) is avail-
able as a template in the "Help" section of Polymath. Only


1.2 -

' 1.18

2 1.17


0 500 1000 1500 2000 2500 3000 3500
Time (s)

Figure 3. Change of the liquid level in the third tank
after introducing a 10% step change in its
inlet flow rate (ode45 Algorithm, RelTol = 10-3).

the function name has to be added (see line 1) and the initial
values of the variables had to be copied from the Polymath
generated model (see lines 3 and 4). The ode45 library func-
tion of MATLAB can be used to solve the problem. This
function is based on essentially the same algorithm as the
RKF45 routine of Polymath.
MATLAB solution of the program in Table 3 yields a simi-
lar solution to that obtained with Polymath, except that the
liquid level in the 3rd tank (h,) exhibits oscillatory behavior
(see Figure 3). These oscillations indicate that there may be
difficulties with the MATLAB solution as well, and suggest
that the error tolerances should be examined.

In the Polymath RKF45 algorithm, the integration step size
is controlled by two parameters: the relative error tolerance
(RelTol) and the minimum number of reporting points (RP).
The default values of these parameters can be found under the
Polymath icon of "Setup preferences and numerical param-
eters," which is under Ordinary Differential Equations and
RKF45: RelTol = 106 and RP = 100. The RKF45 program
includes an adaptive step size control algorithm that monitors
the estimate of the integration error and adjusts the step size of
the integration in order to keep the error below the specified
threshold value. The step size is increased only as long as at
least RP full steps are carried out inside the integration inter-
val. Therefore, a solution of a lower accuracy can be obtained
if the error tolerance is increased, and, at the same time, the
number of reporting points is reduced. For instance, if the
error tolerance is increased to RelErr = 104 and the number
of reporting points is reduced to RP = 30, difficulties
are encountered earlier. At some time t < 2400, the
error message ("negative square root") is issued by
the program. With adjustment of RelErr = 1010 while
holding RP = 30, however, no error message is issued
when the integration is carried out up to tf = 3600. This
implies that a correct solution has been obtained. In
the MATLAB ode45 function, the integration step size
is adjusted so as to keep the estimated error below the
value of e , which is a function of the specified relative
error tolerance, RelTol, and the specified absolute error
tolerance, AbsTol.

, = max(|h, RelTol, AbsTol) (9)

The default values of the ode45 function for the
absolute error tolerance and for the relative error toler-
ance are AbsTol = 106 and RelTol = 10, respectively.
Since the values of h1, h2, and h3 are of the order of
magnitude of 1 (one), the AbsTol is inactive and the
e values are of the order of 103. With the default er-
ror tolerances in MATLAB, the oscillatory solution
(shown in Figure 3) is obtained. Reducing the rela-

Chemical Engineering Education


MATLAB Program for the "Train of Interactive Tanks" Exercise
No. Commands and % Comments
1 function TankTrain
2 clear, clc, format short g, format compact
3 tspan = [0 3600]; % Range for the independent variable

4 yO = [1.807; 1.68; 1.162]; % Initial values for the dependent variables
5 disp(' Variable values at the initial point );
6 disp([' t = ' num2str(tspan(l))]);
7 disp(' y dy/dt ')
8 disp([y0 ODEfun(tspan(1),y0)]);
9 [t,y]=ode45(@ODEfun,tspan,y0);
10 for i= :size(y,2)
11 disp([' Solution for dependent variable y' int2str(i)]);
12 disp([' t y' int2str(i)]);
14 plot(t,y(:,i));
15 title([' Plot of dependent variable y' int2str(i)]);
16 xlabel(' Independent variable (t)');
17 ylabel([' Dependent variable y' int2str(i)]);
18 pause
19 end
20 %------- ----- ---
21 function dYfuncvecdt = ODEfun(t,Yfuncvec);
22 hl =Yfuncvec(l);
23 h2 =Yfuncvec(2);
24 h3 =Yfuncvec(3);
25 Ao = pi / 100; %Area of drainage orifice at bottom of tanks (m^2)
26 g = 9.81;
27 ql = .05; %Volumetric flow rate into Tank i(m^3/s)
28 q2 = .05; %Volumetric flow rate into Tank 2(m^3/s)
29 if (t <= 50), q3 = .05; else q3 = .05 + .05 * .1; end %Volumetric flow rate into Tank 3(m^3/s)
30 H1 = 4; %Height of Tank l(m)
31 H2 = 2.5; %Height of Tank 2 (m)
32 H3 = 2.5; %Height of Tank 3 (m)
33 D1 = 4; %Upper diameter of Tank l(m)
34 D2 = 1; %Lower diameter of Tank 2 (m)
35 if (hl < h2), q4 = 0 - (Ao * sqrt(2 * g * (h2- hl))); else q4 =Ao * sqrt(2 * g * (h - h2));end %Volumetric flow rate out of Tank
1 (m^3/s)
36 if (h2 < h3), q5 = 0 - (Ao * sqrt(2 * g * (h3 - h2))); else q5 = Ao * sqrt(2 * g * (h2 - h3)); end %Volumetric flow rate out of Tank
2 (m^3/s)
37 q6 = Ao * sqrt(2 * g * h3); %Volumetric flow rate out of Tank 3 (m^3/s)
38 if (h >= H1), dhldt = 0; else dhldt = (ql - q4) / (pi 4 * (D2 + (D1 - D2) / H1 * hi) A 2); end %Liquid level in Tank 1 (m)
39 if (h2 >= H2), dh2dt = 0; else dh2dt =(q2 + q4 - q5) (pi * exp(2 * h2) ); end %Liquid level in Tank 2 (m)
40 %Liquid level in Tank 3
41 if (h3 >= H3), dh3dt = 0; else dh3dt =(q5 + q3 - q6) (h3 A 2); end
42 dYfuncvecdt = [dhldt; dh2dt; dh3dt];

Vol. 42, No. 1, Winter 2008

tive error tolerance to RelTol = 106 yields a more accurate,
non-oscillatory solution. Since the latter solution is more ac-
curate by three orders of magnitude, the difference between
the solutions can serve as a measure of the "error" in the
less accurate solution. This "error" in h3 is shown in Figure
4 (note that a few points of larger errors, in the region where
the liquid level changes more rapidly, are not shown). Note
that the maximal errors are of the order of 103 (while h311.2,
see Table 2). Thus, the integration algorithm indeed keeps the
error below the desired relative error tolerance.
The error in the slope dh3/dt is shown in Figure 5 (a few
points of larger errors are not shown in this plot). The maximal
absolute errors are of the order 104, but since the value of
dh3/dt ~ 4x 106, the relative error is of the order of 25. Thus, a
reasonably low relative error of 103 in the variable (h) value
yields absurd slope values, with up to 2500% error.
It is interesting to investigate the reasons for these huge dif-
ferences in the relative error values. In Table 4, the values of
some of the variables are presented at t = 3600 as calculated
by the ode45 function with the two error tolerances. Observe
that there are four correct decimal digits in the values of h1,
q +q3, and q6 when the higher error tolerance of 103 is used.
Introducing these numbers into Eq. (3) to calculate dh3/dt


requires subtraction of two numbers that are very close. This
operation causes the relative error to increase by seven orders
of magnitude!
We conclude that the erroneous oscillations in the ode45
solution and the "negative square-root" error message of the
RKF45 algorithms are caused by the lack of consideration
of the slope value errors by the step-size control logic of the
two algorithms.


The cause of the erroneous oscillations and error mes-
sages was identified in the previous section. To complete the
exercise, however, the stiffness of the problem must also be
There are several possible ways to find the elements of the
af/ah matrix and to calculate the eigenvalues of this matrix
at various times. One option is to use the symbolic manipu-
lation options of MATLAB. The program for derivation of
the elements of the Jacobian matrix and calculation of the
eigenvalues at t = 3600 s is shown in Table 5. First, the vec-
tor of functions (Eq. 7) is defined using symbolic variables
(lines 2 through 4). Note that

Comparison of Some Variables Values as calculated
By the ode45 Function With Two Error Tolerances
Relative t(s) h3 (m) q5+q3 (m3/s) q6 (m3/s) dh3/i
106 3600 1.2003444895 0.1524588034 0.15245874 4.41426
103 3600 1.2000989262 0.1524832835 0.15244314 2.78699
% Error 0.020 -0.016 0.010 63036


0.001 -




-0.001 -

-0.0015 -
0 500 1000

Time (s)

the expressions for q4, q ,
and q6 (as functions of hl,
h2, and h3) were introduced
dt into the functions assuming
that hl e-08 time period. The symbolic
-5 differentiation is carried out
in lines 5 through 9, and
partial derivative expressions
are stored in the matrix df.
Numerical values are introduced into
the constants (lines 10 and 11) and the
final values of h, h , and h,, from Table
2 are introduced into these variables in
line 12. Finally, the numerical values
of the terms of the Jacobian matrix
(line 13) and the eigenvalues of this
matrix (line 15) are calculated. The
eigenvalues are: ?h = 0.37297e-4; X2 -
0.04478; and 3 = -0.106922. Introduc-
ing these values into the stiffness-ratio
equation yields:

SR = = 286.67
As the stiffness ratio is less than
1000, the system is not considered to
be stiff according to the stiffness-ratio

Chemical Engineering Education

Figure 4. Error in the liquid level h, (ode45 Algorithm, RelTol = 10o-).


The exercise presented here demonstrates a systematic
approach for identification of causes for a failure in an ODE
integration algorithm used to reach a correct solution of
a particular problem. The procedure starts with checking
whether a similar algorithm in a different package yields the
same solution, and whether the default algorithm parameters
are appropriate. Then, the stiffness of the problem is assessed,
and a comparison of solutions obtained by stiff and non-stiff
algorithms is carried out.
In courses where students
routinely use numerical software
for solving ODE's, carrying out
the two parts of the exercise can MATLAB Prog
be very beneficial. They demon- No. Commands
state that results of ODE solvers 1 syms ql q2 q3
must be verified. The reduction 2 f(l)= (ql- (Ao
of error tolerances and problem 3 f(2)= (q2 + (Ao
solution by different packages 4 (
4 f(3)= ((Ao * sq:
are very efficient in detecting in-
accurate solutions. The Polymath 5 for i=1:3
package, in particular, is well- 6 df(i,1)=diff(]
suited for such tests as it enables 7 df(i,2)=diff(
exporting the model equation to 8 df(i,3)=diff(
other packages. 9 end
For students with appropriate 10 Ao = pi /100;
numerical analysis background, H=4;
11 H1=4; H2 = 2
this two-part exercise provides
a good demonstration of the 12 hl=1.8656; h2=
concepts of error estimation, step 13 dfnumeric=sul
size control, error propagation 14 format long g
and stiffness in systems of ODE's 15 Lambda=eig(df
using a real-life example. This
exercise demonstrates that a complete 2.00E-04
analysis and validation of the results
are important in problem solving and 1.50E-04
for maintaining confidence in a particu-
lar software package. While students 1.00E-04
may use software packages for routine *
problem solving with a limited knowl- 5.00E-05
edge of the numerical methods utilized *
by the package, there is also a need to 0.OOE+00
include numerical methods within the *
curriculum. This knowledge is needed 2 -5.00E-05
when difficulties are encountered during
numerical problem solving. -1.00E-04

The authors wish to thank Auburn -2.00E-04
University chemical engineering stu- 0
dent Robert Edwards for first bringing
this problem to our attention.

Vol. 42, No. 1, Winter 2008

1. Cutlip, M.B., J.J. Hwalek, H.E. Nuttall, M. Shacham, J. Brule, J. Wid-
man, T. Han, B. Finlayson, E.M. Rosen, and R. Taylor, "A Collection of
Ten Numerical Problems in Chemical Engineering Solved by Various
Mathematical Software Packages," Comput. Appl. Eng. Educ., 6(3),
2. Enright, W.H., "The Design and Implementation of Usable ODE
Software," Numerical Algorithms, 31, 125 (2002)
3. Brauner, N., M. Shacham, and M.B. Cutlip, "Computational Results:
How Reliable Are They? A Systematic Approach to Modal Validation,"
Chem. Eng. Ed., 30(1), 20 (1996)
4. Ashurst, WR., and R.J. Edwards, personal communications (2006)
5. Rice, R.G., and D.D. Do, Applied Mathematics and Modeling for
Chemical Engineers, John Wiley (1995) 7

gram for Calculating the Eigenvalues of the Of/Oh Matrix

g hi h2 h3 HI Ao D1 D2 g
* sqrt(2 * g * (hi - h2)))) / ((pi 4) * ((D2 + ((Dl - D2) HI) * h) A 2));
* sqrt(2 * g * (hi - h2)))- (Ao * sqrt(2 * g * (h2 - h3)))) / (pi * exp(2 * h2));
rt(2 * g * (h2 - h3))) + q3 - ((Ao * sqrt((2 * g * h3))))) / (h3 A 2);


g =9.81; q = 0.05; q2 = 0.05; q3=0.05 + 0.05 * 0.1;
.5; H3=2.5; D =4; D2 =1;
1.7367; h3=1.2282;


ure . rror n e qu eve s ope: (o e gor m, e o = .

500 1000 1500 2000 2500 3000
Time (s)
ht l;i ;il Irl l1 dh /dt A d+ (- Ar ih RlTr l 1 n-3)

University of Florida Home Page
© 2004 - 2011 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Powered by SobekCM