Ch9.5 Heat Conduction Through a Wall

• We now formulate a differential equation for the simplest heat conduction problem.
• This involves the conduction of heat through a slab of some material in just one direction.
• We develop a differential equation for the equilibrium temperature reached as time increases.

• We consider a thin section, or slice, of material from \( x \) to \( x+ \Delta x \) as our compartment.
• The surface area of the face of this slice is \( A \).
• We must account for all inputs and outputs of heat for this section that contribute to rate of change of heat.

• Rate of heat conduction in at \( x \):

\[ \begin{Bmatrix} \mathrm{rate \, heat } \\ \mathrm{conducted \, in~at~} x \end{Bmatrix} = J(x)A \]

• Here's why:
• The units for rate of heat conducted in should be Joules/sec.
• Heat flux \( J(x) \) is rate of heat flow (Watts) per \( m^2 \):
• Unites of \( J(x)A = \left(\frac{W}{m^2}\right)(m^2) = W = \frac{Joules}{sec} \)

• Similarly for rate of heat leaving section, with

\[ \begin{Bmatrix} \mathrm{rate ~ heat ~ conducted } \\ \mathrm{out~at~} x + \Delta x \end{Bmatrix} = J(x+ \Delta x)A \]

• For equilibrium temperatures, LHS of word equation is zero:

\[ \begin{Bmatrix} \mathrm{rate \, of \, change} \\ \mathrm{of \, heat ~in ~section } \end{Bmatrix} = 0 \]

• Thus RHS also equals zero, with

\[ J(x)A - J(x+ \Delta x)A = 0 \]

• Divide both sides by \( (- A \Delta x) \) to obtain

\[ \frac{J(x+ \Delta x) - J(x)}{\Delta x} = 0 \]

• For thermal equilibrium,

\[ \frac{J(x+ \Delta x) - J(x)}{\Delta x} = 0 \]

• Letting \( \Delta x \) go to zero, we obtain

\[ \frac{dJ}{dx} = 0 \]

• Using Fourier's law, express in terms of temperature:

\[ \frac{d}{dx} \left(-k \frac{dU}{dx} \right) = 0 \]

• Our ODE is fairly simple:

\[ \frac{d^2 U}{dx^2} = 0 \]

• The analytic solution is found in Ch11.2, where it is viewed as a boundary value problem.
• We can solve numerically with MATLAB, Python, R, or Maple.
• More complicated differential equations can occur when we include volumetric heat sources, as is the case when heat is produced by an electric current or chemical reaction.

• Suppose we have a 2nd order IVP:

\[ y'' + ay' + cy = 0, ~y(0) = y_{10}, ~y'(0) = y_{20} \]

• Need first order slope functions for RK4.
• Second order IVP is split into system of first order ODEs
• Clever: Let \( y_1 = y, ~y_2 = y_1' \), and note \( y'_2 = y'' \).
• Then IVP above is equivalent to system below:

\[ \begin{align*} y_1' & = y_2,~y_1(0) = y_{10} \\ y_2' & = -ay_2-by_1, ~y_2(0) = y_{20}\\ \end{align*} \]

• Our IVP is

\[ y'' = 0, ~y(0) = y_{10}, ~y'(0) = y_{20} \]

• From Fourier's law, \( J = -ky' \)
• From previous slide, \( y_1 = y, ~y_2 = y_1' \). Then

\[ \begin{align*} y_1' &= y_2,~~ y_1(0) = y_{10} \\ y_2' &= 0, ~~ y_2(0) = y_{20} = (-1/k)J(0) \end{align*} \]

• Our goal is to find \( y_1 = y \), but we also find \( y_2 \) along the way.
• In this case, \( y_2 = (-1/k)J \), but generally it is a helping function (two equations, two unknowns).

• Our system of first order IVPs
  

\[ \begin{align*} y_1' &= y_2,~~ y_1(0) = y_{10} \\ y_2' &= 0, ~~ y_2(0) = y_{20} = (-1/k)J(0) \end{align*} \]

• The following values are inferred from the MATLAB and Maple code listings and surrounding descriptions in book.
• \( L = 1~ m \)
• \( k = 1 \frac{Joules}{m \cdot s \cdot C} \)
  
• \( U(0) = y_1(0) = y_{10} = 10 ~C \)
• \( U'(0) = y_2(0) = y_{20} = (-1/k)J(0) = -1 \frac{C}{m} \)

Ch95Ex1 <- function(L,n) {
#Ch9.5 Example
#Perform Rk4 for temperature y1 and
#scaled heat flux y2 = (-1/k)J(x)

#L is the width of slab or plate: [0, L]
#n is the number of time steps
h = L/n  #This is the time step size

#System Parameters
x0 <- 0
y10 <- 10  #initial value for y1 = y(0)
y20 <- -1  #initial value for y2 = y'(0)
k <- 1     #k value for conductivity (heat flux)

#System of ODEs
a <- 0
b <- 0
f1 <- function(y2) {y2}   
f2 <- function(y1,y2) {-a*y2-b*y1}   

#Initialize vectors
x <- rep(0, n)
y1 <- rep(0, n)
y2 <- rep(0, n)
x[1] <- x0
y1[1] <- y10
y2[1] <- y20

#Runge-Kutta Loop (Generate vectors x, y1 and y2)
for(i in 1:n) {
a1 <- h*f1(y2[i]);       #f1= slope of y1
a2 <- h*f2(y1[i],y2[i]); #f2= slope of y2
b1 <- h*f1(y2[i]+0.5*a1); #Half-step predictor
b2 <- h*f2(y1[i]+0.5*a2,y2[i]+0.5*a2);    
c1 <- h*f1(y2[i]+0.5*b1); #Half-step predictor
c2 <- h*f2(y1[i]+0.5*b2,y2[i]+0.5*b2);   
d1 <- h*f1(y2[i]+c1);     #Full-step predictor
d2 <- h*f2(y1[i]+c2,y2[i]+c2);     
y1[i+1] <- y1[i]+(1/6)*(a1+2*b1+2*c1+d1);  
y2[i+1] <- y2[i]+(1/6)*(a2+2*b2+2*c2+d2); 
x[i+1] <- x[i] + h }

• Rewrite 2nd order IVP as system of two first order IVPs.

\[ y'' + ay' + cy = 0, ~y(0) = y_{10}, ~y'(0) = y_{20} \]

• Always: Let \( y_1 = y, ~y_2 = y_1' \), and note \( y'_2 = y'' \).
• Then IVP above is equivalent to system below:

\[ \begin{align*} y_1' & = y_2,~y_1(0) = y_{10} \\ y_2' & = -ay_2-by_1, y_2(0) = y_{20}\\ \end{align*} \]

• Rewrite 2nd order IVP as system of two first order IVPs.

\[ y'' - y' - 2y = 0, ~y(0) = 4, ~y'(0) = 2 \]

• Let \( y_1 \) = _____ , \( y_2 \) = _____ , and note \( y'_2 = y'' \).
• Then IVP above is equivalent to system below:

\[ \begin{align*} y_1' & = \mathrm{Your~answer~here} \\ y_2' & = \mathrm{Your~answer~here} \\ \end{align*} \]

• Rewrite 2nd order IVP as system of two first order IVPs.

\[ y'' - y' - 2y = 0, ~y(0) = 4, ~y'(0) = 2 \]

• Let \( y_1 = y, ~y_2 = y_1' \), and note \( y'_2 = y'' \).
• Then IVP above is equivalent to system below:

\[ \begin{align*} y_1' & = y_2, ~y_1(0) = 4 \\ y_2' & = y_2 + 2y_1, ~y_2(0) = 2\\ \end{align*} \]

• Rewrite 2nd order IVP as system of two first order IVPs.

\[ y'' +3y' + 2y = 0, ~y(0) = 2, ~y'(0) = -3 \]

• Let \( y_1 \) = _____ , \( y_2 \) = _____ , and note \( y'_2 = y'' \).
• Then IVP above is equivalent to system below:

\[ \begin{align*} y_1' & = \mathrm{Your~answer~here} \\ y_2' & = \mathrm{Your~answer~here} \\ \end{align*} \]

• Rewrite 2nd order IVP as system of two first order IVPs.

\[ y'' +3y' + 2y = 0, ~y(0) = 2, ~y'(0) = -3 \]

• Let \( y_1 = y, ~y_2 = y_1' \), and note \( y'_2 = y'' \).
• Then IVP above is equivalent to system below:

\[ \begin{align*} y_1' & = y_2,~y_1(0) = 2 \\ y_2' & = -3y_2-2y_1, ~y_2(0) =-3\\ \end{align*} \]