**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 T_{0}. At time t, two sides are changed to a constant temperature T_{1}. So, you want to find how the temperature inside the slab varies with time.

So, redefine the variables as follows:

T = (T – T_{o})/(T_{1} – T_{0}), dimensionless temperature

τ = (αt)/L^{2}, 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.