Modeling and Simulation in Python

Chapter 15

License: Creative Commons Attribution 4.0 International



In [1]:

    
# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *

The coffee cooling problem

I'll use a State object to store the initial temperature.



In [2]:

    
init = State(T=90)

And a System object to contain the system parameters.



In [3]:

    
coffee = System(init=init,
                volume=300,
                r=0.01,
                T_env=22,
                t_0=0,
                t_end=30,
                dt=1)









    Out[3]:







  
    
      
      values
    
  
  
    
      init
      T    90
dtype: int64
    
    
      volume
      300
    
    
      r
      0.01
    
    
      T_env
      22
    
    
      t_0
      0
    
    
      t_end
      30
    
    
      dt
      1

The update function implements Newton's law of cooling.



In [4]:

    
def update_func(state, t, system):
    """Update the thermal transfer model.
    
    state: State (temp)
    t: time
    system: System object
    
    returns: State (temp)
    """
    r, T_env, dt = system.r, system.T_env, system.dt
    
    T = state.T
    T += -r * (T - T_env) * dt
    
    return State(T=T)

Here's how it works.



In [5]:

    
update_func(init, 0, coffee)

Here's a version of run_simulation that uses linrange to make an array of time steps.



In [6]:

    
def run_simulation(system, update_func):
    """Runs a simulation of the system.
    
    Add a TimeFrame to the System: results
    
    system: System object
    update_func: function that updates state
    """
    init = system.init
    t_0, t_end, dt = system.t_0, system.t_end, system.dt
    
    frame = TimeFrame(columns=init.index)
    frame.row[t_0] = init
    ts = linrange(t_0, t_end, dt)
    
    for t in ts:
        frame.row[t+dt] = update_func(frame.row[t], t, system)
    
    return frame

And here's how it works.



In [7]:

    
results = run_simulation(coffee, update_func)

Here's what the results look like.



In [8]:

    
plot(results.T, label='coffee')
decorate(xlabel='Time (minutes)',
         ylabel='Temperature (C)')

And here's the final temperature:



In [9]:

    
coffee.T_final = get_last_value(results.T)
T_final = get_last_value(results.T)









    Out[9]:





72.2996253904031

Encapsulation

Before we go on, let's define a function to initialize System objects with relevant parameters:



In [10]:

    
def make_system(T_init, r, volume, t_end):
    """Makes a System object with the given parameters.

    T_init: initial temperature in degC
    r: heat transfer rate, in 1/min
    volume: volume of liquid in mL
    t_end: end time of simulation
    
    returns: System object
    """
    init = State(T=T_init)
                   
    return System(init=init,
                  r=r, 
                  volume=volume,
                  temp=T_init,
                  t_0=0, 
                  t_end=t_end, 
                  dt=1,
                  T_env=22)

Here's how we use it:



In [11]:

    
coffee = make_system(T_init=90, r=0.01, volume=300, t_end=30)
results = run_simulation(coffee, update_func)
T_final = get_last_value(results.T)









    Out[11]:





72.2996253904031

Exercises

Exercise: Simulate the temperature of 50 mL of milk with a starting temperature of 5 degC, in a vessel with the same insulation, for 15 minutes, and plot the results.

By trial and error, find a value for r that makes the final temperature close to 20 C.



In [12]:

    
# Solution

milk = make_system(T_init=5, t_end=15, r=0.133, volume=50)
results = run_simulation(milk, update_func)
T_final = get_last_value(results.T)









    Out[12]:





20.00135627897414



In [13]:

    
plot(results.T, label='milk')
decorate(xlabel='Time (minutes)',
         ylabel='Temperature (C)')



In [ ]:

	T
0	90
1	89.32
2	88.6468
3	87.9803
4	87.3205
5	86.6673
6	86.0207
7	85.3804
8	84.7466
9	84.1192
10	83.498
11	82.883
12	82.2742
13	81.6714
14	81.0747
15	80.484
16	79.8991
17	79.3201
18	78.7469
19	78.1795
20	77.6177
21	77.0615
22	76.5109
23	75.9658
24	75.4261
25	74.8919
26	74.3629
27	73.8393
28	73.3209
29	72.8077
30	72.2996