6.4: Three Dimensional Cake Baking
( \newcommand{\kernel}{\mathrm{null}\,}\)
In the rest of the chapter we will extend our studies to three dimensional problems. In this section we will solve the heat equation as we look at examples of baking cakes.
We consider cake batter, which is at room temperature of Ti=80∘F. It is placed into an oven, also at a fixed temperature, Tb=350∘F. For simplicity, we will assume that the thermal conductivity and cake density are constant. Of course, this is not quite true. However, it is an approximation which simplifies the model. We will consider two cases, one in which the cake is a rectangular solid, such as baking it in a 13′′×9′′×2′′ baking pan. The other case will lead to a cylindrical cake, such as you would obtain from a round cake pan.
This discussion of cake baking is adapted from R. Wilkinson’s thesis work. That in turn was inspired by work done by Dr. Olszewski,(2006) From baking a cake to solving the diffusion equation. American Journal of Physics 74(6).
Assuming that the heat constant k is indeed constant and the temperature is given by T(r,t), we begin with the heat equation in three dimensions,
∂T∂t=k∇2T.
We will need to specify initial and boundary conditions. Let Ti be the initial batter temperature, T(x,y,z,0)=Ti.
We choose the boundary conditions to be fixed at the oven temperature Tb. However, these boundary conditions are not homogeneous and would lead to problems when carrying out separation of variables. This is easily remedied by subtracting the oven temperature from all temperatures involved and defining u(r,t)=T(r,t)−Tb. The heat equation then becomes
∂u∂t=k∇2u
with initial condition
u(r,0)=Ti−Tb.
The boundary conditions are now homogeneous. We cannot be any more specific than this until we specify the geometry.
We will consider a rectangular cake with dimensions 0≤x≤W,0≤y≤L, and 0≤z≤H as show in Figure 6.4.1. For this problem, we seek solutions of the heat equation plus the conditions
u(x,y,z,0)=Ti−Tb,u(0,y,z,t)=u(W,y,z,t)=0,u(x,0,z,t)=u(x,L,z,t)=0,u(x,y,0,t)=u(x,y,H,t)=0.

Solution
Using the method of separation of variables, we seek solutions of the form
u(x,y,z,t)=X(x)Y(y)Z(z)G(t).
Substituting this form into the heat equation, we get
1kG′G=X′′X+Y′′Y+Z′′Z.
Setting these expressions equal to −λ, we get
1kG′G=−λ and X′′X+Y′′Y+Z′′Z=−λ
Therefore, the equation for G(t) is given by
G′+kλG=0.
We further have to separate out the functions of x,y, and z. We anticipate that the homogeneous boundary conditions will lead to oscillatory solutions in these variables. Therefore, we expect separation of variables will lead to the eigenvalue problems
X′′+μ2X=0,X(0)=X(W)=0,Y′′+v2Y=0,Y(0)=Y(L)=0,Z′′+κ2Z=0,Z(0)=Z(H)=0.
Noting that
X′′X=−μ2,Y′′Y=−v2,Z′′Z=−κ2,
we find from the heat equation that the separation constants are related,
λ2=μ2+v2+κ2.
We could have gotten to this point quicker by writing the first separated equation labeled with the separation constants as
1kG′G⏟−λ=X′′X⏟−μ+Y′′Y⏟−v+Z′′Z⏟−κ.
Then, we can read off the eigenvalues problems and determine that λ2=μ2+v2+ κ2.
From the boundary conditions, we get product solutions for u(x,y,z,t) in the form
umnℓ(x,y,z,t)=sinμmxsinvnysinκℓze−λmnkt,
for
λmnl=μ2m+v2n+κ2ℓ=(mπW)2+(nπL)2+(ℓπH)2,m,n,ℓ=1,2,…
The general solution is a linear combination of all of the product solutions, summed over three different indices,
u(x,y,z,t)=∞∑m=1∞∑n=1∞∑ℓ=1Amnlsinμmxsinvnysinκℓze−λmnkt,
where the Amn ’s are arbitrary constants.

We can use the initial condition u(x,y,z,0)=Ti−Tb to determine the Amnℓ's. We find
Ti−Tb=∞∑m=1∞∑n=1∞∑ℓ=1Amnlsinμmxsinvnysinkℓz.
This is a triple Fourier sine series.
We can determine these coefficients in a manner similar to how we handled double Fourier sine series earlier in the chapter. Defining
bm(y,z)=∞∑n=1∞∑ℓ=1Amnlsinvnysinkℓz,
we obtain a simple Fourier sine series:
Ti−Tb=∞∑m=1bm(y,z)sinμmx.
The Fourier coefficients can then be found as
bm(y,z)=2W∫W0(Ti−Tb)sinμmxdx.
Using the same technique for the remaining sine series and noting that Ti−Tb is constant, we can determine the general coefficients Amnl by carrying out the needed integrations:
Amnl=8WLH∫H0∫L0∫W0(Ti−Tb)sinμmxsinvnysinkℓzdxdydz=(Ti−Tb)8π3[cos(mπxW)m]W0[cos(nπyL)n]L0[cos(ℓπzH)ℓ]H0=(Ti−Tb)8π3[cosmπ−1m][cosnπ−1n][cosℓπ−1ℓ]=(Ti−Tb)8π3{0,for at least one m,n,ℓ even,[−2m][−2n][−2ℓ],for m,n,ℓ all odd.
Since only the odd multiples yield non-zero Amnℓ we let m=2m′−1, n=2n′−1, and ℓ=2ℓ′−1 for m′,n′,ℓ′=1,2,…. The expansion coefficients can now be written in the simpler form
Amn=64(Tb−Ti)(2m′−1)(2n′−1)(2ℓ′−1)π3.
Substituting this result into general solution and dropping the primes, we find
u(x,y,z,t)=64(Tb−Ti)π3∞∑m=1∞∑n=1∞∑ℓ=1sinμmxsinvnysinkℓze−λmntkt(2m−1)(2n−1)(2ℓ−1),
where
λmnℓ=((2m−1)πW)2+((2n−1)πL)2+((2ℓ−1)πH)2
for m,n,ℓ=1,2,…..
Recalling that the solution to the physical problem is
T(x,y,z,t)=u(x,y,z,t)+Tb,
we have the final solution is given by
T(x,y,z,t)=Tb+64(Tb−Ti)π3∞∑m=1∞∑n=1∞∑ℓ=1sinˆμmxsinˆvnysinˆκℓze−ˆλmnℓkt(2m−1)(2n−1)(2ℓ−1)

We show some temperature distributions in Figure 6.4.3. Since we cannot capture the entire cake, we show vertical slices such as depicted in Figure 6.4.2. Vertical slices are taken at the positions and times indicated for a 13′′×9′′×2′′ cake. Obviously, this is not accurate because the cake consistency is changing and this will affect the parameter k. A more realistic model would be to allow k=k(T(x,y,z,t)). However, such problems are beyond the simple methods described in this book.
In this case the geometry of the cake is cylindrical as show in Figure 6.4.4. Therefore, we need to express the boundary conditions and heat equation in cylindrical coordinates. Also, we will assume that the solution, u(r,z,t)=T(r,z,t)−Tb, is independent of θ due to axial symmetry. This gives the heat equation in θ independent cylindrical coordinates as
∂u∂t=k(1r∂∂r(r∂u∂r)+∂2u∂z2),
where 0≤r≤a and 0≤z≤Z. The initial condition is
u(r,z,0)=Ti−Tb,
and the homogeneous boundary conditions on the side, top, and bottom of the cake are
u(a,z,t)=0,u(r,0,t)=u(r,Z,t)=0.

Solution
Again, we seek solutions of the form u(r,z,t)=R(r)H(z)G(t). Separation of variables leads to
1kG′G⏟−λ=1rRddr(rR′)⏟−μ2+H′′H⏟−v2.
Here we have indicated the separation constants, which lead to three ordinary differential equations. These equations and the boundary conditions are
G′+kλG=0,ddr(rR′)+μ2rR=0,R(a)=0,R(0) is finite, H′′+v2H=0,H(0)=H(Z)=0.
We further note that the separation constants are related by λ=μ2+v2.
We can easily write down the solutions for G(t) and H(z),
G(t)=Ae−λkt
and
Hn(z)=sinnπzZ,n=1,2,3,…,
where v=nπZ. Recalling from the rectangular case that only odd terms arise in the Fourier sine series coefficients for the constant initial condition, we proceed by rewriting H(z) as
Hn(z)=sin(2n−1)πzZ,n=1,2,3,…
with v=(2n−1)πZ.
The radial equation can be written in the form
r2R′′+rR′+μ2r2R=0.
This is a Bessel equation of the first kind of order zero which we had seen in Section 5.5. Therefore, the general solution is a linear combination of Bessel functions of the first and second kind,
R(r)=c1J0(μr)+c2N0(μr).
Since R(r) is bounded at r=0 and N0(μr) is not well behaved at r=0, we set c2=0. Up to a constant factor, the solution becomes
R(r)=J0(μr).
The boundary condition R(a)=0 gives the eigenvalues as
μm=j0ma,m=1,2,3,…,
where j0m is the mth roots of the zeroth-order Bessel function, J0(j0m)=0.
Therefore, we have found the product solutions
Hn(z)Rm(r)G(t)=sin(2n−1)πzZJ0(raj0m)e−λnmkt,
where m=1,2,3,…,n=1,2,…. Combining the product solutions, the general solution is found as
u(r,z,t)=∞∑n=1∞∑m=1Anmsin(2n−1)πzZJ0(raj0m)e−λnmkt
with
λnm=((2n−1)πZ)2+(j0ma)2,
for n,m=1,2,3,….
Inserting the solution into the constant initial condition, we have
Ti−Tb=∞∑n=1∞∑m=1Anmsin(2n−1)πzZJ0(raj0m).
This is a double Fourier series but it involves a Fourier-Bessel expansion. Writing
bn(r)=∞∑m=1AnmJ0(raj0m),
the condition becomes
Ti−Tb=∞∑n=1bn(r)sin(2n−1)πzZ.
As seen previously, this is a Fourier sine series and the Fourier coefficients are given by
bn(r)=2Z∫Z0(Ti−Tb)sin(2n−1)πzZdz=2(Ti−Tb)Z[−Z(2n−1)πcos(2n−1)πzZ]Z0=4(Ti−Tb)(2n−1)π
We insert this result into the Fourier-Bessel series,
4(Ti−Tb)(2n−1)π=∞∑m=1AnmJ0(raj0m),
and recall from Section 5.5 that we can determine the Fourier coefficients Anm using the Fourier-Bessel series,
f(x)=∞∑n=1cnJp(jpnxa),
where the Fourier-Bessel coefficients are found as
cn=2a2[Jp+1(jpn)]2∫a0xf(x)Jp(jpnxa)dx.
Comparing these series expansions, we have
Anm=2a2J21(j0m)4(Ti−Tb)(2n−1)π∫a0J0(μmr)rdr.
In order to evaluate ∫a0J0(μmr)rdr, we let y=μmr and get
∫a0J0(μmr)rdr=∫μma0J0(y)yμmdyμm=1μ2m∫μma0J0(y)ydy=1μ2m∫μma0ddy(yJ1(y))dy=1μ2m(μma)J1(μma)=a2j0mJ1(j0m).
Here we have made use of the identity ddx(xJ1(x))=J0(x) from Section 5.5.
Substituting the result of this integral computation into the expression for Anm, we find
Anm=8(Ti−Tb)(2n−1)π1j0mJ1(j0m).
Substituting this result into the original expression for u(r,z,t), gives
u(r,z,t)=8(Ti−Tb)π∞∑n=1∞∑m=1sin(2n−1)πzZ(2n−1)J0(raj0m)e−λmmktj0mJ1(j0m).
Therefore, T(r,z,t) is found as
T(r,z,t)=Tb+8(Ti−Tb)π∞∑n=1∞∑m=1sin(2n−1)πzZ(2n−1)J0(raj0m)e−λnmktj0mJ1(j0m),
where
λnm=((2n−1)πZ)2+(j0ma)2,n,m=1,2,3,….
We have therefore found the general solution for the three-dimensional heat equation in cylindrical coordinates with constant diffusivity. Similar to the solutions shown in Figure 6.4.3 of the previous section, we show in Figure 6.4.6 the temperature evolution throughout a standard 9′′ round cake pan. These are vertical slices similar to what is depicted in Figure 6.4.5.


Again, one could generalize this example to considerations of other types of cakes with cylindrical symmetry. For example, there are muffins, Boston steamed bread which is steamed in tall cylindrical cans. One could also consider an annular pan, such as a bundt cake pan. In fact, such problems extend beyond baking cakes to possible heating molds in manufacturing.