Engineering With Python Rotating Header Image

Unsteady Heat Flow in a Slab – Simple Parabolic Partial Differential Equation (PDE)

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 Carnahan, Luther and Wilkes, “Applied Numerical Methods”, Wiley NY 1969 pg 434.  Very similar to an earlier problem that we did from Cutlip and Shacham.

Problem Description:

Assume an infinite parallel-sided slab with a constant thermal diffusivity, α, that is initially at a uniform temperature T0.  At time t, two sides are changed to a constant temperature T1.  So, you want to find how the temperature inside the slab varies with time.

So, redefine the variables as follows:

T = (T – To)/(T1 – T0), dimensionless temperature

τ = (αt)/L2, dimensionless time

X = x/L, dimensionless length

Parabolic Partial Differential Equation:

∂T/∂τ = ∂²T/∂X²

Boundary Conditions:

τ = 0: T = 0 for 0 ≤ X ≤ 1

τ > 0: T = 1 at X = 0 and X = 1

Solution Approach:

In this case, we’ll use an ‘explicit approach’ and replace the differentials with selected finite difference forms.  For a detailed explanation, go to Finite Difference Method at wikipedia.  Now, would you really use this approach to solve a simple PDE?  Probably not.  There are other approachs you would probably use first.  I like the method of lines approach where the spatial variable is replaced with an appropriate finite difference form and then use an established initial value solver.  That way you can develop a ‘template’ for solving ODE’s/PDE’s that is similar and quick.  Also, there are python packages such as ‘escript finley’ from the Earth Systems Science Computational Centre (ESSCC) at the University of Queensland which are very interesting and deserve a post.  There are non-python based, like Elmer Multiphysics, FreeFem and PDESolutions to name a few.  If you are aware of others, let me know and we’ll discuss them.  However, I want to stick to ‘open source’, or available in a ‘free’ version.  Now, let’s look at the programming solution.

Programming Solution:

#!usr/bin/env python
#========================================================================
# Example 7.1 Carnahan, Luther and Wilkes, "Applied Numerical Methods"
# pg.434
# Unsteady-State Heat Conduction in an Infinite, Parellel-Sided Slab
# Solved using an explicit method
#========================================================================
#Import Packages for Solution
#========================================================================
from scipy import *
from pylab import *
from matplotlib import *
#========================================================================
#Input constants
#========================================================================
dtau = 0.001 # Set dimensionless time increments
dx = 0.05    # Set dimensionless length increments
Tmax = 0.95  # Set maximum dimensionless temperature
M = 21       # Counter for length discretization
#========================================================================
#Calculate parameters
#========================================================================
dx = 1.0/(M-1)
dx_x = 1.0/(M-1)
ratio = dtau/(dx**2)
const = 1.0 - 2.0*ratio
#========================================================================
#Set counters to zero
#========================================================================
i = 0
tau = 0.0
#========================================================================
# Set up arrays for solution and print
#========================================================================
Tnew = zeros(M, dtype = float)
T = zeros(M, dtype = float)
T[0] = 1.0
T[-1] = 1.0
print "T initial = ", T
#========================================================================
# I just pick 400 on trial and error for the total array
#========================================================================
T_sol = zeros((400,M), dtype = float)
T_sol[:,0] = 1.0
T_sol[:,-1] = 1.0
#========================================================================
# While loop to iterate until mid-point temperature reaches Tmax
#========================================================================
while T[10] < Tmax:
    i = i + 1
    tau = tau + dtau
#========================================================================
# Calculate new tempertures
#========================================================================
    for j in range(1,M-1):
        Tnew[j] = ratio*(T[j-1] + T[j+1]) + const*T[j]
#========================================================================
# Substitute new Temperatures in array for T
#========================================================================
    for k in range(1,M-1):
        T[k] = Tnew[k]
        T_sol[i,k] = T[k]

print "Tau and T_final =", tau, T_sol[i,:]
#========================================================================
# Set up array for spatial values of x to plot
#========================================================================
x = [i*dx_x for i in range(M-1)]
x.append(1.0)
#========================================================================
# Plot the solutions
#=======================================================================
plot(x,T_sol[50,:])
plot(x,T_sol[100,:])
plot(x,T_sol[150,:])
plot(x,T_sol[250,:])
plot(x,T_sol[i,:])
#legend(['Tau = 0.5','Tau = 0.1','Tau = 0.15','Tau = 0.25',
#'Tau = final time'])
title('Normalized Slab Temperatures')
grid()

Now, since my last post oh so many months ago (about three), I’ve been experimenting with different ways to present the python code without using a screen shot.  I’m still learning, and  I’ve started using a ‘syntax highlighter’ above.  It presents  the line numbers, and it seems to be keeping the indentation.  You should be able to ‘highlight’ the code segment and paste into your editor (say SciTE), then save and run.  I’m also putting more comments into the code so that they would go along with the paste.  I’m always interested in feedback, let me know what you think.

Problem Solution:

Here is the plot of the solution, not much to say about it.  It compares favorably to the one on Carnahan.

T initial =  [ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. 0.  0.  1.]
Tau and T_final = 0.327,  [ 1.          0.99218968  0.98457167  0.97733356  0.97065357  0.96469619
0.9596081   0.9555146   0.95251648  0.95068755  0.95007287  0.95068755
0.95251648  0.9555146   0.9596081   0.96469619  0.97065357  0.97733356
0.98457167  0.99218968  1.        ]

So, the output above is the initial conditions and the final answer (including the time step) from IPython.   I’ll post the conventional solution for the elliptic PDE next.  I haven’t had a lot of success with the nitric acid flowsheet.  That one is probably better solved in a spreadsheet.  If anyone has a good Python solution for that one, send it to me and I’ll let you guest post.  Thanks for reading my blog.

Solving a Binary Batch Distillation – Solution

Solve the Problem:

Let’s run the program and see what the plot looks like:

henleySeader362FinalPlot

Well, not much more to say.  From the plot, x_w(t = 5 hours) is about 0.35.  Henley and Seader calculated xw (t = 5 hours) = 0.351.  I’d say pretty close.  It’s easy to calculate W(t).   It is 100 kg – 10 kg/hr*5 hr = 50 kg.  So, this one is pretty easy.  I’m not sure what to do next, I’ll think about it.  Any suggestions?

Solving a Binary Batch Distillation – 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, list out your equations and problem data)
  2. Remember [Accumulation = In – Out + Source/Sink]
  • Think about what you need to do and the answer you want
  • You need to solve for an initial value ordinary differential equation, so you’ll need an ODE solver
  • 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

henleySeader362Startup

Input the Constants:

henleySeader362InputConst

Define the Function for the ODE Solver:

henleySeader362FuncDef

Set up to Plot the Solution:

henleySeader362Plot

This was pretty bare bones.  Pretty straightforward.  Let’s press on to the solution.

Solving a Binary Batch Distillation – 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

www.wikipedia.org - Batch Distillation

www.wikipedia.org - Batch Distillation

Engineering, Vol. 11, No. 6, pp. 723-736, but they took it from Henley and Rosen, “Materal and Energy Balance Computations”, Wiley, NY 1969 pg. 362.  Very similar to an earlier problem that we did from Cutlip and Shacham.

Problem Description:

We need to find the composition of a ‘batch distillation‘ as a function of time.  Essentially, we need to solve the following ODE:

dxw/dt = -[D/W(t)]*(y – xw)

t = time, hours

xw = liquid composition, dimensionless

W(t) = -D*t + W(0), total amount remaining in ‘pot’, kg

D = overhead rate, kg/hr

y = α*xw/[1 + xw *(α – 1)], vapor composition, dimensionless

α = relative volatility, dimensionless

Initial Conditions:

D = 10.o kg/hr

W(0) = 100 kg

What is the composition at t = 5 hr?   Not so bad, we just need to solve a single initial value ODE.  Piece ‘o cake.

Solution of Equilibrium Equations – Solution

Solve the Problem:

Let’s see what the output of the program looks like:

henleyRosen2FinalAnswer

As a check, I added up the mole fractions to ensure they add up to 1.0.  If they didn’t, something is wrong.  Let’s compare to the what Rosen and Adams calculated.

CH4 R&A = 0.226  EWP = 0.226

H2O R&A = 0.110  EWP = 0.110

CO R&A = 0.0997  EWP = 0.0997

H2 R&A = 0.512  EWP = 0.512

CO2 R&A = 0.053  EWP = 0.053

Tot R&A = 15.843  EWP = 15.843

Hmm, I would say the answers are the same.  They converged in 20 iterations with their spreadsheet.  We took a different approach and used a different solver and we got the same answer, so odds are is that it is correct.  Let’s do a batch distillation for our next problem and see if it is like the one we have already solved.  See you then.

Solution of Equilibrium Equations – 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, list out your equations and problem data)
  2. Remember [Accumulation = In – Out + Source/Sink]
  • Think about what you need to do and the answer you want
  • You need to solve for the roots of several equations, so you need a solver 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

henleyRosen1ImportPackages

Input Problem Parameters and Initial Conditions:

henleyRosen2InitialCond

Note, that x_int, the inital ‘advancement of reaction’, is input as a list with the A1,1 = Rxn 1 and A1,2 = Rxn 2.

Define the Equilibrium Reactions in a Function:

henleyRosen2SolveFunction

Note that the ‘fsolve’ function is included here.  It uses the function, ‘equil’ and the initial conditions for x.

Defining the Final Mole Amounts and Print (Get ready… it looks complicated):

henleyRosen2ConcentrationPrint

Now, this looks pretty complex… and it is.  It is really bookkeeping.  Keeping track of all of the compounds and how they are changing with respect to how the reaction proceeds.  Also, since the lines were long I used the ‘\’ continuation operator to break the command into two lines.  You have to be really careful and for several equations, it can be very long.  Is that what graduate students are for?  Sorry, I suffered as a grad student, it must be in the guild rules.  Stay around and we’ll post the answer in the next post.  You’ll say, we went to all that trouble for that!!

Solution of Equilibrium Equations – 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 Rosen, “Materal and Energy Balance Computations”, Wiley, NY 1969 pg. 376.  I will also be using information from C. G. Hill, “An Introduction to Chemical Engineering Kinetics and Reactor Design”, John Wiley and Sons, NY 1977, Chapter 2

Problem Description:

Rosen and Adams were using the spreadsheet to solve for the ‘roots’ of a series of equilibrium material balances.  Typically, you have equilibrium coefficients (Kp) that you use to set up the equations for solution.  If you are lucky, you also know how the K-values change with temperature and pressure.   They didn’t really describe much beyond what components were involved and that they were going to use a damped Newton-Raphson method.  They used the spreadsheet to set up the data, take the first derivative and then iterate to solve for the roots.

Essentially, you’ll be solving for the roots of the equilibrium material balance equations:

CH4 + H2O = CO + 3H2 (Reaction 1)

CO + H2O = CO2 + H2 (Reaction 2)

Here, we are writing the equilibrium mole number in terms of an “equilibrium degree of reaction advancement_i)”.  Essentially, it is based on the stoichiometry and is defined as (Hill, Chapter 2):

ni = ni0 + [νi*xi]rxn1 + [νi*xi]rxn2 + …..

ni = equilibrium number of moles component i

ni0 = initial number of moles component i

νi = stoichiometric coefficients of component i in rxn 1

xi = reaction ‘advancement’ coefficient of rxn 1

and the others defined accordingly.

Equilibrium Number of Moles:

n_CH4 = n_CH4_int – x_1
n_H2O = n_H2O_int – x_1 – x_2
n_CO = n_CO_int + x_1 – x_2
n_H2 = n_H2_int +3*x_1 + x_2
n_CO2 = n_CO2_int + x_2

Total Number of Moles (add up all of the components):

n_tot = (n_CH4_int + n_H2O_int + n_CO_int + n_H2_int + n_CO2_int) + 2*x_1

Remember Partial Pressures (Use CH4 as an example):

yCH2 = (n_CH4/n_tot)*P

P = pressure, atm

other components defined the same way

Equilibrium Reactions Written to Solve for Roots:

f1 = (y_CO*P)*(y_H2*P)**3 – k_1*(y_CH4*P)*(y_H2O*P)
f2 = (y_CO2*P)*(y_H2*P) – k_2*(y_CO*P)*(y_H2O*P)

Problem Data:

Kp1 = 0.54 atm

Kp2 = 2.49 atm

P = 1.0 atm

Initial Number of Moles:

n_CH4_int = 6.0 moles

n_H2O_int = 5.0 moles

n_CO_int = 0.0 moles

n_CO2_int = 0.0 moles

n_H2_int = 0.0 moles

Okay, whew… I think that is all we need to start setting up to solve the problem.  See you at the next post.

Calculate a Four Component ‘Flash’ – Solution

Solve the Problem:

Let’s see what the output of the program looks like:

henleySeader279Output

We can see that the mole fractions are propane = 0.071, butane = 0.183, pentane = 0.310 and hexane = 0.435 and they (whew!) add up to 1.0 which is a good sign.  We calculate α = 0.1219 and Rosen and Adams calculate α = 0.1219.  I would say that is pretty close, even for government work!  I don’t know which I’ll post next, come back and find out.  Thanks for visiting my blog.  See you and take care.

Calculate a Four-Component ‘Flash’ – 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 solve for the root of an equation, so you need a solver 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

henleySeader279ImportPackages

Input the Data for Solving the Problem (Scipy arrays again, oh boy!):

henleySeader279InputData

Note the x = [].  That initializes a list named ‘x’ so that we can fill it with data later.  Also, I have a list called ‘name’.  More on that in a moment.

Set up Function for Solution:

henleySeader279SetUpFunction

Here we have set up the function for solution.  It is a ‘brute force’ approach.  Each portion is spelled out with the component name.  You could probably streamline it with a ‘for’ loop or something similar.  I’ll leave that as an ‘exercise for the reader’.  For me, doing it this way takes a little longer up front, but when I come back days, weeks, months later, I can easily figure out what the heck I was doing.  Also, note the data in Scipy arrays and the ‘sum’ is the return value.

Solve the Equation and Print the Solution:

henleySeader279SolveandPrintSol

Note the  x.append().  That fills our ‘x’ liste with the individual values that we can then print and sum.  The sum is to ensure the values add up to 1.0 as a check on the solution.  If they don’t add up to 1.0 exactly, you don’t have the right answer and can’t pass Go or collect $200.  Now, we have our program, we can proceed to the solution.  Stay tuned…

Calculate a Four-Component ‘Flash’ – 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. 279.

Problem Description:

Rosen and Adams were using the spreadsheet to solve for the ‘roots’ of the flash equation.  You calculate a flash to calculate a temperature or pressure at the top or bottom of a column, or perhaps the feed condition or a knock-out drum or if you like to solve equations.  They didn’t really describe much beyond what components and that they were going to use Newton’s method.  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, take the first derivative and then iterate to solve for the root.

Essentially, you’ll be solving for the roots of the equilibrium flash equation:

f(α) = ∑(zi*(1 – Ki)/(1 + α*(Ki – 1)))

at a specified temperature and pressure.

Ki = yi/xi equilibrium constant of component i (Dimensionless)  Also, assume K is not a function of composition for this problem, typically they are, but not this time…

yi = mole fraction of component i in vapor phase (Dimensionless)

xi = mole fraction of component i in liquid phase (Dimensionless)

zi = feed mole fraction of component i  (Dimensionless)

α = fraction of feed that is vaporized (it is β in the wikipedia article on ‘flash evaporation’)

Problem Data:

K_propane = 4.2

K_butane = 1.75

K_pentane = 0.74

K_hexane = 0.34

Temperature and pressure are constant

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