BackAcademic

Dec 11, 2024

Approximating Drink Cooling with PDEs

#Math

Have you ever wondered how long it takes for your favorite drink to cool down in the fridge? While intuition might suggest that smaller bottles cool faster or that turning up the fridge temperature helps, we can actually quantify the cooling process mathematically. This post delves into approximating the cooling time using partial differential equations (PDEs) to model heat transfer.

This article assumes you have a basic understanding of PDEs, including what they represent and how they describe changes in physical systems. If you're unfamiliar with these concepts, you might want to explore an introduction to PDEs first, as we'll be diving straight into applying the heat transfer equation in cylindrical coordinates.

The inspiration for this post comes from a recent course I took on numerical methods, where I discovered how elegant and fascinating such approximations can be. Here, we'll approximate the drink as a homogeneous, cylindrical object, focusing on radial heat transfer and ignoring the effects along the height (z-direction) for simplicity. Let's dive into the mathematics behind cooling your drink!

The Heat Transfer Equation

To model the cooling process, we use the heat equation, which describes how heat diffuses through a material. In cylindrical coordinates (), the equation can be written as:

Here, is the temperature, is time, is the radial distance, is the angular position, is the height, and is the thermal diffusivity of the material. For our simplified model, we assume:

The drink is radially symmetric, so,
The height is negligible, so,

With these assumptions, the equation simplifies to:

This reduced form focuses solely on radial heat transfer.

Our Specific Problem

To model the problem we want to solve, namely the time it will take our drink to cool down in the fridge, we need to define a few things:

The type of liquid
The initial and boundary conditions

The Liquid

Let’s say we want to figure out how long it will take a glass of water to cool down in the fridge. For this, we will need to define the thermal diffusivity, , for water.

Thermal diffusivity can be calculated with:

where, is the thermal conductivity, is the density, and is the specific heat capacity at constant pressure. We’ll use SI units for everything.

For water, these values are:

Therefore, the thermal diffusivity is:

💡For the smart cookies, technically, this diffusivity is only accurate at a specified temperature (25 ) because density changes with temperature. However, since density doesn’t change that much it is ok for our crude approximation.

Initial and Boundary Conditions

In order to solve the PDE, we need initial and boundary condtions:

Initial Condition: To start, the water is at room temperature,
Boundary Condition: Suddenly, when placed in the fridge, the edge of the cup () is at the fridge temperature,
Boundary Condition: At the centre of the cup, there is no heat flux, thus,

💡These are very heavily simplified conditions, in reality, the cup would never instantly be at the as soon as it is placed in the fridge, but I am just showcasing a very crude approximation.

Solving

Solving this PDE might be possible analytically, however, it is likely quite difficult. For our purposes, solving this numerically will be easier.

To solve this PDE numerically, we’ll implement a technique called the Method of Lines.

The Method of Lines (MOL) is a technique that transforms our PDE into a system of ordinary differential equations (ODEs) by discretizing the spatial derivatives while keeping time continuous. This allows us to use well-established ODE solvers to find our solution. Let's break down how we'll apply this method to our cooling problem.

Discretizing the spatial derivative
Creating our system of ODEs
Solving with established ODE solvers

Step 1: Discretizing

To discretize our spatial derivatives, we’ll apply something called Finite Differences.

💡Read more about finite differences here.

First, we'll divide our spatial domain (the radius of our cup) into N equally spaced points. At each point, we'll approximate our spatial derivatives using central differences. For the first derivative with respect to r, we use:

And for the second derivative:

Step 2: System of ODEs

From this, we can transform our original PDE into a discretized system of ODEs:

As you can see, each of the values have a subscript which means we can represent these as a matrix in code.

The represents the spacing between our values and is the value at that given index.

The represents our system of ODEs which we can solve to find our solution of values for any given time (in the specified solving time range).

Step 2.2: Resolving Issues

You may have noticed a problem with this when it comes to actually creating the system in code. We reference, in a few places. However, when , this point will be undefined. Similarly, when , the point will be undefined. Therefore, we need to handle these cases separately.

For the case, it is quite simple. We won’t solve for the last point because, if you recall, we defined one of the boundary conditions at as a constant value. The temperature at which occurs at will always be the same, .

For , we can recall our second boundary condition:

Therefore, at , our finite difference for the first derivative becomes:

From this, we can derive our expression for the ODE at

Step 3: Solving

We can use established ODE solvers to solve this system of ODEs. For the purpose of this article, I will use Python and solve_ivp function from SciPy.

For the specifics of this problem, let’s say the radius of our cup is

Implementation

I will break down the Python implementation into bite-sized parts, however, I will also provide the entire code file at the end of the post.

Constants and Imports

To start, let’s define our constants and imports:

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

k = 0.6 # W/m K, thermal conductivity
rho = 997 # kg/m3, density
cp = 4182 # J/kg K, specific heat capacity
alpha = k / (rho * cp) # thermal conductivity
L = 0.04 # m, radius of the cup
T_water = 298 # K, water start temperature
T_fridge = 277, # K, fridge temperature

N = 50 # number of spatial points we are solving for
r = np.linspace(0, L, N) # vector of r values
dr = r[1] - r[0] # spatial step size (Delta r)

Defining the System

The solve_ivp function takes three required arguments:

fun(t,y): a function that describes the RHS of the system. The time derivative of the state y, at time t where y can be an array (system).
t_span: a two member sequence (array, tuple, etc. containing two values) which represents the interval of integration (t0, tf). i.e. the start and end time for our solution
y0: The initial state of y. i.e. the initial condition.

To start, let’s define the system (fun(t,y)):

def dTdt(t, T):
    dT = np.zeros(T.shape)
    dT[0] = 2 * alpha * ((T[1] - T[0]) / dr ** 2)
    for i in range(1, N - 1): # iterating to N - 1 because we ignore the last point (constant)
        dT[i] = alpha * ((T[i + 1] - T[i - 1]) / (r[i] * 2 * dr)  + (T[i + 1] - 2 * T[i] + T[i - 1]) / (dr ** 2))

    return dT

We start by handling the special case of , after that we populate the rest of our ODE system. We only go until because, as mentioned previously, the temperature is fixed at the outer boundary, hence the change over time of the temperature, , remains zero.

Solving

Now, we are ready to solve. We need to define the start and end times for our solution. We’ll start by solving for one hour (). We also need to define our initial state. Our vector is all initially at except for the last value which is at :

t_span = (0, 60 * 60) # solving for 1 hour to start

T_initial = T_water * np.ones(N) # initially at 25 C
T_initial[-1] = T_fridge # r = L is at 4 C

From this, we can call solve_ivp

sol = solve_ivp(
        dTdt, 
        t_span, 
        T_initial, 
        t_eval=np.linspace(t_span[0], t_span[1], 60),
        method="BDF"
    )

💡The t_eval parameter is used to explicitly define how many time points we want to solve for. If it is not provided, it is chosen through the solve_ivp algorithm. We define it to be 60 equally spaced points from 0 to 3600 (1 hour). One for each minute.

💡The method parameter is used to set the algorithm used for solving the system of ODEs. For this problem, I chose the BDF (backward differentiation) method which is an implicit method, providing more stable results. Using the default solver, RK45 (Runge-Kutta 4th order) would cause the solution to display some instability.

Visualizing Results

Finally, we need to actually find the answer to our problem. We can do this through visualizing the results. We’ll plot the temperature as a function of radius for various times. To start, we’ll plot the distribution every 10 minutes to see how the temperature is changing.

for i in range(0, len(sol.t), 10): # iterating over time points by 10
    # plotting temperature (converted to C) vs. radius at time, i, 
    plt.plot(r, sol.y[:, i] - 273, label=f"t = {int(sol.t[i] / 60)} min")

plt.legend()
plt.grid(True)
plt.xlabel("Radius (m)")
plt.ylabel(r"Temperature ($^\circ C$)")
plt.show()

Generated plot of temperature vs. radius at 10 minute intervals for 1 hour.

💡I did some customization to make my plot colours a little different, however, I did not include this as it is irrelevant.

As we can see, the temperature is gradually getting closer and closer to the the fridge temperature. I will consider the drink cooled down when the centre of the drink has also reached . Currently, even after 1 hour, the centre of the drink is still at around . This means we need to leave it in the fridge for longer. Let’s run the solver for 2 hours now and see what happens:

t_span = (0, 120 * 60) # solving for 2 hours now

💡Leaving our code as is except for this change will lead to the plots being generated at intervals of 20 minutes now. I will leave it as an exercise to figure out how to make it 10 minute intervals.

Generated plot of temperature vs. radius at 20 minute intervals for 2 hours.

At 2 hours (120 min), we can see that the temperature distribution is a flat line around . This means that it takes 2 hours for the entire drink to fully cool down in the fridge. Pretty neat, right?

Animation

For fun, I also created an animation showing how the temperature changes over time. I won’t include the code for this as it would make this article unnecessarily long, however, I think it’s pretty cool.

Temperature vs. radius animated over 2 hours at intervals of 20 minutes.

Conclusion

In this article, we explored how to use partial differential equations and numerical methods to model the cooling of a drink in a fridge. Through our implementation in Python, we were able to determine that it takes approximately 2 hours for a room temperature drink to fully cool down to fridge temperature. This mathematical approach not only satisfies our curiosity but also demonstrates the power of PDEs in solving real-world problems. So, next time you come across someone saying “when am I going to use this complex math in real life?”, you can show them this article!

Full Code

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

plt.rcParams.update({
  "text.usetex": True, # latex rendering for plots
})

k = 0.6 # W/m K, thermal conductivity
rho = 997 # kg/m3, density
cp = 4182 # J/kg K, specific heat capacity
alpha = k / (rho * cp) # thermal conductivity
L = 0.04 # m, radius of the cup
T_water = 298 # K, 25 C water start temperature
T_fridge = 277 # K, 4 C fridge temperature

N = 50 # number of spatial points we are solving for
r = np.linspace(0, L, N) # vector of r values
dr = r[1] - r[0] # spatial step size (Delta r)

def dTdt(t, T):
    dT = np.zeros(T.shape)
    dT[0] = 2 * alpha * ((T[1] - T[0]) / dr ** 2)
    for i in range(1, N - 1): # iterating to N - 1 because we ignore the last point (constant)
        dT[i] = alpha * ((T[i + 1] - T[i - 1]) / (r[i] * 2 * dr)  + (T[i + 1] - 2 * T[i] + T[i - 1]) / (dr ** 2))

    return dT

t_span = (0, 120 * 60) # solving for 2 hours now

T_initial = T_water * np.ones(N) # initially at 25 C
T_initial[-1] = T_fridge # r = L is at 4 C

sol = solve_ivp(
        dTdt, 
        t_span, 
        T_initial, 
        t_eval=np.linspace(t_span[0], t_span[1], 61),
        method="BDF"
    )

for i in range(0, len(sol.t), 10): # iterating over time points by 10
    # plotting temperature vs. radius at time, i, 
    plt.plot(r, sol.y[:, i] - 273, label=f"t = {int(sol.t[i] / 60)} min")

legend = plt.legend()
legend.get_frame().set_facecolor("black")
plt.grid(True)
plt.xlabel("Radius (m)")
plt.ylabel(r"Temperature ($^\circ C$)")
plt.show()