Frolian's blog standing on the shoulders of PyGiants

• About
• Archive

Solving the Heat Equation like in Harvard's Cook my Meat App

Tue, 14 Jan 2020

Food Numerical Physics Physics

Two years ago, I've had the pleasure to take a Harvard class entitled Science & Cooking: From Haute Cuisine to Soft Matter Science (part 1) on edx, a MOOC platform.

This has lead me to do some fun kitchen science, in particular trying to cook an egg using math formulas (spoiler: I was not successful).

One of the tools introduced by the class is the Cook my Meat App. The app allows you to specify a cooking protocol for meat modelled as a 1D medium. Ever since seing this app, I thought it would be nice to replicate it with my own tools and learn from it. This post is the attempt to do just this.

I'll try to go from the theory (the heat equation in 1D) to the implementation using the Crank-Nicolson time stepping method, in Python. In the end, we'll try to compare our results with these from the app (and they'll hopefully be similar).

I'll also apply these simulations to some other cooking problems.

The 1D heat equation and the Crank-Nicolson method

The heat equation using the Crank-Nicolson discretization

The problem tackled by the Cook my Meat app is that of cooking a steak. The steak has thickness $L$ and can be heated on one side or both sides by a heat source, while being in contact with the outside air on the other side. It also starts at a given initial temperature.

Our problem then consists of solving for the temperature distribution as a function of time. The temperature evolution is governed by the 1D heat equation with a constant diffusion constant, written as

$$ \frac{\partial T}{\partial t} = D \frac{\partial^2 T}{\partial x^2} $$

The physical interpretation of this equation is that heat flows through space along the $x$ coordinate at a rate that is proportional to the diffusion constant $D$ and to the local second derivative with respect to space of temperature.

The Crank-Nicolson method consists, when discretizing the continuous temperature field to a $N + 1$ nodes grid (with notation $T_i^n$, where $i$ represents the space dimension and $n$ the time dimension), to use the following equations

$$ \frac{T_i^{n+1} - T_i^{n}}{\Delta t} = \dfrac{1}{2} \left [ D \frac{T_{i+1}^{n} - 2 T_{i}^{n} + T_{i-1}^{n}}{\Delta x^2} + D \frac{T_{i+1}^{n+1} - 2 T_{i}^{n+1} + T_{i-1}^{n+1}}{\Delta x^2}\right ] \, \forall i \in [1, N-1] $$

Rearranging terms to have time $n+1$ on the left and time $n$ on the right and using the notation $\sigma = \frac{D \Delta t}{2 \Delta x ^2}$, this equation becomes

$$ - \sigma T_{i-1}^{n+1} + (1 + 2 \sigma) T_{i}^{n+1} - \sigma T_{i+1}^{n+1} = - \sigma T_{i-1}^{n} + (1 - 2 \sigma) T_{i}^{n} - \sigma T_{i+1}^{n} $$

This is the general term in the medium that will be used to update temperatures across time and space.

Boundary conditions

However, our problem still lacks boundary conditions. The two edges can have either a fixed temperature or a zero heat flux condition, depending on whether there is heat transfer (through a pan for instance, fixed temperature) or not (contact with air, zero heat flux).

To write the boundary conditions properly, we will use a little trick, which I read about here, based on the above equation. Let's tackle the left boundary first. Writing the equation for $i=0$, we obtain:

$$ - \sigma T_{-1}^{n+1} + (1 + 2 \sigma) T_{0}^{n+1} - \sigma T_{1}^{n+1} = - \sigma T_{-1}^{n} + (1 - 2 \sigma) T_{0}^{n} - \sigma T_{1}^{n} $$

The peculiar thing here is the appearance of a virtual node of index -1.

Since we wish to have a zero heat flux condition at the left boundary, we would write $\frac{\partial T}{\partial x} |_{x=0} = 0$ in the continuous system.

In the discrete system, we can use the virtual node and a finite difference to write

$$ \frac{T_{0}^{n+1} - T_{-1}^{n+1}}{\Delta x} = \frac{T_{0}^{n} - T_{-1}^{n}}{\Delta x} = 0 $$

This allows us to obtain

$$ (1 + \sigma) T_{0}^{n+1} -\sigma T_{1}^{n+1} = (1 - \sigma) T_{0}^{n} - \sigma T_{1}^{n} $$

In the case of a fixed temperature on the left, denoted $T_l$, we can just write $T_{0}^{n+1} = T_l$ and put it on the first line of the system.

Doing the same thing at the right boundary yields the following formula in the case of zero heat flux right boundary

$$ -\sigma T_{N}^{n+1} (1 + \sigma) T_{N+1}^{n+1} = -\sigma T_{N}^{n} (1 - \sigma) T_{N+1}^{n} $$

and $T_{N}^{n+1} = T_r$ on the last line of the linear system.

Assembling the whole system

Finally, the whole system becomes

$$ A \underline{T}^{n+1} = B \underline{T}^{n} + C $$

To see what the $A$, $B$ and $C$ matrices look like in an exemple, we can write them out in the case of a fixed temperature on the left and a zero heat flux boundary condition on the right

$$ \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0\\\\ -\sigma & 1+2\sigma & -\sigma & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & -\sigma & 1+2\sigma & -\sigma & \cdots & 0 & 0 & 0 & 0 & 0 \\\\ 0 & 0 & \ddots & \ddots & \ddots & \ddots & 0 & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\sigma & 1+2\sigma & -\sigma \\\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\sigma & 1+\sigma \end{bmatrix} \begin{bmatrix} T_0^{n+1} \\\\ T_1^{n+1} \\\\ T_2^{n+1} \\\\ \vdots \\\\ T_{N}^{n+1} \\\\ T_{N+1}^{n+1} \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0\\\\ \sigma & 1-2\sigma & \sigma & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & \sigma & 1-2\sigma & \sigma & \cdots & 0 & 0 & 0 & 0 & 0 \\\\ 0 & 0 & \ddots & \ddots & \ddots & \ddots & 0 & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \sigma & 1-2\sigma & \sigma \\\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \sigma & 1-\sigma \end{bmatrix} \begin{bmatrix} T_0^{n} \\\\ T_1^{n} \\\\ T_2^{n} \\\\ \vdots \\\\ T_{N}^{n} \\\\ T_{N+1}^{n} \end{bmatrix} + \begin{bmatrix} T_l \\\\ 0 \\\\ 0 \\\\ \vdots \\\\ 0 \\\\ 0 \end{bmatrix} $$

Assuming the matrix $B$ is invertible, the solution at timestep $n+1$ is obtained by

$$ \underline{T}^{n+1} = A^{-1} (B \underline{T}^n + C) $$

Since material parameters are constant, we can compute $B^{-1} A$ only once and then apply it repeatedly, which should be computationally efficient.

Implementation using Python

The units we will use for the implementation of our simulations are as follows: millimeters, seconds, ° Kelvin.

In this unit, the heat diffusion constant of water (to which steak and other food is assimilated) is $D = 0.14$ (in mm^2/sec).

import holoviews as hv
hv.extension('bokeh', logo=False)
hv.opts.defaults(hv.opts.QuadMesh(colorbar=True, width=500, height=300, tools=['hover'], cmap='viridis'))

import numpy as np
from enum import Enum

KELVIN  = 273.15   

class BC(Enum):
    """A class for holding boundary conditions."""
    HEATING_LEFT = 0
    HEATING_RIGHT = 1
    HEATING_BOTH = 2
    NO_HEATING = 3

def assemble_A_matrix(N, D, dx, dt, boundary_condition):
    sigma = D * dt / (2 * dx**2)
    A = np.zeros((N+1, N+1))
    for i in range(N+1):
        if i == 0:
            # imposed temperature
            if boundary_condition == BC.HEATING_LEFT or boundary_condition == BC.HEATING_BOTH:
                A[i, i] =  1
            # zero heat flux    
            elif boundary_condition == BC.HEATING_RIGHT or boundary_condition == BC.NO_HEATING:
                A[i, 0] = 1 + sigma
                A[i, 1] = - sigma
        elif i==N:
            # imposed temperature
            if boundary_condition == BC.HEATING_RIGHT or boundary_condition == BC.HEATING_BOTH:
                A[i, i] = 1
            # zero heat flux
            elif boundary_condition == BC.HEATING_LEFT or boundary_condition == BC.NO_HEATING:
                A[i, i-1] = - sigma
                A[i,   i] = 1 + sigma
        else:
            A[i, i-1] =       - sigma
            A[i,   i] = 1 + 2 * sigma
            A[i, i+1] =       - sigma
    return A

def assemble_B_matrix(N, D, dx, dt, boundary_condition):
    sigma = D * dt / (2 * dx**2)
    B = np.zeros((N+1, N+1))
    for i in range(N+1):
        if i == 0:
            # imposed temperature
            if boundary_condition == BC.HEATING_LEFT or boundary_condition == BC.HEATING_BOTH:
                B[i, 0] =   0
            # zero heat flux    
            elif boundary_condition == BC.HEATING_RIGHT or boundary_condition == BC.NO_HEATING:
                B[i, 0] = 1 - sigma
                B[i, 1] =     sigma
        elif i==N:
            # imposed temperature
            if boundary_condition == BC.HEATING_RIGHT or boundary_condition == BC.HEATING_BOTH:
                B[i, i] = 0
            # zero heat flux
            elif boundary_condition == BC.HEATING_LEFT or boundary_condition == BC.NO_HEATING:
                B[i, i-1] =     sigma
                B[i,   i] = 1 - sigma
        else:
            B[i, i-1] =         sigma
            B[i,   i] = 1 - 2 * sigma
            B[i, i+1] =         sigma
    return B

def assemble_C_matrix(N, D, dx, dt, boundary_condition, Tl, Tr):
    sigma = D * dt / (2 * dx**2)
    C = np.zeros((N+1,))
    if boundary_condition == BC.HEATING_LEFT or boundary_condition == BC.HEATING_BOTH:
        C[0] = Tl
    if boundary_condition == BC.HEATING_RIGHT or boundary_condition == BC.HEATING_BOTH:
        C[N] = Tr
    return C

def bc_type_checker(Tl, Tr, Tair):
    EPS = 0.001
    if abs(Tl - Tair) < EPS and abs(Tr - Tair) < EPS:
        return BC.NO_HEATING
    elif abs(Tl - Tair) < EPS:
        return BC.HEATING_RIGHT
    elif abs(Tr - Tair) < EPS:
        return BC.HEATING_LEFT
    else:
        return BC.HEATING_BOTH

class CrankNicolsonSolver:
    def __init__(self, x, T_initial, dt, D, verbose=False):
        self.t = 0
        self.dt = dt
        self.x = np.asarray(x).copy()
        self.T = np.asarray(T_initial).copy()
        self.N = self.x.size - 1
        self.dx = self.x[1] - self.x[0]
        self.bc_is_set = False
        self.D = D
        self.last_bc = None
        self.verbose = verbose
        
    def set_boundary_conditions(self, Tl, Tr, Tair):
        if self.last_bc is not None:
            if self.last_bc == (Tl, Tr, Tair):
                return
        bc = bc_type_checker(Tl, Tr, Tair)
        if self.verbose: print(f'[Set_BC] bc found: {bc}')
        A = assemble_A_matrix(self.N, self.D, self.dx, self.dt, bc)
        self.B = assemble_B_matrix(self.N, self.D, self.dx, self.dt, bc)
        self.C = assemble_C_matrix(self.N, self.D, self.dx, self.dt, bc, Tl, Tr)
        if self.verbose: print(f"[Set_BC] inverting matrix A with conditioning {np.linalg.cond(A)}")
        self.Ainv = np.linalg.inv(A)
        self.bc_is_set = True
        self.last_bc = (Tl, Tr, Tair)
        
    def step(self):
        assert self.bc_is_set
        self.T = np.dot(self.Ainv, np.dot(self.B, self.T) + self.C)
        self.t += self.dt
        
    def n_steps(self, n):
        for i in range(n):
            self.step()
            

To check that our solver is working, let's look at the solver output after a couple of steps of a sample problem.

x = np.linspace(0, 30, num=200) # 30 mm steak
T_initial = np.ones_like(x) * (23 + KELVIN)  # °K
T_hot = 45 + KELVIN
T_cold = 23 + KELVIN
T_air = 23 + KELVIN
dt = 0.01
D = 0.14
solver = CrankNicolsonSolver(x, T_initial, dt, D)
solver.set_boundary_conditions(T_hot, T_cold, T_air)

T0 = solver.T.copy() - KELVIN
solver.n_steps(5000)
T1 = solver.T.copy() - KELVIN
solver.n_steps(10000)
T2 = solver.T.copy() - KELVIN

hv.Curve((x, T0), label='T0').redim(y='temperature').opts(padding=0.05) * hv.Curve((x, T1), label='T1') * hv.Curve((x, T2), label='T2')

This looks like what we are expecting: heat slowly starts diffusing from the left boundary to the right.

With the solver in place, we can use it to simulate some cooking. First, we're going to make some functions that allow us to visualize the heat diffusion process.

Then, we'll adapt our solver to inputs from the Cook my Meat App.

Visualizing the heat flow on reference cases

To verify that our solver is doing the right thing, we're going to simulate some reference cases:

• fixed temperature left / zero flux on the right
• zero flux left / zero flux right with bumpy initial conditions

def make_snapshots(solver_params, boundary_conditions, tmax, n_snapshots, M=20):
    time_grid = np.linspace(0, tmax, num=n_snapshots)
    dt = time_grid[1] - time_grid[0]
    params = solver_params.copy()
    params['dt'] = dt / M
    solver = CrankNicolsonSolver(**params)
    solver.set_boundary_conditions(**boundary_conditions)
    snapshots = []
    for _ in range(time_grid.size):
        solver.n_steps(M)
        snapshots.append(solver.T.copy())
    return time_grid, snapshots

nx = 50
solver_params = dict(x=np.linspace(0, 30, num=nx), 
                     T_initial=np.ones_like(np.linspace(0, 30, num=nx)) * (23 + KELVIN), 
                     dt=.0001, 
                     D=0.14)

boundary_conditions = dict(Tl=150 + KELVIN, 
                           Tr=23 + KELVIN, 
                           Tair=23 + KELVIN)


tmax = 4 * 60
n_snapshots = 64
time_grid, snapshots = make_snapshots(solver_params, boundary_conditions, tmax, n_snapshots)

hv.QuadMesh((time_grid, solver_params['x'], np.array(snapshots).T - KELVIN), 
            kdims=['time', 'x'], 
            vdims='temperature')

The above behaviour is what we can expect in this situation: heat slowly diffuses from the hot boundary, getting slower as time advances.

solver_params = dict(x=np.linspace(0, 30, num=nx), 
                     T_initial=np.ones_like(np.linspace(0, 30, num=nx)) * (23 + KELVIN) + 10 * np.exp(-0.01 * (np.linspace(0, 30, num=nx) - 10)**2), 
                     dt=.0001, 
                     D=0.14)

boundary_conditions = dict(Tl=23 + KELVIN, 
                           Tr=23 + KELVIN, 
                           Tair=23 + KELVIN)

time_grid, snapshots = make_snapshots(solver_params, boundary_conditions, tmax, n_snapshots)

hv.QuadMesh((time_grid, solver_params['x'], np.array(snapshots).T - KELVIN), 
            kdims=['time', 'x'], 
            vdims='temperature')

Here again, this behaviour is what we expect: a bump in temperature slowly diffuses to the surrounding medium making it more homogeneous as time progresses.

Using cook my meat app inputs

Now that the basic machinery is in place, we can turn to replicating the data from the Cook My Meat app. It relies on a simple text format to describe cooking procedures.

The default recipe proposed by the app is as follows:

3cm Steak starts at 23°C
150°C and 23°C for 4:00
23°C and 150°C for 4:00
23°C and 23°C for 5:00

It consists of a header and an initial temperature followed by boundary temperatures and the duration of cooking with these in place.

We can build a simple parser using regular expressions that extracts the information from such an input.

import re

four_minutes_a_side = """3cm Steak starts at 23°C
150°C and 23°C for 4:00
23°C and 150°C for 4:00
23°C and 23°C for 5:00"""

def parse_input(recipe):
    """Parses the Cook my meat text format and returns a dictionary with data properly formatted."""
    header, *lines =  recipe.split('\n')
    thickness_cm, T_initial = list(map(float, re.findall("([\d.]+)cm Steak starts at ([-\d.]+)°[C]", header)[0]))
    T_initial += KELVIN
    steps = []
    for line in lines:
        Tleft, Tright, mins, secs = re.findall("([-\d.]+)°[C] and ([-\d.]+)°[C] for ([\d]+):(\d{2})", line)[0]
        steps.append((float(Tleft) + KELVIN, float(Tright) + KELVIN, 60 * int(mins) + int(secs)))
    parsed_params = dict(thickness_mm=thickness_cm * 10, initial_temperature_C=T_initial, steps=steps)
    return parsed_params

parse_input(four_minutes_a_side)

{'thickness_mm': 30.0,
 'initial_temperature_C': 296.15,
 'steps': [(423.15, 296.15, 240),
  (296.15, 423.15, 240),
  (296.15, 296.15, 300)]}

Now, we can program a simulation that runs each substep of the cooking while using the final temperature as input for the next simulation step.

def griddify(parsed_params, fine_time_grid):
    """Transforms the input temperature data to a regular time grid."""
    data = np.array(parsed_params['steps'])
    temps, t = data[:, :2], data[:, 2]
    t = np.cumsum(t)
    index = 0
    boundary_conditions = []
    for tt in fine_time_grid:
        if tt > t[index]:
            index += 1
            if index >= data.shape[0]:
                index -= 1
        boundary_conditions.append(temps[index])
    return boundary_conditions

def cook_recipe(recipe, n_snapshots=100, M=100, nx=100):
    """Takes a Cook my meat recipe as input, simulates it and returns temperature snapshots."""
    parsed_params = parse_input(recipe)
    total_time = sum([step_time for Tl, Tr, step_time in parsed_params['steps']])
    coarse_time_grid = np.linspace(0, total_time, num=n_snapshots)
    coarse_time_grid_indices = np.arange(n_snapshots) * M
    fine_time_grid = np.linspace(0, total_time, num=n_snapshots * M)
    dt = fine_time_grid[1] - fine_time_grid[0]
    boundary_conditions = griddify(parsed_params, fine_time_grid)
    x = np.linspace(0, parsed_params['thickness_mm'], num=nx)
    Tair = parsed_params['initial_temperature_C']
    T = np.ones_like(x) * Tair
    solver = CrankNicolsonSolver(x, T, dt, D=0.14)
    snapshots = []
    for ind, (t, (Tl, Tr)) in enumerate(zip(fine_time_grid, boundary_conditions)):
        solver.set_boundary_conditions(Tl, Tr, Tair)
        solver.step()
        if ind in coarse_time_grid_indices:
            T = solver.T.copy()
            snapshots.append(T)
    return x, coarse_time_grid, snapshots

x, time_grid, snapshots = cook_recipe(four_minutes_a_side, n_snapshots=50, M=500, nx=50)

hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), 
            kdims=['time', 'x'], 
            vdims='temperature').redim.range(temperature=(23, 100))

The above plot is qualitatively similar to the original app.

We can compare the final temperature from our solver with the final output of the app (obtained with some Javascript hacking):

x_app = np.linspace(0, 30, num=30)
final_temps_app = np.array([86.18586243519505,
  85.86449361301695,
  85.25487377408177,
  84.37978835482498,
  83.27157213371525,
  81.97050758592101,
  80.52286390721143,
  78.97869610729124,
  77.38953387284053,
  75.80609063774105,
  74.27611499608037,
  72.84249054685091,
  71.54166834143024,
  70.40249052953442,
  69.44543686657468,
  68.68229960704663,
  68.11626877794897,
  67.74239024927502,
  67.54834421381581,
  67.51548195134583,
  67.62005391346293,
  67.83456168692639,
  68.12916948444723,
  68.473116548245,
  68.83607928479354,
  69.1894401963763,
  69.50742898563637,
  69.76810900571606,
  69.954189125051,
  70.05364687780765])

hv.Curve((x_app, final_temps_app[::-1]), label='Cook my meat') * hv.Curve((x, snapshots[-1] - KELVIN), label='ours')

There seems to be a little offset between them.

I can think of several reasons for which the resultats are not in very good agreement:

• time and space grids are not the same in the two versions
• I used a method to project our system to a grid of times instead of applying a constant dt=1 throughout the simulation

However, the obtained temperatures still look quite close and the spatial shape is also quite resembling. So I'd say this validates our replication attempt.

Trying the different cooking methods

Let's now try the other default recipes featured in the app using our solver.

every_15_secs = """3cm Steak starts at 23°C
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15
150°C and 23°C for 0:15
23°C and 150°C for 0:15"""

x, time_grid, snapshots = cook_recipe(every_15_secs)

hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), 
            kdims=['time', 'x'], 
            vdims='temperature').redim.range(temperature=(23, 100))

sear_then_cook_low = """3cm Steak starts at 23°C
230°C and 23°C for 0:15
23°C and 230°C for 0:15
110°C and 110°C for 2:30
23°C and 23°C for 15:00"""

x, time_grid, snapshots = cook_recipe(sear_then_cook_low)

hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), 
            kdims=['time', 'x'], 
            vdims='temperature').redim.range(temperature=(23, 100))

sous_vide_liquid_nitrogen = """3cm Steak starts at 23°C
53°C and 53°C for 60:00
-200°C and -200°C for 0:30
200°C and 200°C for 2:00"""

x, time_grid, snapshots = cook_recipe(sous_vide_liquid_nitrogen)

hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), 
            kdims=['time', 'x'], 
            vdims='temperature').redim.range(temperature=(23, 100))

All of these cooking methods look quite similar to these produced by the Cook my meat app. Let's then turn to some analysis using our tool.

Post-processing to get the meat state

As in the app, let's first add some post-processing that allows to determine the meat quality instead of the temperature.

This essentially consists of monitoring the maximum temperature seen at a given location as a function of time.

This in turn determines the state of the proteins according to this table:

| Cooking state | Protein and temperature of change | |--------------- |--------------------------------------- | | Raw | No protein denaturization below 40°C | | Rare | Myosin denatures 40-55°C | | Medium Rare | Glycogen denatures 55-60°C | | Medium | Myoglobin denatures 60-70°C | | Well | Actin denatures 70-120°C | | Browned | Browning reactions 120-180°C | | Charred | Charring at >180°C |

We can implement this using the matplotlib color bar logic which already has a function for mapping from values to integers matplotlib.colors.BoundaryNorm.

from matplotlib import colors

bounds = [40, 50, 55, 60, 70, 120, 180]
norm = colors.BoundaryNorm(boundaries=bounds, ncolors=6, clip=False)

print(norm(39), norm(41), norm(51), norm(57), norm(61), norm(71), norm(121), norm(181))

PROTEIN_STATES = {-1: 'Raw', 
                  0: 'Rare',
                  1: 'Medium Rare',
                  2: 'Medium',
                  3: 'Well',
                  4: 'Browned',
                  5: 'Charred'}

def snapshots_to_protein_state(snapshots_K):
    snapshots_C = np.array(snapshots) - KELVIN
    protein_states = []
    for time_index in range(1, np.array(snapshots_C).shape[0] + 1):
        max_temps = np.max(snapshots_C[:time_index, :], axis=0)
        proteines_state = norm(max_temps)
        protein_states.append(proteines_state)
    return protein_states


x, time_grid, snapshots = cook_recipe(four_minutes_a_side, n_snapshots=50, nx=100)
hv.QuadMesh((time_grid, x, np.array(snapshots_to_protein_state(snapshots)).T),
            kdims=['time', 'x'], 
            vdims='protein_state', label='four minutes on each side').options(cmap='reds_r')

-1 0 1 2 3 4 5 6

layouts = []
for recipe, label in zip([every_15_secs, sear_then_cook_low, sous_vide_liquid_nitrogen],
                         ['flip every 15 seconds', 'sear then cook low', 'sous vide + liquid nitrogen']):
    x, time_grid, snapshots = cook_recipe(recipe)
    qm = hv.QuadMesh((time_grid, x, np.array(snapshots_to_protein_state(snapshots)).T),
            kdims=['time', 'x'], 
            vdims='protein_state', label=label).options(cmap='reds_r')
    layouts.append(qm)
hv.Layout(layouts).cols(1)

Comparison of cooking procedures

While discussing this post with colleagues, I found out that the solver we've developed is a nice tool to make hypotheses and understand the physics of cooking. In particular, I found two experiments interesting.

Why the liquid nitrogen matters

If we have another look at this recipe

3cm Steak starts at 23°C
53°C and 53°C for 60:00
-200°C and -200°C for 0:30
200°C and 200°C for 2:00

we could wonder if the liquid nitrogen is really needed. Does this really make a difference if we leave it out?

Let's see:

sous_vide_liquid_nitrogen_modified = """3cm Steak starts at 23°C
53°C and 53°C for 60:00
200°C and 200°C for 2:00"""

x, time_grid, snapshots = cook_recipe(sous_vide_liquid_nitrogen)
qm1 = hv.QuadMesh((time_grid, x, np.array(snapshots_to_protein_state(snapshots)).T),
            kdims=['time', 'x'], 
            vdims='protein_state', label='sous vide with liquid nitrogen').options(cmap='reds_r')

x, time_grid, snapshots = cook_recipe(sous_vide_liquid_nitrogen_modified)
qm2 = hv.QuadMesh((time_grid, x, np.array(snapshots_to_protein_state(snapshots)).T),
            kdims=['time', 'x'], 
            vdims='protein_state', label='modified sous vide without liquid nitrogen').options(cmap='reds_r')

hv.Layout([qm1, qm2]).cols(1)

(qm1[3750, :] + qm2[3750 - 30, :]).cols(1)

What we do see from the above comparison is that the cooling acts as a sort of "heat buffer" that prevents the searing at the end from adding too much heat to the steak center and thus overcooking much of the inner part. So in a sense, the liquid nitrogen method is quite optimal, the core being medium rare cooked and at the same having a "crust" which is slightly more cooked.

Several hours low-temperature roasts

A good example of a slow-roast recipe can be found here: https://www.thespruceeats.com/prime-rib-roast-slow-method-995284.

What is interesting is that the recipe has quite a lot of steps:

• searing the roast
• cooking at low temperature (200°F/93°C) until it reaches 128°F/53°C
• letting it sit until it reaches 120°C/48°C

I tried to transcribe the recipe using the Cook my meat format:

slow_roast = """12cm Steak starts at 23°C
200°C and 23°C for 3:00
23°C and 200°C for 3:00
93°C and 93°C for 100:00
23°C and 23°C for 25:00"""

x, time_grid, snapshots = cook_recipe(slow_roast)
qmt = hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), 
            kdims=['time', 'x'], 
            vdims='temperature', label='slow roast (temperature)').redim.range(temperature=(30, 70))

qmp = hv.QuadMesh((time_grid, x, np.array(snapshots_to_protein_state(snapshots)).T),
            kdims=['time', 'x'], 
            vdims='protein_state', label='slow roast (protein state)').options(cmap='reds_r')

hv.Layout([qmt, qmp]).cols(1)

As can be seen, the resting part is really important too since it allows the heat to change the proteins from rare to medium-rate.
Now the interesting part is to play with the recipe. Let's say we cook the steak 10°C higher. What happens?

slow_roast_10deg_higher = """12cm Steak starts at 23°C
200°C and 23°C for 3:00
23°C and 200°C for 3:00
103°C and 103°C for 100:00
23°C and 23°C for 25:00"""

x, time_grid, snapshots = cook_recipe(slow_roast_10deg_higher)
qmt = hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), 
            kdims=['time', 'x'], 
            vdims='temperature', label='slow roast (temperature)').redim.range(temperature=(30, 70))

qmp = hv.QuadMesh((time_grid, x, np.array(snapshots_to_protein_state(snapshots)).T),
            kdims=['time', 'x'], 
            vdims='protein_state', label='slow roast (protein state)').options(cmap='reds_r')

hv.Layout([qmt, qmp]).cols(1)

Too bad, the roast is now medium instead of medium rare!

An interesting application of our simulation tool is that we can use it to trace the right cooking time as a function of steak thickness. We will do this by a simple grid search.

from tqdm.notebook import tqdm

thicknesses = np.arange(8.5, 17, 0.5)
durations = np.linspace(30, 215, num=30)
results = []
for thickness in tqdm(thicknesses):
    for duration in durations:
        recipe = """{}cm Steak starts at 23°C
200°C and 23°C for 3:00
23°C and 200°C for 3:00
103°C and 103°C for {:.0f}:00
23°C and 23°C for 25:00""".format(thickness, duration)
        x, time_grid, snapshots = cook_recipe(recipe)
        protein_state = snapshots_to_protein_state(snapshots)[-1]
        val = np.sum(protein_state == 1) / protein_state.size
        results.append((thickness, duration, val * 100))

HBox(children=(FloatProgress(value=0.0, max=17.0), HTML(value='')))

hv.HeatMap(results, kdims=['thickness', 'duration'], vdims='percent_medium_rare').opts(tools=['hover'], width=500)

The heatmap is quite intuitive: there is a sweet spot in terms of cooking time where the roast becomes just right and if you pass it, it's overcooked. Below the sweet spot you get a little of the meat cooked the way you want it, but not the most of it. Also, an interesting observation from this plot is that the thicker your piece of meat, the wider the margin of error you have in terms of duration for cooking your roast right.

Conclusion

In this blog post, we went from implementing the heat equation with several boundary conditions to checking its results to doing some kitchen science. I hope you've enjoyed reading about this as much as I had pleasure developing it. Now back to cooking eggs :)

This post was entirely written using the Jupyter Notebook. Its content is BSD-licensed. You can see a static view or download this notebook with the help of nbviewer at 20200109_CookMyMeatClone.ipynb.

---

This article has been tagged with the following terms:

Food Numerical Physics Physics

---

• 
• 
• 
• 
• 

---

Comments