Heat Conduction

1 Heat Conduction

1

2

HEAT CONDUCTION

T

1

2

3

4

5

6

7

8

air

9

10 T

air

L

Fig. 1.1 Computational grid for discretized continuum problem. The slab is of thickness L, we break the bar into 10 elements each of equal width. The ends of the bar are exposed to the ambient air. The elements 2-9 transfer heat by conduction with their nearest neighbor only. The two end elements transfer heat with their interior neighbor and the outside air. The slab extends to infinity in the vertical directions and into the page.

1.1 CONDUCTION IN A CONTINUUM In many cases we want to understand the thermal behavior of a body that cannot be characterized by a single temperature. We can use the concepts of the previous sections to develop a model for such a system. We return to our first example where we take a body out of the hot environment and place it in cool air, however now the body is large and the single temperature representation is poor. What do we do? To make our task simpler lets look at the problem in one-dimension. The “slab” is very large in two dimensions, with respect to its thickness. Both sides of the slab are exposed to the air and we want to model temperature gradients inside the slab. Instead of representing the body as one thermal mass, we will represent it be a collection of 10 coupled thermal masses. The system model is shown in Figure 1.1. As we have done in the previous sections, we can conduct an energy balance on each element. The amount of energy contained in any element is
 

 (1.1)

 where the variables are the same as we have used in previous sections. Weassume that the size, mass, and specific ) that flows from element to element heat of each element is constant. For an interior element, the heat flux (   is given by Fourier’s conduction law    !"  #  # %$ is the distance between the two points, where  is $ is the number of elements. Likewise the heat flux leaving the element i to    '& (  #  

#  is the thermal conductivity of the material and

# 

 12.  3&4"   *)+     /.  0!"  # #  5 +-,

(1.2)

where

(1.3)

Performing an energy balance on the interior element we obtain

(1.4)

where A is the #  cross-sectional area normal to the direction of heat flow. Since the mass of the element is dependent on its size, , we can rewrite the equation as

  /. 0 "  9/. '& (  #  .76/ ) +  #  #  5 (1.5) +8, 6 ; and 6/) specific heat do not change where is the material density. Assuming that the density, thermal conductivity, from element to element and using the = definition of thermal diffusivity as =:   (this property is tabulated for most materials, e.g coper = 1.17 < ?> and building insulation 0.001 < ?  > )) we obtain +   09@    '& 5 #  (1.6) +8, :

CONDUCTION IN A CONTINUUM

3

  are constant this equation is valid for any interior node. The property : is
 a material Since the material properties #  : . -> . Therefore we can define an #elemental   # heat conduction time constant property that has units of %$ we can write the conduction time constant Since the length total length by #  : of. the element is related to the   09@   '& as
  +

+8,



$


 

5

!"  /.  "   

  #  . 6 *)4+  7. #  5 +8,  "  which is rewritten as   "  5 +  

  $
   $
  +8,  6 *  ) # where
 . Likewise for the other side, node $ , the energy balance yields (    & "  5 +  

  $
   $ +-,
 

(1.7)

At the edges exposed to the air the energy balance becomes, for the first node for example,

1.1.1

(1.8)

(1.9)

(1.10)

Modeling and implementation

The governing differential equations for the system equations 1.7, 1.9, and 1.10. This system of ordinary differential equations can be solved using many of the techniques describes in the section on differential equations. Below we provide a sample program that will perform a simulation of the slab with the two ends exposed to the exterior. function cond1d() N = 20; %% number of points tau in = 1; %% tau in the interior tau out = 1; %% tau to the outside world y0 = ones(N,1); %% initial values options = []; [t,Y] = ode45(@eqns,[0 5],y0,options,N,tau in,tau out); %% Plotting the solution. for i =1:length(t) plot(Y(i,:)); axis([0 length(Y(1,:)) 0 1.1]); pause(0.1); end %% Governing differential equations function dy = eqns(t,y,N,tau in,tau out) dy = zeros(length(y),1); %% interior points dy(2:end-1) = y(1:end-2) - 2*y(2:end-1) + y(3:end); dy(2:end-1) = Nˆ2*dy(2:end-1)/tau in; %% boundaries dy(1) = -N*y(1) /tau out - Nˆ2*(y(1) -y(2) )/tau in; dy(end) = -N*y(end)/tau out - Nˆ2*(y(end)-y(end-1))/tau in;

The results of this model for different parameters are shown in Figure 1.2. In this figure we show the result for various time constants. The results can be understood quite intuitively. The time constant for convection controls  how long it takes for the boundary to take on the temperature of the air. The time constant of the interior,
 , determines how long it takes temperature changes to propagate across the bar. When the internal time constant is very short compared to the convective scale then we would expect that the interior of the bar can respond to changes at the boundary. Under these conditions we would expect that the bar would be at nearly uniform temperature. This fact is seen in Figure 1.2-iv. When the time scale of heat transfer to the air is very short compared to the time scale of the interior, then we would expect that the boundary temperature can change significantly, but the interior is slow to respond. This fact can be seen in 1.2-ii and 1.2-iii where significant thermal gradient develop in the bar. We can easily change the boundary conditions in the simulation. If we would like to maintain a fixed temperature on one wall, we simply set the temperature of that wall in the initial condition and change the boundary condition

Plots every 0.05 secs time

0.5

0 0 1

0.5 Length Plots every 0.05 secs time

0.5

0 0

1

0.5 Length

Plots every 0.05 secs time

0.5

0 0

1

Temp (normalized)

Temp (normalized)

1

Temp (normalized)

HEAT CONDUCTION

Temp (normalized)

4

1

0.5 Length

1

Plots every 0.05 secs time

0.5

0 0

1

0.5 Length

1

Fig. 1.2 Solution to the 1-D conduction problem for various parameters. Going left to right then down the parameters are: i) 
 ,   , ii) 
 ,   , iii) 
 ,    , and iv) 
 ,    . The time between plots of the temperature profile is noted as well as the direction of time. The bar always starts at a uniform temperature of one.

in the governing differential equations to be dy(1) = 0. If we would like to assume that the wall is adiabatic, that means that no heat can escape, then we simply set dy(1) = - (y(1) -y(2) )/tau in. 1.1.2

Continous equations



0 of the bar at one of the nodes, , in the simulation. The temperature at a neighboring Consider the temperature node to the right, , can be approximated by the Taylor Series.      +  #   +    # @   +  #   +!  # ! + +  +     +  ! "   +  #   +   +  + 

#  @

Likewise, the temperature to the left of node i may be written as

3&4  

+#  #   +!  # ! +     +  ! "  

  '& @   @+    # @   @4+!  # ! +  + ! "  0!9@   '&  +   @ +!  #   #  +  +! " 

(1.11)

(1.12)

Adding these two expressions together results in

which can be written as

(1.13)

(1.14)

Since the spacing between nodes is small we can assume that the term with the fourth derivative is negligible and 09@   3&4   the equation becomes

# 

$

+  +

(1.15)

TWO-DIMENSIONAL CONDUCTION

5

If we refer back to our numerical approximation to the conduction in the bar we find that the discrete equation is an  approximation to the continuous equation

,



:



 5

(1.16)

This result is the one-dimensional heat transfer equation which you will encounter in later courses.

1.2 TWO-DIMENSIONAL CONDUCTION The same arguments that applied to our derivation of the discrete one-dimensional heat conduction model apply equally well to the two-dimensions. In two dimensions we simply need to track the amount of heat that is conducted out of both the vertical and horizontal sides of a box. The program provided below is set to initialize the 2-D body at a uniform temperature of one and then hold the outer edge at zero. function T = cond2d() N = 20; %% number of points = 1; %% tau in the interior tau in tau out = .1; %% tau to the outside world y0 = ones(N,N); %% initial values y0(1,:) = zeros; %% set the temperature on the boundary to zero y0(end,:) = zeros; y0(:,1) = zeros; y0(:,end) = zeros; options = []; [t,Y] = ode45(@eqns,[0 1],y0(:),options,N,tau in,tau out); %% animated plot for i =1:10:length(t) T = Y(i,:); T = reshape(T,N,N); contourf(T,10) pause(0.01) end %% system of equations function dy = eqns(t,y,N,tau in,tau out) y = reshape(y,N,N); dy = zeros(size(y)); %% derivative w.r.t X and Y dy(2:end-1,2:end-1) = y(1:end-2,2:end-1) - 2*y(2:end-1,2:end-1) + y(3:end,2:end-1)... + (y(2:end-1,1:end-2) - 2*y(2:end-1,2:end-1) + y(2:end-1,3:end)); dy(2:end-1) = Nˆ2*dy(2:end-1)/tau in; %% boundary conditions dy(1,:) = 0; dy(end,:) = 0; dy(:,1) = 0; dy(:,end) = 0; dy = dy(:);

%% set the 2-d array into a 1-D array

6

HEAT CONDUCTION

Fig. 1.3 Two-dimensional temperature contours after 1 second of simulation with 

   and    . All figures have an initial condition of 1 and air temperature of 0. The figures going from left to right then down are i) all sides exposed to the air, ii) one side held at the initial temperature and all others exposed to the air, iii) adjacent sides held at the initial temperature, and iv) opposite sides held at the initial temperature.