0

I'm using a system of ode's to model coffee bean roasting for a class assignment. The equations are below.

Systems of Equations

The parameters (other than X_b and T_b) are all constants.

When I try to use odeint to solve this system, it gives a constant T_b and X_b profile (which conceptually doesn't make sense).

Below is the code I'm using

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

# Write function for bean temperature T_b differential equation
def deriv(X,t):
    T_b, X_b = X 
    dX_b = (-4.32*10**9*X_b**2)/(l_b**2)*np.exp(-9889/T_b)
    dT_b = ((h_gb*A_gb*(T_gi - T_b))+(m_b*A_arh*np.exp(-H_a/R_g/T_b))+ 
    (m_b*lam*dX_b))/(m_b*(1.099+0.0070*T_b+5*X_b)*1000)
    return [dT_b, dX_b]

# Establish initial conditions
t = 0 #seconds
T_b = 298 # degrees K
X_b = 0.1 # mass fraction of moisture

# Set time step
dt = 1 # second

# Establish location to store data
history = [[t,T_b, X_b]]

# Use odeint to solve DE
while t < 600:
  T_b, X_b = odeint(deriv, [T_b, X_b], [t+dt])[-1]
  t += dt
  history.append([t,T_b, X_b])

# Plot Results
def plot_history(history, labels):
  """Plots a simulation history."""
  history = np.array(history)
  t = history[:,0]
  n = len(labels) - 1
  plt.figure(figsize=(8,1.95*n))
  for k in range(0,n):
    plt.subplot(n, 1, k+1)
    plt.plot(t, history[:,k+1])
    plt.title(labels[k+1])
    plt.xlabel(labels[0])
    plt.grid()
  plt.tight_layout()


plot_history(history, ['t (s)','Bean Temperature $T_b$ (K)', 'Bean Moisture Content $X_b$'])
plt.show()

Do you have any ideas why the integration step isn't working?

Thank You!!

Richard
  • 56,349
  • 34
  • 180
  • 251
L Buck
  • 3
  • 1
  • 2
  • Your equations' implementation looks wrong. Shouldn't you have, for example, `dX_t` instead of `dX_b` and `np.exp(-9889/(T_b+273.15))` instead of `np.exp(-9889/T_b)`? I'm expecting a very clean correspondence between your written equations and their implementation in code, but I'm not seeing that. – Richard Apr 25 '18 at 15:20
  • There is some difference in notations between the equations (my apologies) because they come from different literature sources X_b and X are the same, and the T_b+273.15 is already accounted for previously. – L Buck Apr 25 '18 at 15:25
  • Having the correct equations would be a better way to get help. I've also reformatted your code and added the necessary library includes so that it could run. Unfortunately, you leave many constants undefined, so I get error messages (`NameError: name 'l_b' is not defined`). Please fix this. – Richard Apr 25 '18 at 15:36

1 Answers1

0

You're repeatedly solving the system of equations for only a single timepoint.

From the odeint documentation, the odeint command takes an argument t which is:

A sequence of time points for which to solve for y. The initial value point should be the first element of this sequence.

Since you pass [t+dt] to odeint, there is only a single timepoint so you get back only a single value which is simply your initial condition.

The correct way to use odeint is similar to the following:

output = odeint(deriv, [T_b, X_b], np.linspace(0,600,600))

Here output, again according to the documentation is:

Array containing the value of y for each desired time in t, with the initial value y0 in the first row.

Richard
  • 56,349
  • 34
  • 180
  • 251
  • Ah, thank you! I can tell this is calculating the function correctly now. – L Buck Apr 25 '18 at 15:55
  • Great! If you feel this answered your question, feel free to either upvote it using the arrow to the left of the answer or accept it by clicking the outline of a checkmark. – Richard Apr 25 '18 at 16:01
  • Is there a way to plot the every value that is output? I.e. the value of X_b and T_b at every time step? – L Buck Apr 25 '18 at 16:05
  • @LBuck: Yes, you can do that pretty much as you were. The `np.linspace` yields a list of times against which you can plot the values in `output`. – Richard Apr 25 '18 at 16:18