Towards a Mathematical and Numerical Model for Cooking Bratwurst

To make progress on this problem, we need first to develop a mathematical model of what's going on, and then apply a numerical method to approximate the solution to that mathematical model.

Mathematical Model

We have a chance at this if we are allowed to make several simplifying assumptions: Since we are assuming the brat is a perfect cylinder, we can reduce the problem from three space dimensions to one, by considering the temperature on a line from the center axis of the brat to the outer edge. Let u(x,t) represent the temperature at point x at time t. We will approximate the solution for x in the interval [0,R], where x=0 corresponds to the center point and x=R is the outer edge. And we will do this from time t=0 until the coldest spot in the meat (let's assume this is at x=0) achieves the desired minimum temperature Tdone.

Basic heat-flow modeling tells us that the following PDE must be true at each point in time and space:

    ut = alpha2 * uxx,

where the subscripts indicate partial derivatives, and alpha is the diffusion coefficient.

In order to find a unique solution u(x,t) to this problem we need to specify initial and boundary conditions. In this case, the initial condition is simply u(x,0) = Tinit. The boundary condition on the left endpoint is ux(0,t)=0. This is reasonable since, by symmetry, there should be no heat flow across this central point. The boundary condition on the right endpoint is derived by accounting for heat exchange between the meat and the outside world. Because the brat is rotating, the temperature of the `outside world' varies with time --- the temperature varying from the temperature of the grill to the temperature of the air above the grill, and back again. To be specific, a reasonable boundary condition on the right edge would be

ux(R,t) = C * (Text(t) - u(R,t)),

where C is a known constant and we'll take

Text(t) = Tave + Tdiff * cos(2*pi*t/Troll)

where

Tave = (Tgrill + Troom) / 2
Tdiff = (Tgrill - Troom) / 2

This leaves us with several more constants that need values. Below are some reasonable choices. These are scaled so that time is in minutes, distance is in centimeters, and temperature is in degrees C.

alpha = 1.236
Tinit = -4.0
Tgrill = 220.
Troom = 40.
Tdone = 90.
Troll = .2
C = 0.105 / 0.452 = (heat exchange coefficient) / (thermal conductivity)

Numerical Method

The Crank-Nicolson scheme described in Section 13.1 is almost exactly what we need. But this code (available in the on-line repository associated with the book) assumes that the boundary conditions are u(x,t) = 0 at the left and right endpoints; i.e., the unknown is known to be zero at both ends, for all time.

So what we need to do is modify the parabolic2 code to reflect our boundary conditions. This means to add two equations and two unknowns: the equations corresponding to our boundary conditions, and the unknowns corresponding to the temperature at the left and right endpoints. It also means that the equations corresponding to points 1 and n-1 need to be modified to include terms involving unknown values at point 0 (the left end) and point n (the right end), respectively.

We need some discrete equations to approximate the boundary conditions: