chemical engineering education
Visit
us
on the
Web
at
Jennifer
Sinclair
Curtis
... 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 PengRobinson Equation of State Using Mathematica (p. 47)
Binous
Random Thoughts: Student Ratings of Teaching: Myths, Facts, and Good Practices (p. 33)
Felder, Brent
Editorial: Engineers Deserve A Liberal Education (inside front cover)
Kiilg
Celebrating ChE In Song (p. 52)
Lane
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
V
Tufts University
AIL
40
w
editorial
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 educationa 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
careersto 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
cultures.
* 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
inundergraduateengineeringeducation.Buttheproblem
is that a nominally fouryear 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 premedicaleducation mayactuallyincreasethe
number and diversity of engineering students.
Some preengineering 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 entrylevel 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.Surelyengineeringthe problemsolv
ing professioncan find the way to make the change.
My class at Yale had its 50th reunion this year. In the
reflective 50yearclass 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
LABORATORY
Feature
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 largescale 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 laboratorytested 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 studentlearning opportunities (e.g., treatment of reac
tionrate 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
material.
* 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
References.
EDITORIAL AND BUSINESS ADDRESS:
Chemical Engineering Education
Department of Chemical Engineering
University of Florida * Gainesville, FL 32611
PHONE and FAX: 3523920861
email: cee@che.ufl.edu
EDITOR
Tim Anderson
ASSOCIATE EDITOR
Phillip C. Wankat
MANAGING EDITOR
Lynn Heasley
PROBLEM EDITOR
James O. Wilkes, U. Michigan
LEARNING IN INDUSTRY EDITOR
William J. Koros, Georgia Institute of Technology
PUBLICATIONS BOARD
* CHAIRMAN
John P. O'Connell
University of Virginia
* VICE CHAIRMAN *
C. Stewart Slater
Rowan University
* MEMBERS
KristiAnseth
University of Colorado
Jennifer Curtis
University of Florida
Rob Davis
University of Colorado
Pablo Debenedetti
Princeton University
Diane Dorland
Rowan
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
> DEPARTMENT
10 Chemical Engineering at Tufts University
Nakho Sung, and Daniel Ryder
> EDUCATOR
2 Jennifer Sinclair Curtis of the University of Florida
Nicholas Peppas, Gintaras Reklaitis, and Erik Ydstie
> RANDOM THOUGHTS
33 Student Ratings of Teaching: Myths, Facts, and Good Practices
Richard M. Felder and Rebecca Brent
> CLASSROOM
23 Active Problem Solving and Applied Research Methods in a Graduate
Course on Numerical Methods
Eric L. Maase and Karen A. High
> CURRICULUM
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
> LABORATORY
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
> TEACHING TIPS
52 Celebrating ChE In Song
Alan M. Lane
> CLASS AND HOME PROBLEMS
47 Applications of the PengRobinson 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 00092479) 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 326116005. 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 326116005. 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
NICHOLAS PEPPAS
The University of Texas
GINTARAS REKLAITIS
Purdue University
ERIK YDSTIE
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 Floridaas 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.
THE EARLY YEARS
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 postDepression 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 dedicationbordering on stubborn
nessto 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 19yearold with golden locks,
holding an alreadyfilledout 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 MidWestern 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 graduateher 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 twoyear 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 nonFickian 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
Florida.
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 pseudoconvective 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,404409 (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 Pennsylvaniawhere 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.
GRADUATE SCHOOL AND FAMILY
Jennifer started at Princeton University in 1983, while
Gavin was employed a state awayat 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 publishedamong othersin a seminal
issue of the AIChE Journal in 1989 ["GasParticle Flow in
a Vertical Pipe with ParticleParticle Interactions," AIChE J.
(1989) 35, 14731486]. 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 SloanKettering, 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
4
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."
AN ACADEMIC CAREER IN THE SERVICE OF
STUDENTS
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 "preexisting 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 locationPittsburghwhere 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 homeJennifer, 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 CARNEGIE MELLON
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 firstrate 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
Grossmann.
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
concernedGavin 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
ingdecision 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 mid1990s, 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 Arizonawhose
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 preexisting condition. And
that's when Gavin and Jennifer remembered their roots.
THE PURDUE YEARS
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 dualacademic
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 cotaughi 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 cando style, reenergizing the staff, recruiting new fac
ulty, and rallying the firstyear 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.
` FORGING ON
Jennifer's optimistic outlook prevailed. In spite of her many
administrative duties and the neverending 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 challengesto 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.
AN INTERNATIONALLY KNOWN RESEARCHER
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 gassolid 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
softwarepackage vendors. Both companiesFluent and
CFXhave adapted her models for dilute and densephase
gassolid 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 velocimetera purchase that helped her propel her
early research program further toward the leading edge of
science.
GRADUATE STUDENTS WHO DID THEIR THESES WITH JENNIFER CURTIS
1. John Doney ChE, M.S., Carnegie Mellon Uni
versity, "Heat Transfer in Dilute ParticleLaden
Turbulent Flows," 1993
2. Eduardo Bolio, ChE, Ph.D., Carnegie Mel
lon University, "Dilute Turbulent GasSolid
Flows with ParticleParticle Interactions,"
1994
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 GasSolid 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 GasSolid 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 GasSolid 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 BlakePowell, ChE, M.S., Purdue
University, "Comparison ofAdHoc 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 GasSolid Flow
in Risers," 2000
12. Edward (Nick) Jones, ChE, Ph.D., Purdue
University, "The Effect of PSD on GasSolid
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 ParticleLaden 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
souriRolla)
16. Kunn Hadinoto, ChE, Ph.D., Purdue Uni
versity, "DensePhase 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 (coadvised student)
18. Bill Ketterhagen, ChE, Ph.D., Purdue
University, "DEM Simulations of Granular
Discharge from Hoppers," 2006 (coadvised
student)
19. Robert Hamilton, ChE, Ph.D., Purdue Uni
versity, "Evolving Particle Size Distribution
in the Flow of GasSolid Mixtures," (1999 )
(coadvised 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 LiquidSolid
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  ) (coadvised 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
multiyear 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 EnergyOffice 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 SocietyPetroleum
Research Fund (ACSPRF), as a member of an NSForga
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 cochair, 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.
AN EDUCATOR, ADVISOR, AND MENTOR
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 universitywide Teaching for
Tomorrow Program. In 2001, Jennifer served, along with
then Purdue Engineering Dean Linda Katehi, on a 10mem
ber invited NSF panel on the Future of Engineering Educa
tion. From 20032006,
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
twoweek 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 Henthornnow an assistant professor in chemical
engineering at the University of MissouriRollaremarks,
"I think the thing that impressed me most about Jennifer
was her leadership styleshe 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."
A NEW LIFE AT THE UNIVERSITY OF FLORIDA
In 2001 Jennifer met Barry Curtis, another Purdue chemical
engineer, who had graduated two years before her (in 1981), hav
ing done the fiveyear coop 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
tripto 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 Managementthe same school from which her father graduated
24 years earlierand 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.
BY NAKHO SUNG AND DANIEL RYDER
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. Fulltime 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 hightech 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
curricula.
For the past 15 years, the department has been assembling
a core group of biological engineering faculty. Graduate
research is divided within three main areasbiological 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.
TUFTS SCIENCE AND TECHNOLOGY CENTER
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, ElectroOptic 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.
CHEMICAL AND BIOLOGICAL ENGINEERING
The department has a longstanding history relating to the
development of educational programs in bioengineering. In
1986, the then Chemical Engineering Department initiated
a Massachusettssponsored continuing education program
to train technical staff for biotechnology business. The Bio
technology Engineering Certificate program, which began
as a series of handson, introductory, technical courses and
remains a vital link to local industry, has led to the develop
ment of a graduatedegree program. Today, the department
offers graduate degrees in both chemical engineering and
biotechnology engineering, including courseonly M. Eng.
degrees attractive to working professionals, as well as tradi
tional thesisbased M.S. and Ph.D. programs.
In 2004 the department hosted a series of NSFsponsored
workshops focused on the development of chemical and bio
logical engineering curricula. Participants from academia and
local bioindustries 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
upperlevel course offerings have been supplemented with
applied biologybased electives (see Table 2). We are currently
redesigning the senior laboratory courses to comprise indus
trysponsored research projects in biotechnology, metabolic
engineering, energy sustainability, systems engineering, and
advanced materials.
FACULTY
Current Faculty
The current department faculty includes nine tenuretracked
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
FlytzaniStepha
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 subppm levels
via reversible adsorption. FlytzaniStephanopoulos'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 initiativethe 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; plantwide simulations, scaleup, 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 structurefunction relationships,
including those of fibrous proteins (collagens, silks) and
polysaccharides. They apply physiological and genetic engi
12
:hBE faculty. Top, LR: Christos Georgakis, Blaine Pfeifer,
bum Lee, Jerry Meldon, Daniel Ryder, Hyunmin Yi.
LR: Ramin Haghgooie, Nakho Sung, Eric Anderson.
neering approaches to control biopolymer chemistry in order
to then explore selfassembly 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
handson experience when the
department's distillation column
needed repair. Below, Danny
Pierre, a graduate student in Pro
fessor FlytzaniStephanopoulos'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
obesityrelated 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 wellcontrolled 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
microfluidics and stem cell biology to generate flowthrough
incubators that afford spatial control over the chemical en
vironment of the cultured cells. This incubator is used for
a number of ongoing studies involving heterogeneous cell
components (e.g., adipocytepreadipocyte coculture) 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 bioprocesses 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 FlytzaniStephanopoulos's catalysis lab. Right, graduate student Lisa
Schupmann works with a multireactor system in Professor Georgakis's Systems Research Institute.
transplant genetic material responsible for an important thera
peutic product into a convenient and processfriendly bacterial
microorganism (such as Escherichia coli) for eventual product
scaleup 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: Solgel 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 IIIV nitride
semiconductors.
Professor and Chair Nakho Sung's research in the area
of structurepropertyprocessing relationships of polymers
and composites includes development of in situ monitoring
techniques for characterization of reactions in polymers and
composites via fiberoptic 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 distalend
and evanescent fiber optic probes. UV reflection spectroscopy,
coupled with bifurcated fiberoptic 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 longterm stability against tempera
ture and moistureinduced 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
lowcost membranes for fuel cell application and hydrogen
purification. Current focus is on anionic membranes suitable
for lowtemperature 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 hybridizationbased
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 channelcan 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 highthroughput
biosensors and BioMEMS devices for environmental and
biological threat detection. In closely related collaboration
with FlytzaniStephanopoulos, they are also developing na
Chemical Engineering Education
A closeup of the Tufts Chemical and Biological Engine
Xray diffractometer.
noscale gold particlebased 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 highthroughput biosensing platforms and implant
able biophotonic devices for biomedical and environmental
applications. They exploit the stimuliresponsive 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 thermoresponsive morphology transition of common
biopolymers such as gelatin and agarose can lead to efficient
means for consistent manufacturing of nanometerscale sur
face diffraction gratings under mild processing conditions.
They are currently working on building allbiopolymeric
implantable waveguides, photonic bandgap crystal fiberbased
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 "FibraCel"
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
Fellowship.
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 yearsincluding
10 years as chairis a testament to his
dedication to the profession.
TABLE 1
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
ACKNOWLEDGMENTS
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
TABLE 2
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
CHEMICAL KINETICS, HEAT TRANSFER,
AND SENSOR DYNAMICS REVISITED
in a Simple Experiment
MARIA. E. SAD, MARIO R. SAD, ALBERTO A. CASTRO, TERESITA F. GARETTO
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
resources).
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 semiindustrial equipment to computational simulations
of real processes.
CONTENTS: THEIR RELEVANCE
AND DESCRIPTION
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 noninstanta
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.
THEORETICAL FUNDAMENTALS
AND MATHEMATICAL MODEL
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
model.
This set of differential equations can be solved using an
appropriate numerical method (Euler, RungeKutta, 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,).(Tt)= 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 firstorder 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
EXPERIMENTAL
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 coldstorage 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 (1020 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.
C EXPERIMENT #1
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).
370
360
350
2
. 340
3
 330
I
320
310
300
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 wellmarked 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]
( EXPERIMENT #2
Natural cooling of the vessel with no chemical
reaction
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 353363 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.
( EXPERIMENT #3
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 (343353 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 25 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 temperaturetime experimental
values, by way of numerical integra
tion of the differential equations (with
methods adequate for "stiff' systems)
using any parameter search engine
(13)
, 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
directoriented free applications (such as Scientist from
Micromath) to solve the problem of statistical determination
370 2.5
In [(TTmo)/(TTm)] = 0.001009.0 + 0.0034
360 , R2 0.9987
A.U/V.Cpv= 0.001009  2
350 
1.5 F7
340
E
" 330
E . 1
l _ =
320
, * *  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 
LJ
2I
. 330  . .
S 4 ^
320
E
) : c 
I�
: In[(TR  TF)/(TR TF)]= 0.1259.08
310 R2 = 0.9973
0 1/ = 0.9973s, 2
300
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 temperaturetime 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.
CONCLUSIONS
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.
NOMENCLATURE
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]
REFERENCES
1. Feisel, L.D., and G.D. Peterson, "The Challenge of the Laboratory
in Engineering Education," J. Eng. Ed., 91(4), 367 (2002)
TABLE 1
Comparison of Parameter Values
PARAMETER VALUE OBTAINED VALUE OBTAINED VALUE OBTAINED VALUE REPORTED IN
MODEL I MODEL II INDEPENDENT THE LITERATURE
EXPERIMENT
A 2.1010 2.1010  2.1010  6.85.10"
E (Cal/mol) 8120 8238  82389200
U (KCal/s m2 K) 0.0031 0.0035 0.0036 0.0020.02
T (s 1)  0.1259 0.1259
Vol. 42, No. 1, Winter 2008
292
288
284
a.
E 280
I
 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 ContinuousFlow 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
(1516), 2479 (2006)
10. Bhanot, S., S. Kumar, and M.K. Vasantha, "Online 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, "EndPoint 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 FedBatch 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 PeroxideSulfiteThiosulfate
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 2729 (2005)
19. Solodov A., and V. Ochkov, "Stiff Differential Equations" in Differen
tial Models. An Introduction with Mathcad, Springer Berlin Heidelberg
(2005)
20.
21. Perry R.H., and C.H. Chilton (Eds.), Chemical Engineers'Handbook,
5th Ed., McGrawHill 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
classroom
  ^ K.___________________________
ACTIVE PROBLEM SOLVING
AND APPLIED RESEARCH METHODS
In a Graduate Course on Numerical Methods
ERIC L. MAASE AND KAREN A. HIGH
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 problembased 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
lemsolving skills. The researchers developed workshops to
explicitly develop skills deemed appropriate for improving
students' problemsolving skill and confidence.
Prince,31] in his article on active learning, provides additional
motivation for incorporating problembased learning in the
classroom. One benefit that he shows is that PBL typically
involves significant amounts of selfdirected 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 preEn
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.
MODELING COURSE
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 (IVPODEs)
(k) Boundary Value Problems (BVPODEs)
(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 graduatelevel 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 problemsolving principles using Visual Basic for
Applications (VBA) linked through Microsoft Excel and
MATLAB have developed during previous semesters as
24
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 methodsboth 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
phenomena.
AN ACTIVE LEARNING ENVIRONMENT
In modification to the course the instructors introduced
a semesterlong 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 openended. Defining the problem and the over
all goals of the group projects is an aspect of the students'
learning.
Other educational objectives inherent in the problem state
ment include:
* Consider different types/scales of phenomena (micro
scopic, multiple gradient, maximum gradient, macro
scopic).
* 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
coordinates.
* 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
ingand in some cases, quite eagerto 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.
PROJECT GROUPS
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.
INTRODUCTION TO LABORATORY METHODS
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).
TABLE 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 studyproviding 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 openended 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
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 justintime
teaching[91 or lecturebydemand 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
26
aims of encapsulated drug delivery. This case study is dis
cussed after the necessary computer instruction (VBA/Excel
and MATLAB) is complete.
Case Study #2Modeling and Numerical
Simulation
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 realworld 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
TABLE 2
Laboratory Safety/Practices for 307 EN CHE 5743 Class Lab
1. Safety glasses (eye protection) are encouraged
2. Move deliberatelyno 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)
Signature
TABLE 3
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 12 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 neededyour 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 coauthor. 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
compartments.
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
Peripheral
c, A4 Figure 2.
= kpCp+kcC Case Study #2:
d1 Drug delivery
model.
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 STUDY USE
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 modelsin tractable mathematical
terms able to capture essential behaviors of the actual system
of interestwas 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.
PROJECT REPORTING
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.
FINAL REPORTS AND PRESENTATIONS
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 WORK AND COURSE
ASSESSMENT
Student Project Examples
One group undertook an investigation of the design of a
cost effective controlledrelease drug delivery system for
Digoxin.1101 In their efforts they proposed a design  a four
layer delayedrelease pillintended to sustain optimal blood
plasma concentration over a 24hour period; the proposed
TABLE 4
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
experiments?
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 inclass 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 nextlevel 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
effectiveness."
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
goals?
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.
fourlayer 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 endofcourse 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
sections.
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
Point.
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
procedures.
TABLE 5
Student Designs: FourLayer 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 problembased learning environment.
Q11. I can search the literature for relevant background
material for projects.
Q12. I understand the components of a research report /
proposal.
Q13. I know what goes into a developing a good presenta
tion.
The students were given a 5point 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 problembased 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
coauthors 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
fortably"
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 freeform
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
lembased 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
CONCLUDING COMMENTS
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, impracticalas 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.
REFERENCES
1. Smith, K., S. Sheppard, D. Johnson, and R. Johnson, "Pedagogies
of Engagement: ClassroomBased Practices," J. Eng. Ed., 94(1), 87
(2005)
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
(1970)
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, JustInTime
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)
TABLE 7
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
APPENDIX
The following are freeform 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
dissolution."
"I never knew how to model physical phenomenon before I
enrolled in this course. It helped me improve my computer
skills."
"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 almostreal
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 reallife situations."
"This seemed a bit vague .1i ....,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
skills."
"Not much, I was good at them."
"Gave me a better understanding of Runge Kutta and other
methods."
"The project gave opportunity to solve numerically the
ODEs."
"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
reallife situations."
"Really helped us to model realtime problems.
"I felt like the project was very numericalmethod 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
skills."
Q5. Describe your ability to develop a costeffective 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 costeffective encapsulated drug consider
ing mathematical models, economic evaluation models."
"I couldn't actually develop an exact costeffective 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 costeffective encapsulated drug."
"The costeffective 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
course."
"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...
STUDENT RATINGS OF TEACHING:
MYTHS, FACTS, AND GOOD PRACTICES
RICHARD M. FIELDER
North Carolina State University
REBECCA BRENT
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
quality.
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 multisection 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 falseit 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/felderpublic>.
Rebecca Brent is an education consultant
specializing in faculty development for ef
f fective university teaching, classroom and
computerbased simulations in teacher
education, and K12 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
teaching.
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, upperlevel courses
tend to get higher ratings than lowerlevel 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 2address 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
semester.
* Use a few global or summary items with Likertscale
(15) 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.
REFERENCES
1. Felder, R.M., "WhatDo They Know, Anyway?"( ,.... 3 /., 2( ,),
134135 (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, 315 (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., 326328 (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), 6978 (2007)
4. Greenwald, A.G., and G.M. Gillmore, "Grading Leniency is a Remov
able Contaminant of Student Ratings," American Psychologist, 52(11),
12091217 (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, 107118 (1980)
6. Felder, R.M., and R. Brent, "How to Evaluate Teaching," Chem. Eng.
Ed., 38(3), 200202 (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
curriculum
0
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*
Kansas State University * Manhattan, KS 66506
I n their contribution, MartinezUrreaga 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 reactionspecifically, 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 firstorder, 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 processnetwork synthesis
Urreaga and his collaborators[11 to enhance its usefulness, based on process graphs (Pgraphs). 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 birthdeath 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 birthdeath 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 Songtien Chou obtained his B.S. in chemical engineering from National
tuations. Second, the master equation of the birthdeath 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 timedriven approach[14 University in Taiwan; he was the department chair between 2000 and
5] in addition to the eventdriven 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,
YungKang City, Tainan Hsien, 71003 Taiwan � Copyright ChE Division of ASEE 2008
Vol. 42, No. 1, Winter 2008 3.
MartinezUrreaga and his collaborators11l on the basis of the
Gillespie algorithm.l0, 17 18] The eventdriven 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 eventdriven approach, the timedriven 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 eventdriven
approach is exceedingly complex to determine. The means
and variances of the random variables characterizing the
birthdeath process have been computed in implementing the
eventdriven and timedriven 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 MartinezUrreaga and his
collaborators.[1]
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.
MODEL FORMULATION
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)
k2
Clearly, the population of A decreases and that of B increases
when one molecule of A transforms into one molecule of
B with forward reactionrate 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
reactionrate constant k2. The system constitutes, therefore, a
special class of Markov processes, the birthdeath 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 birthdeath 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)
(2)
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}.
MASTER EQUATION
The derivation of the master equation of a stochastic process
has been extensively elaborated in our previous contribution
to this journal.12s5 For the birthdeath process of interest, the
probability balance around any arbitrary state nA of the inde
pendent random variable, NA(t), gives rise to1" 12]
d
P(nA;t)= W+ (nA 1;t)p(n 1;t)
dt
+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 firstorder
rate law, the intensity of birth, W+(nA; t), and the intensity of
death, W (nA; t), have been defined as"'
dn
W+ (nA;t) = = k2nB = k2 (nTO nA) (7)
dt
dn,
W (nA;t) dnA kn, ,A
dt
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 nonlinear 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]
MEANS AND VARIANCES
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 1exp[(k +k2)t]}
2 (k +k2)
[nA k + k2) + k2 exp[(k+ k2) t + klk (12)
"TO
For the case, nA0 = n, this equation becomes[301
S(kt)= k 1 {1exp(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),
1/2
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(kik) + k exp[(k +k2)t +kk2 = (t) (15)
From this expression, oB(t) is given by
1/2
TB (t)= nTO {1exp[(k k2)t}/2
(k, + k2)
2 22
[T ( k2) +k2 exp[(k +k2)t+ klk2
Vol. 42, No. 1, Winter 2008
TA(t)
Randomly or
stochastically
behaving
processes or
systems are
ubiquitous in
the chemical
or allied
industries.
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 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))
1/2
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
nTO
exp (k, + k2)t]
kk 1/2
kL k2
k2)t]+k2
k2')t /2
1/2
no (nAk k k) +k exp[ (k1 +k2)t +kk2
TO17)
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 exp2(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, MartinezUrreaga 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).
STOCHASTIC SIMULATION
The master equation of the birthdeath 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 3336] 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 eventdriven approach10 1523] and the other resorting to the
timedriven 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.
EventDriven Approach
The eventdriven 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 eventdriven 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 MartinezUrreaga and his collaborators, '1
has been revised and augmented to include the evaluation of
the means, variances, and standard deviations of the birthdeath
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, datarecording 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;
(20)
where W [nAZ(t);t] = k2[n TnA(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:
z
A.z A nAZ (0)
Z=1
Z
B,z ()  nBZ ()
z=1
c. Store the sum of squares of realizations at 6:
z
4A,Z Zn) n2
z=1
z
Zn,z 0) nB ,z (0)
Z=1
d. Store the square of sum of realizations at 6:
Z
4'BZK0) ZnB
AZ(BI,Z
z=l
[ A,Z (0)2
B,Z(0) 2
e. Compute the sample means at 6[13 19 20]:
1 z 1
mAZ AZ AAZZ6A
Z z1 Z
1 z 1
mBZ )z () 
MBZ()ZEB,Z() =BZ()
Z Zl Z
(21)
(22)
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(lu)
 W [nA, (t); t+W_[n,(t);t]
sz (e)
 � nAZ AZZ)  Z AZe)1 )
z 1) z znz O (Z 1) Z
(23)
sB,z (0)
Z1 Z
(Z n ))
 nz(O)
z Z1
(Z 1) z z (
(24)
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
eventdriven approach as presented above.
TimeDriven Approach
As briefly indicated at the outset of this contribution, the timedriven
40
approach[14 15] differs from the eventdriven 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 birthdeath process of inter
est in particular by the timedriven 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:
1
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:
z
Z
A.Z (t)  nAZ (t)
Z=1
Z
c. Store the sum of squares of realiza
tions at t:
z
QA.Z (t) Z T Z (t)
z=1
Z
^B.Z (t) n n.z (t)
z=1
Chemical Engineering Education
d. Store the square of sum of realizations at t:
z
Z
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)
71 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 timedriven approach as
presented above.
(30)
f. If 1 < Z < Z,, then compute sample the variances and standard
deviations at t:
z (tt)

Z1) nZzt(t)
sB.z (t)
n zt)
S(zi { B
(Z  1) z t1 )
1 z 2 I
Z 1aTz1 = (Z1)L
 i nB (t)
Z z1
1 )
_ 1) {^B.z(t)
(zi~l 
Z Az(t) (31)
Z1
I Z (t) (32)
Z
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  nA. (t)
Otherwise, no event occurs; thus,
nAz (t) nz (t At)
nB.Z (t) nT  nAz (t)
Vol. 42, No. 1, Winter 2008
RESULTS AND DISCUSSION
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 MartinezUrreaga and
his collaborators."1 Thus, k1 is 4  103 s 1
and k2 is 1  103 s1; 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 G1 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
tinezUrreaga 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 birthdeath
process is equivalent to the deterministic
model.
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 G1 as the stan
dard deviation envelopes, [(N (t)) � o(t)]
for N (t) and [(NB(t)) � o,(t)] for NB(t),
41
0 200 400 600 800 1000
t,s
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 G1:
As indicated by MartinezUrreaga 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 birthdeath 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
G2 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 G2 in
Appendix G; they are essentially identical to [(NA(t)) o,(t)]
42
0 200 400 600 800 1000
torO,s
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 eventdriv
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 G3 and G4 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 timedriven approach. For the first set of
nAO and n B, the results with Z = 2 are presented in Figure
G3, and those with Z = 50 are presented in Figure G4; they
are compared with the solutions of the deterministic and
stochastic models. Note that mA.z(t) and mBz (t) in Figure G4
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 G5 in Ap
pendix G.
Table 1 lists the computational times, required by
the Monte Carlo method via the eventdriven and
timedriven 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 eventdriven approach was shorter
than that required by the timedriven approach.
Figures 2 and 3 in conjunction with Figures G2
through G5 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 eventdriven and timedriven approaches, are
further illustrated with an example of photoelec
trochemical disinfection of bacteria whose rate
of decaying is assumed to be of the firstorder, 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 birthdeath
processes in which the intensity of birth is absent;
hence, it is termed the puredeath process.[10 11 25]
Vol. 42, No. 1, Winter 2008
In modeling the photoelectrochemical disinfection of bacteria as
a puredeath 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 puredeath 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)
dt
(35)
which is the master equation of the process11, 25]; see Appendix A. On
the basis of the firstorder rate law, the intensity of death, W (nc;t) ,
in the above expression is
dnc
W (nc;t) dnc  k3nc
dt
where k3 is the firstorder rate constant. For the puredeath 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)
and
7c (t)= nco exp(k3t)[1 exp(k3t)],
(38)
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
(39)
(40)
TABLE 1
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 EventDriven
and TimeDriven Approaches
Initial values Eventdriven approach* Timedriven 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
Notes:
* In the eventdriven approach, the sample means, variances, and standard
deviations have been recorded every two seconds, i.e., AO = 2 s.
� In the timedriven 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.
(36)
or, from Eq. (39),
(c (t)= lexp (k 3t)12 (Nc (t))1/2 (41)
From Eqs. (37) and (40), the coefficient of variation, CVc(t),
is computed as
1 1exp(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 puredeath 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 MartinezUrreaga 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 eventdriven and timedriven 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.
CONCLUDING REMARKS
The contribution of MartinezUrreaga 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 firstorder, i.e., linear, rate law as
a birthdeath process. The master equation arising from the
model has rendered it possible to analytically obtain the ex
44
pectedmeans and variances of the process. Second, the master
equation of the process has been stochastically simulated
through the Monte Carlo method via the timedriven approach
in addition to the eventdriven 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 firstorder 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.
ACKNOWLEDGMENT
This is contribution No. 071235, Department of Chemical
Engineering, Kansas Agricultural Experiment Station, Kansas
State University, Manhattan, KS 66506.
NOMENCLATURE
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)
10
t or 0, min
Figure 4. Temporal evolution of the mean
and standard deviation envelope near the tailend
of photoelectrochemical disinfection of a population
ofE. coli per unit volume: The average of 50
Monte Carlo simulations by resorting to the
eventdriven and timedriven 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
k,
k,
k,
NA(t)
NB(t)
Nc(t)
nA
nB
nc
nAO
nBO
nTO
nco
(NA(t))
(NB(t))
(Nc(t))
p(nA; t)
forward reactionrate constant, s 1
reverse reactionrate constant, s 1
firstorder 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)
E[NA(t)]
E[NB(t)]
E[Nc(t)]
probability of nA molecules being present at time
t
p(nc; t) probability of nc live bacteria being present at
time t
t
U
u
Var[NA(t)]
Var[NB(t)]
Var[Nc(t)]
W (nA; t)
time
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 birthdeath process in
state nA at time t
W (nA; t) intensity of death for the birthdeath process in
state n at time t
W (nc; t) intensity of death for the puredeath process in
state nc at time t
Greek Letters
0 datarecording time in the eventdriven 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 eventdriven approach
REFERENCES
1. MartinezUrreaga, J., J. Mira, and C. GonzalezFernandez, "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, "BirthDeath 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, "SurfaceRenewal 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 ContinuousFlow 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 Processesan Introduction for Physical
Scientists, Academic Press, San Diego, pp. 48, 60, 226, 328, 330, 375,
380(1992)
11. van Kampen, N.G., Stochastic Processes in Physics and ( I ...
NorthHolland, Amsterdam, pp. 5558, 9697, 134136, 139, 163
(1992)
12. Oppenheim, I., K. E. Shuler, and G. H. Weiss, Stochastic Processes in
Chemical Physics: The Master Equation, The MIT Press, Cambridge,
MA, pp. 5361 (1977)
13. Sobol', I.M., A Primerfor the Monte Carlo Method, CRC Press, Boca
Raton, FL, pp. i1, ix, 1516, 42 (1994)
14. Rod, V., and T. Misek, "Stochastic Modeling of Dispersion Formation
in Agitated LiquidLiquid Systems," Transactions of the Institution of
Chemical Engineers, 60, 48 (1982)
15. Rajamani, K., WT. Pate, and D.J. Kinneberg, "TimeDriven and
EventDriven Monte Carlo Simulations of LiquidLiquid Dispersions:
A Comparison," Ind. Eng. Chem. Fund., 25, 746 (1986).
16. Kendall, D.G., "AnArtificial Realization of a Simple 'BirthandDeath'
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
SteadyState 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 EnzymeSubstrate 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 NonLinear
Sieving Kinetics," Chem. Eng. Comm., 61, 59 (1987)
27. Fox, R.O., and L.T. Fan, "Application of the Master Equation to
45
Coalescence and Dispersion Phenomena," Chem Eng Sci., 43, 655
(1988)
28. Chou, S., L.T. Fan, A. Argoti, R. VidalMichel, 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, 8991, 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. 97102 (1989)
32. Pang, WK., PK. Leung, W.K. Huang, and W Liu, "On Interval
Estimation of the Coefficient of Variation for the ThreeParameter
Weibull, Lognormal, and Gamma Distribution: A SimulationBased
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
(2005)
41. Fan, L.T., H.T. Chen, and D. Aldis, "An Adaptive Random Search
Procedure for LargeScale Industrial and Process Systems Synthesis,"
Proceedings of the Symposium on Computers in Design and Erection
of Chemical Plants, Aug. 31Sep. 4, Karlovy Vary, Czechoslovakia,
pp. 279291 (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. 2223 (1999)
45. McQuarrie, D.A., Mathematical Methods for Scientists and Engineers,
University Science Books, Sausalito, CA, pp. 4, 1053 (2003)
APPENDICES
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
THE PENGROBINSON
EQUATION OF STATE
Using Mathematica
HOUSAM BINOUS
National Institute of Applied Sciences and Technology
Consider a single component characterized by its critical
pressure, temperature, acentric factor, Po T , and o. The
PengRobinson equation of state[2 4] is given by:
RTa
P = (1)
(Vb) [V(V+b)+b(Vb)]
where
R T
b = 0.07780 R
Pc
(RTo )2
a= 0.45724 L 1+m  )
Pc
T
Tc
m= 0.37464 + 1.54226  0.26992 2. (5)
The PengRobinson 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 multicomponent mixture:
3 + ( B)Z2 +(A 3B2 2B)Z + (AB +B2
* Tunis, Tunisia
where
c c c c
A=EEy.yA, or EExxA,
I1 ji1 I1 j1i
A = AAA 0 (1k )
C C
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 fulltime 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
(6)
*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 doublespaced pages and should be accompanied
by the originals of any figures or photographs. Please submit them to Professor James O. Wilkes
(email: wilkes@umich.edu), Chemical Engineering Department, University of Michigan, Ann
Arbor, MI 481092136.
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)
v
where the fugacity coefficients are:
(Zv 1)  In(ZvB) A
B 241B
v, = exp 2 yA B 1 B1 (13)
1   In +(+52)
A B Z  B
The vaporphase and liquidphase fugacity coefficients,
which measure deviation from ideality, are defined from the
liquidphase and gasphase 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 gasphase
compressibility factor, Zv, with its liquidphase 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
enthalpy:
HD = RT(Z  1) 1+
24B RT
P
L +(I + 2)B T d A(RT)2 (RT)2' (
LoZ+ 1 )B dT P PA (15)
ADIABATIC FLASH CALCULATIONS FOR
MULTICOMPONENT MIXTURES
Problem statement
A quaternary mixture, at 485 psia and 100 �F, is composed
of 0.41% hydrogen, 5.71% methane, 70.97% benzene, and
48
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 vaporphase fraction, the temperature
and other relevant variables.
Solution
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 builtin function FindRoot and are
given by:
Eq[i_]:= y[i]
EQ[i_]:= x[i] == z[i]
4
== K[i] x[i]
/(1+ vap (K[i]  1));
4
sEq[1],Eq[2],Eq[3], y[i]
sol = FindRootEQI,EQ[2,EQ[3
EQ[I1],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,
1I1
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},
{x[2],0.03},{x[3],0.7},{x[4],0.25},
,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 vaporphase fraction vap and temperature are 0.0367
and 560.297 �R, respectively.
SOLUBILITY OF METHANOL IN NATURAL GAS
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.
Solution
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 builtin function FindRoot to solve
13 nonlinear algebraic equations simultaneously. The Math
ematica commands, which allow the determination of the 13
unknowns, are the following:
o
E
0
D 0.001
0.00001
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,
T=
419.67,437.67,455.67,473.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,
I=1
Lz[i] (K[i] 1)/(1 + vap (K[i]  1))== 0,
I1
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[i1]},
{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
methanol.
HIGHPRESSURE CHEMICAL EQUILIBRIUM
COMPUTATION
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 gasphase
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(2X) 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.
Solution
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 builtin 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}
{y[2],0.2},{y[1],0.1},{y[3],0.4},{X,0.6},
{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
pressures.
CONCLUSION
We have shown through simple examples how one can
take advantage of the numerical and graphical capabilities
of Mathematica to perform flash calculations for multicom
ponent mixtures, solubility determination, and highpressure
chemical equilibrium computations. Similar computations can
be performed very easily using Matlab. 7
NOMENCLATURE
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.
1.0
0
S0.8
0.4
900� R
0.2 \ 1080� R
1260� R
1440� R
0.0
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
HD:
k:
C:
z:
R:
a:
v:
vap
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
REFERENCES
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 FluidPhase Equilibria, 3rd Ed., PrenticeHall,
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, EquilibriumStage Separation Opera
tions in Chemical Engineering, Wiley, NY (1981)
7.
do?objectType= author&objectId=1093893> 7
Vol. 42, No. 1, Winter 2008
This onepage 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.
CELEBRATING ChE IN SONG
ALAN M. LANE
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 endofsummerlab 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 andif you dareto use in your classroom. I
believe there are many wouldbe 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 email (alane@eng.ua.edu), and you just
might find it posted on this Web site and played on chemical
engineers' iPods around the globe!
REFERENCES
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 doublespaced pages and should be accompanied
by the originals of any figures or photographs. Please submit them to Professor James O. Wilkes
(email: wilkes@umich.edu), Chemical Engineering Department, University of Michigan, Ann
Arbor, MI 481092136.
CAN I TRUST THIS SOFTWARE PACKAGE?
An Exercise in Validation
of Computational Results
MORDECHAI SHACHAM
BenGurion University of the Negev * Beer I,. .., 84105, Israel
NEIMA BRAUNER
TelAviv University * TelAviv 69978, Israel
W. ROBERT ASHURST
Auburn University * Auburn, Alabama 36849
MICHAEL B. CUTLIP
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 BenGurion 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 TelAviv 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 TelAviv University. Her research interests include
numerical methods or programming. In most cases, it is suf hydrodynamics and transport phenomena in twophase 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 nanoscale, and the influence of surface chemical treatments
The failure of a software package to reach a correct solu on micro and nanoscale 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 cochaired the ASEE Summer School for chemical engineering faculty in
Copyright ChEDivi2002. He is also coauthor 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
tion.
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.
ANALYZING THE BEHAVIOR OF THE
RUNGEKUTTA ALGORITHM FOR AN
"INTERACTING TANKS" SIMULATION
PROBLEM
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)
dh2
q2 + q4  q5 = Te2h2 (2)
dt
dh3
q5 + q3  q6 = h 3 (3)
dt
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: #).
TABLE 1
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
(m^2)
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, .....
corn>)
2MATLAB is product of The Math Works, Inc. ()
Chemical Engineering Education
The system is started at steadystate 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'
TABLE 2
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.39E05 8.66E06 9.39E05 1.17E05
dlh/dt 6.22E06 6.22E06 3.06E05 1.16E05
dh/dt 1.11E04 5.39E07 2.71E04 7.08E06
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 RungeKutta
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.
RECOMMENDED STEPS AND
ADDITIONAL INFORMATION
FOR SOLUTION
1. Most software packages include some
variation of the RungeKutta 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
dh,
dt
f= dh (7)
dt
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. (13),
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.
SOLVING THE PROBLEM WITH A DIFFERENT
SOFTWARE PACKAGE
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 ifelse
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.21
1.2 
S1.19
E
Cl,
' 1.18
J
_1
2 1.17
1.16
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 = 103).
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.
CHECKING THE ROLE OF THE TRUNCATION
ERROR TOLERANCE
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
4000
TABLE 3
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)]);
13
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,
nonoscillatory 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
I TABLE
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 squareroot" error message of the
RKF45 algorithms are caused by the lack of consideration
of the slope value errors by the stepsize control logic of the
two algorithms.
CHECKING THE ROLE OF
THE PROBLEM STIFFNESS
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
checked.
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
Tolerance
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.0015
0.001 
0.0005
0
0.0005
0.001 
0.0015 
0 500 1000
1500
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
e08 time period. The symbolic
5 differentiation is carried out
'e05
in lines 5 through 9, and
6.0
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.37297e4; X2 
0.04478; and 3 = 0.106922. Introduc
ing these values into the stiffnessratio
equation yields:
SR = = 286.67
X,
As the stiffness ratio is less than
1000, the system is not considered to
be stiff according to the stiffnessratio
criterion.
Chemical Engineering Education
Figure 4. Error in the liquid level h, (ode45 Algorithm, RelTol = 10o).
CONCLUSIONS
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 nonstiff
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 twopart 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 reallife example. This
exercise demonstrates that a complete 2.00E04
analysis and validation of the results
are important in problem solving and 1.50E04
for maintaining confidence in a particu
lar software package. While students 1.00E04
may use software packages for routine *
problem solving with a limited knowl 5.00E05
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.00E05
when difficulties are encountered during
numerical problem solving. 1.00E04
ACKNOWLEDGMENTS 1.50E04
The authors wish to thank Auburn 2.00E04
University chemical engineering stu 0
dent Robert Edwards for first bringing
this problem to our attention.
Vol. 42, No. 1, Winter 2008
REFERENCES
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),
169(1998)
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
TABLE 5
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);
f(i),hl);
f(i),h2);
f(i),h3);
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;
bs(df);
_numeric)
g
ure . rror n e qu eve s ope: (o e gor m, e o = .
59
500 1000 1500 2000 2500 3000
Time (s)
ht l;i ;il Irl l1 dh /dt A d+ ( Ar ih RlTr l 1 n3)
