Engineering With Python Rotating Header Image

‘Graphical’ Evaluation to Calculate the Number of Transfer Units – Solution

Solve the Problem:

Calculate the number of transfer units (Nog)

Run the program and let’s see what the output shows.  It will show up in your IPython window.

Number of Transfer Units based on Gas Phase Resistance

Number of Transfer Units based on Gas Phase Resistance

The answer we get is NOG = 3.52 transfer units.  Henley and Seader got 3.44 transfer units (obtained graphically) and Rosen and Adams calculated 3.73 units (with a spreadsheet).  So, essentially, we all got the same answer within the error of the data and the different methods.  I show the array values so you can compare with Rosen and Adams if you want (and have the article, get it, it’s a good one.)  We like short and sweet, let’s proceed to another problem.  Ciao!

‘Graphical’ Evaluation to Calculate the Number of Transfer Units – Programming Approach

Plan your solution:

  • Draw a picture, in this case, list all of your data
  • Remember the fundamentals and apply
  1. Draw your material or energy balance envelope (If necessary, not in this case)
  2. Remember [Accumulation = In – Out + Source/Sink]
  • Think about what you need to do and the answer you want
  • You need to perform a numerical integration of data, so you need an integration package
  • You’ll want to print the solution to get the information

How to start your program:

  • With Ipython open, open your editor
  • Label your program (you’ll never remember it, go ahead and label it!)
  • Import the packages you’ll need for solving

hsImportPackages

I looked at the Scipy.integrate documentation and found an integration package that could work with unevenly spaced data.  I chose the Romberg integration package, but you could chose one of the others and try them.

hsInputData

I entered the data for x and x* into Scipy arrays.  As we go on, you’ll see why this is such a nice feature.  You can work with the arrays and perform all kinds of stuff easily.  So, why don’t I quit talking and show you.

Manipulate the Data to get ready for the Integration:

hs653ManipulateData

You see, we can use the ‘arrays’ just like variables and work with them easily.  Very cool.

Set up the Integration and Print the Solution:

hs653RomIntegration

So, that’s is pretty easy, don’t you think.  We’ll discuss the solution in the next post.

‘Graphical’ Evaluation to Calculate the Number of Transfer Units – Problem Description

Reference:

This is from E.M. Rosen and R. N. Adams, “A Review of Spreadsheet Usage in Chemical Engineering Calculations”, Computers and Chemical Engineering, Vol. 11, No. 6, pp. 723-736, but they took it from Henley and Seader, “Equilibrium Stage Separation Operations in Chemical Engineering”, Wiley, NY 1981 pg. 653.  I don’t have a copy of the book, so, I may not include everything.  If you do and want to add to this, please do.  Also, why do we chem e’s always have to have such long titles?  And why do we always have to include ‘chemical engineering’ in the title?  Like we’ll forget or something.

Problem Description:

Rosen and Adams were using the spreadsheet to set up an interpolation and ‘table look up’ function with the data.  Usually, you use this approach when you have a ‘packed’ column that doesn’t have ‘discrete’ trays/stages.  They didn’t really describe the column arrangement or the components involved.  Okay, okay, I’ll go the library and get a copy of the Henley and Rosen book so that I know more about the problem.  They used the spreadsheet to set up the data, interpolate and estimate so they had values at even locations.  With that data table, they set up a four-interval Newton-Coates integration to calculate the number of transfer units.  Here is what the data looks like set up using arrays in Scipy.

X and X* Equlibrium Data as Scipy Arrays

X and X* Equlibrium Data as Scipy Arrays

You usually get this type of data from the equilibrium operating diagram for the system you are working with.  X* is the equilibrium concentration.  Consult Henley and Seader or look at a Perry’s Chemical Engineer’s handbook for more detail.  Treybal, “Mass Transfer Operations” is another good text you may want to consult.

Essentially, you’ll be doing and integration as follows:

N(OG) = ∫(1 – x)lm/((1 – x)*(x – x*)) dx

with (1 – x)lm = ((1 – x) – (1 – x*))/(ln((1 – x)/(1 – x*)))

N(OG) = Number of transfer units based on the gas phase resistance to mass transfer

This is what we have, let’s post on our programming approach.

A New Set of Problems to Work…From the Early Days of Personal Computers

What Next?

from www.wikipedia.org - 'Personal Computers'

from www.wikipedia.org - 'Personal Computers'

We just completed working the ten problems from the 1997 ASEE workshop.  We worked a variety of problems that had solutions with other software packages such as Microsoft Excel, MATLAB, Mathematica, Maple and Polymath.  Please compare and let me know what you think.  Now, let’s look at the following journal article for a ‘journey into the past’.  In 1987, as personal computers were getting started, E. M. Rosen and R. N. Adams published an article in ‘Computers and Chemical Engineering’, Vol. 11, No. 6, pp. 723-736 describing solution approaches for the following problems:

  • Numerical integration to calculate the number of transfer units (stages)
  • Four component flash
  • Solve a ‘two reaction’ gas-phase equilibrium
  • Determine a molar volume from the Van Der Waal’s equation
  • Binary differential distillation
  • Steady state heat conduction in a plate
  • Un-steady state heat conduction in a slab
  • Material balance around a 5-unit nitric acid process

It was ‘cutting edge’ for its time and still relevant now as a way to understand basic numerical approaches for solving different types of chemcial engineering problems.  They used the ‘Visicalc’ spreadsheet program to set up and solve these problems.  Gosh, I’m old enough to remember Visicalc…well, that’s another set of posts.  You know, I sold my IBM PC 5150 a couple of years ago on eBay.  Really, I did!  You had to boot from the floppy and use DOS.   Thank God those days are over, for this we will use our IPython/Scipy/Matplotlib set (and maybe we’ll use Sage….) on our nice modern laptop, while listening to Grooveshark….

So, let’s begin…

We will follow the problem pattern that they used.  I urge you to get a copy of the article and read it.  These guys are pretty smart and we can still benefit from their knowledge.  So, let’s start with the next post about integrating to calculate transfer units.  We’ll see how things have changed (or haven’t changed in 20+ years….).

Fitting Polynomials and Correlation Equations to Vapor Pressure Data – Solution

Solve the problem:

1. Fit the Data to a Polynomial

Below is the program output in IPython with values for the various data fits.

prob13IPythonOutputIt is a little confusing, but now you know why in the earlier post that the polynomial coefficients are ‘backwards’.   You see the polynomial defined as p = 0.0007449*x3 + 0.03945*x2 + 1.198*x + 24.46 instead of the common way of p = 24.46 + a1*x + a2*x2 +….

Now let me define these items (we want the lowest sse and the highest R-squared):

sse_3 = 204.733367 ( ‘sum of squared errors’ for the 3rd order polynomial)

R2_vp3 = 0.99961988 (‘R squared’ for 3rd order polynomial)

sse_log = 0.321997 (‘sum of squared errors’ for the Clapeyron equation fit)

R2_log_vp = 0.9915015 (‘R squared’ for the Clapeyron equation fit)

sse_ant = 0.0118258 (‘sum of squared errors’ for the Antione equation fit)

R2_ant_vp = 0.99968788 (‘R squared’ for the Antione equation fit)

Here is the plot for the polynomial fit

prob133rdOrderPolyFit

Not a bad fit from the visual examination and the R squared is good, but the sse is pretty large, at lower temperature and pressure the fit is not quite as good as at higher temperature and pressure.  Also, notice I forgot to put labels on this plot, but from the values Temperature is the y-axis and Pressure is the x-axis.  Let’s look at the Clapeyron equation fit.

2. Fit the Data to the Clapeyron Equation:

prob13clapeyronPlotAgain, I forgot to label the axis, but you can figure out which is which.  Not good practice, but remember, “Do as I say, not necessarily as I do…”.  Here the fit is not as good visually, but the sse is lower which means it fits better (errors cancel out) over the span of the plot.  Also, I used the 1st order polyfit as the easiest way to fit the transformed data.   Now, let’s look at the Antione equation.

3. Fit the Data to the Antoine Equation:

prob13antioneEqnPlotWow, the transformed data visually fits well and the sse is the lowest and R-squared is the best.  Also, I used the scipy.optimize.leastsq package to regress the transformed data.  See how well the data ‘fits’ over the entire span of the data.   Now, we know why we use Antione’s Equation to fit vapor pressure data.

What’s Next?

Well, I’ve completed the 10 problems of the ASEE 1997 problem set.  Was it useful for you?  I learned a fair amount and I hope you did too.  I intend to work a few more problems and then go into using GUI’s and grabbing data from the web.  For my next set of problems, I will work a few from the following journal article by Rosen, E.M. and Adams, R.N., “A Review of Spreadsheet Usage in Chemical Engineering Calculations”, Computers in Chemical Engineering, Volume 11, No. 6, pp. 723-736.  Come back and join me.  Even better, comment freely and share any problems you have worked.  Thanks.

Fitting Polynomials and Correlation Equations to Vapor Pressure Data – Programming Approach

Plan your solution:

  • Draw a picture, in this case, list all of your data and equations
  • Remember the fundamentals and apply
  1. Draw your material or energy balance envelope (If necessary, not in this case)
  2. Remember [Accumulation = In – Out + Source/Sink]
  • Think about what you need to do and the answer you want
  • You need to perform linear regression, so you want a ‘least squares’ package
  • You’ll want to plot the solution to get the information
  • Need a regression solver (scipy.optimize), statistical support(scipy.stats) and plotting package (Matplotlib,pylab)

How to start your program:

  • With Ipython open, open your editor
  • Label your program (you’ll never remember it, go ahead and label it!)
  • Import the packages you’ll need for solving

prob13importPackages

Next, input your vapor pressure data into an array so that you can begin to manipulate the data.  Arrays are ‘objects’ (here is a good link) and are very flexible and useful for managing numbers.

Input your Vapor Pressure Data:

prob13inputArrayData

Set up and try the Polynomial Fit (try 3rd Order):

prob13polynomialFit

Now, you’ll notice that the polynomial coefficients are ‘backwards’ on line 38 above due to the way the coefficient output is put into an array.  I haven’t checked but there may be packages in the stats package that calculate the residual difference squared and R squared fit, but I went ahead and put them in for good measure.

Transform the Data and fit to the Clapeyron Equation:

prob13clapeyronEqn

Note that the lists have defined as ‘arrays’ and set type to ‘float’.  I don’t completely get the Python Types yet, but often I find setting them to float can help.

Transform the Data and fit to the Antione Equation:

prob13antioneEqn

Now, it’s ready to run.  In the next  post, we’ll examine the output and use that to answer the questions.

Fitting Polynomials and Correlation Equations to Vapor Pressure Data – Problem Description

Reference: This is Problem 3 in the ASEE 1997 problem set.   We will be fitting vapor pressure data for benzene to  a polynomial and then to a couple of correlating equations.

Concepts:

  • Put the data into arrays so that we can work with it
  • Use Ipython with Scite editor to edit and save files
  • Import packages into workspace
  • Define constants
  • Use ‘lists’ to work with data sets
  • Use the scipy.stats ‘mean’ package
  • Use the scipy.optimize ‘leastsq’  package
  • Use the matplotlib package to plot the results

Problem Description:

We will fit vapor pressure data for benzene to a polynomial and then perform a linear regression of the data using the Claperyon equation and the Antoine equation.  Both are correlating equations that have been used to fit vapor pressure data for various compounds.

Data and Equations:

Vapor Pressure Data for Benzene

Temperature, T deg C Pressure, P mm Hg

-36.7                                                    1

-19.6                                                     5

-11.5                                                   10

-2.6                                                   20

7.6                                                   40

15.4                                                   60

26.1                                                 100

42.2                                                200

60.6                                               400

80.1                                                760

General Polynomial Expression

P(x) = a0 + a1*x + a2*x2 + a3*x3 + . . . +an*xa

a0 through an are the coefficients calculated by a least squares fit of a selected degree n of the polynomial

Clapeyron Equation

log(P) = A + B/T

T = Absolute Temperature in deg K and A and B are parameters that are determined from a least squares fit of the transformed data

Antoine Equation

log(P) = A – B/(C + T)

T = Absolute Temperature in deg K and A, B and C are parameters that are determined from a non-linear regression of the transformed data.

Questions:

  1. Correlate the data with a couple of different degree polynomials assuming that the absolute temperature is the independent variable and P is the dependent variable.  Determine which polynomial fits the data best.
  2. Correlate the data using the Clapeyron equation.
  3. Correlate the data using the Antione equation.

Now, I saved this one for last since using a spreadsheet (Microsoft Excel, Openoffice Calc, Gnumeric) is easier or use an open source statistical package such as R.  The program for this one is a little messy.  I suspect it can be streamlined, so any suggestions are welcome.

Dynamics of a Heated Tank with Proportional/Integral (PI) Control – Solution

Solve the problem:

Okay, let’s open our IPython window, pull up our script and set Kc = 0.0.  Also, be sure that both the proportional and integral parts of your controller is on.  It doesn’t make a lot of difference since Kc is zero, but it keeps everything straight.  Also, ensure that Ti is ‘on’ for the step and Tr is ‘off’ so it is a constant.

Problem 10 Part 1 Plot

Problem 10 Part 1 Plot

Essentially, you have turned off the controller and therefore the heater.  The tank equilibrates after a little more than a half hour at the feed temperture. The plot of the heater output is what?  It is a straight line at zero.  Now you know what the ‘open loop’ system performancd looks like.  Let’s move on. Continue reading →

Dynamics of a Heated Tank with Proportional/Integral (PI) Control – Programming Approach

Plan your solution:

  • Draw a picture, in this case, list all of your data and equations
  • Remember the fundamentals (material and energy balances, equations of state) and apply
  1. Draw your material or energy balance envelope (Proportional/Integral Contoller)
  2. Remember [Accumulation = In – Out + Source/Sink]
  • Think about what you need to do and the answer you want
  • You need to solve a set of ordinary differentials equations with supporting algebraic equations
  • You’ll want to plot your results

How to start your program (we do this every time…):

  • With Ipython open, open your editor
  • Label your program (you’ll never remember it, go ahead and label it!)
  • Import the packages you’ll need for solving

Import the packages you need for solution:

prob10importPackages

I like to import everything I think I’ll need.  It is annoying enough sometimes working out programming problems.  So, I try to catch it all up front.

Define your constants, your initial conditions and the time interval for your ODE solver:

prob10constantsIniConditions

This time, I put in all of the parameters, initial conditions and the time start and end (I wasn’t lazy like I have been in past posts).  Also, notice that three empty lists t_plot, q_plot and qlim_plot in the constants list.  Later you will see that we define those items inside the function for our ODE solver.  However, remember that variables defined inside a function are only available inside that function.  So, this is one way to bring the information out so that it can be plotted later.  There are other ways, such as declaring a global variable, but this works and it never hurts to learn to work with lists.  Also, note that I put in an ‘on’ and an ‘off’.  This allows to you to turn the proportional and integral parts of your controller ‘on’ and ‘off’ as you need. Continue reading →

Dynamics of a Heated Tank with Proportional/Integral (PI) Control – Problem Description

Reference:

This is problem 10 from an article  “A Collection of Representative Problems in Chemical Engineering for Solution by Numerical Methods”.  These problems came from Session 12 of the ASEE Chemical Engineering Summer School held in Snowbird, Utah in 1997.  In this problem, you will be modeling the temperature of a tank heated with a heater controlled with a PI controller.

Concepts:

  • Solve a system of ordinary differential equations (with a step function) and appropriate initial conditions
  • Use Ipython with Scite editor to edit and save files
  • Import packages into workspace
  • Define constants
  • Define functions for the ODE solver
  • Use the odeint solver with initial conditions and time span array
  • Print the solution

Problem Description:

A well-mixed tank is being used to continuously heat a solution for further processing.  The process set-up is shown below:

Well-mixed Heated Tank with PI Controller

Well-mixed Heated Tank with PI Controller

The temperature in the tank is not high enough or changes enough to worry about properties changing with temperature and the feed rate is constant.  You want to change the temperature of the outlet stream by changing the setpoint on the controller.  The parameters are:

W = Feed rate to tank, kg/min

Cp = Heat capacity of the solution, kj/kg – deg C

rho = Density of the solution, kg/m3

Ti = Inlet temperature, deg C

Tr = Setpoint temperature, deg C

Tm = Measured temperature, deg C

V = Tank volume, m3

q = Heat input, kj/min

Mr. Phelps, your mission, should you decide to accept it, is to figure out how to control the temperature of the outlet stream of the tank by changing the setpoint on the PI controller.  This blog will self-destruct in 30 sec….

Problem Equations and Set-up:

Energy Balance for the tank:

dT/dt = [W*Cp*(Ti – T)/(rho*V*Cp)] + q

Thermocouple Response First Order plus Dead Time:

To(t) = T*(t – τd)

t = Time, min

τd = Dead time constant, min

Pade Approximation for the Dead Time Effect (connect the temperature input to the thermocouple To):

dTo/dt = [T – To – (τd/2)*(dT/dt)]*(2/τd)

Model the Thermocouple Shield and Electronics:

dTm/dt = (To – Tm)/τm

τm = Thermocouple time constant, min

Control the Energy Input to the Tank with a PI Controller:

q = qs + Kc*(Tr – Tm) + (Kc/τi)*∫(Tr – Tm)dt

Kc = Proportional Gain of the Controller

τi = Integral Time Constant, min

Energy Required to Maintain Steady-State:

qs = W*Cp*(Tr – Tis)

Define the Error Term as a First Order ODE:

d(errsum)/dt = Tr – Tm     errsum = 0 at t = 0 steady-state (before you mess it up by changing something)

Substitute into PI Controller Equation:

q = qs + Kc*(Tr – Tm) + (Kc/τi)*(errsum)

Parameters for Solution:

rho*V*Cp = 4000 kj/deg C             W*Cp = 500 kj/min-deg C

Tis = 60 deg C                                        Tr = 80 deg C

tau_d = 1 min                                        tau_m = 5 min

Kc = 50 kj/min-deg C                         tau_i = 2 min

Questions:

  1. Caculate the open loop performance (set Kc = 0) of this system when the system is initially operating at design steady state at a temperature of 80°C, and inlet temperature Ti is suddenly changed to 40°C at time t = 10 min (the dreaded ‘step’ change). Plot the temperatures T, T0, and Tm to steady state, and verify that Padé approximation  for 1 min of dead time is working properly.
  2. Calculate the closed loop performance of the system for the conditions of part 1 and the baseline parameters from the Table. Plot temperatures T, T0, and Tm to steady state.
  3. Repeat part 2 with Kc = 500 kJ/min×°C.
  4. Repeat part 3 for proportional only control action by setting the term Kc/tI = 0.
  5. Put limits on q so that the maximum is 2.6 times the baseline steady state value and the minimum is zero. Demonstrate the system response from baseline steady state for a proportional only controller when the set point is changed from 80°C to 90°C at t = 10 min. Kc = 5000 kJ/min×°C. Plot q and qlim versus time to steady state to demonstrate the limits. Also plot the temperatures T, T0, and Tm to steady state to indicate controller performance