Ch9.5 Heat Conduction Through a Wall

Introduction

  • 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.

title title

Problem Description

title

  • Given plate of width \( L \).
  • One side is hot and the other relatively cold.
  • Heat flows through material in \( x \)-direction towards cooler side.
  • A typical application is the wall of a building where inside and outside are at different temperatures.

General Compartmental Model

  • 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.

title title

General Compartmental Model

  • The only heat input is due to conduction in at \( x \).
  • The only heat output is due to conduction out at \( x+ \Delta x \).

title

title title

Input-Output Word Equation

title

\[ \begin{align*} \begin{Bmatrix} \mathrm{rate \, of \, change} \\ \mathrm{of \, heat ~in ~section } \end{Bmatrix} &= \begin{Bmatrix} \mathrm{rate \, heat } \\ \mathrm{conducted \, in~at~} x \end{Bmatrix} \\ \\ & - \begin{Bmatrix} \mathrm{rate ~ heat ~ conducted } \\ \mathrm{out~at~} x + \Delta x \end{Bmatrix} \end{align*} \\ \]

Variables and Parameters

title

  • To translate word equation into math equation, need to define our variables and parameters.
  • Heat flux is defined as the rate of flow of heat per unit area per unit time.
  • \( J(x) \) = flux of heat at \( x \)
  • \( U(x) \) = temperature at \( x \)
  • \( A \) = surface area of section

Input Word Equation

  • 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} \)

Output Word Equation

  • 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 \]

title title

Thermal Equilibrium

  • 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 \]

Formulating Differential Equation

  • 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 \]

Thermal Equilibrium ODE

title

  • Thermal equilibrium:

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

  • Solving for \( U(x) \),

\[ U(x) = a x + b \]

  • Time \( t \) is not present.
  • Equilibrium distribution of temperature \( U(x) \) as a function of \( x \).

Boundary Value Problem

title

  • Solution to ODE:

\[ U(x) = a x + b \]

  • Use thermometer to find
  • \( U(0)=u_0 \)
  • \( U(L)=u_L \).
  • Then
  • \( a = (u_L - u_0)/L \)
  • \( b = u_0 \)

Solving the ODE

  • 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.

Review of RK4 for 2nd Order IVP

  • 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*} \]

RK4 for Our IVP

  • 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).

Parameters for Numerical Solution

  • 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} \)

Numerical Solution Using R

Ch95Ex1(1,100)

title

Numerical Solution Using R

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)

Numerical Solution Using R

#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

Numerical Solution Using R

#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 }