Dynamic Thermal Systems

MICTR1Heat transfer – Dynamic models1

02.03.2011/SGt

“State space”thermal models

1.Introduction

If we are looking on thermal systems in which the heat capacities are large or the time in consideration is short, we can not assume the system to be independent of time, i.e. there will be a rate of change in temperatures in the system during the process.

Let us look at an example to illustrate, what we are talking about. Figure 1 shows a thick wall made of a material with a low coefficient of heat conductivity, let’ssay rubber. At a given time the temperature on its left side is increased from tinit to ta.

Figure 1: Dynamic heat flow in a flat plate.

At beginning (τ0) the temperature in the plate is the same all the way through the plate, tinit. Later (τ1 and τ2), the temperature raises against the ambient temperature, ta.

At initial time, τ0, it has an overall initial temperature of let say 0°C. It is isolated on its sides but at the one end it is put into contact with hot water at a temperature on 100°C. After a few seconds, the temperature of the plate is still 0°C, except at the one end close to the water, where it heated up a little. After some more time the heat will propagate through the plate as indicated by the different time history curves, see curves for τ1 and τ2.

The actual modeling equations for the system has been develop by Fourier and is given by the second order partial differential equation in space and time

(1)

Here  is the “coef. of thermal conductivity” given by  = k /(cp).

This equation can be solved by – for some situations – analytical methods, but for most practical systems only by numerical methods, like FEM.

2.MultipleLumped Capacitance Modeling

Here we will look at another method, called the lumped parameter model. What we do is to divide the plate into individual elements or lumps in which the temperature is – not in time – considered to the constant. By this method we transform the partial differential equation (6) into a number of ordinary differential equations (ODE’s).

The more number of lumps the higher accuracy, but also more time consuming is the solving process! What is the necessary minimum number of elements?The problem has been solved by Biot. He found out, that it is a question of the relative proportion between the resistancies in the system: Let us explain it by an example: We have a long solid rod made by bad conduction material like rubber, and you put this rod into contact with a hot liquid like water! Then it is reasonably, that the temperature in the far end – away from the hot water - will not begin to raise before after some time. In this system you have a low convective heat resistance Rconv compared to the resistance to heat conduction, Rcond. Then you should imagine another system. You take a short solid rod made of silver (Silver: k = 420 W/(m K)) and put this into contact with hot air at the one end. In this situation you will experience a raising temperature of the silver rod - as before - but you will not see any temperature difference in the rod – it will have (almost) the same temperature in both ends at all times.

Biot defined a ratio called Biot’s number:

[-](2)

where Lc is a characteristic length of the solid material, given by volume divide by surface area, i.e. Lc = V/A . We have:

(3)

If Biot’s number is small (NB < 0.1) then there will be a small temperature difference in the material from one side to the other as with the silver rod. In this case we can consider all the heat energy stored in the rod as a single lump of heat capacity.

If Biot’s number is large (NB > 0.1) we must consider the heat capacity divided into several lumps with individual temperatures in each of them. We often talk about nodes where we say, that the whole heat capacity for each of the element is located in the node.

Example:

Given a thick flat plate as shown on figure 4 we first have to estimate the minimum number of elements in which the plate should be divided.

If we divide the length – or thickness – into n elements we have for the heat resistance between any of them

(4)

The minimum number of elements is the given by the equation NB0.1 or

(5)

or

(6)

The heat capacity of the elements is

(7)

WhereM is the total mass of the plate and cp is the specific heat capacity of the material

Figure 2 shows a model of the heat flow through a thick plane wall. The wall is divided into 5 elements and the temperatures are calculated in 5 notes.

Note-1:

For the temperaturet1, note-1, we have for the heat flows in and out of the element:

In: where (8)

Out: where (14)

The conservation of energy gives:

(9)

which gives

(10)

And then

(11)

Figure2: Multiple lumped model (5 elements) of the thick wall. The arrow (top) shows the second element with the temperature t2

Note-2:

where (12)

Energy balance

(13)

And from here

(14)

Note j = 3 and j = 4 are semilar to note 2.

hvor j=3,4 (15)

Note-5:

In: (16)

Out: where (17)

Conservation of energy gives.

(18)

And then

(19)

3.State space modelling

A special way of presenting the equations is called ”state space”. On vector form it looks like this

The “bold” symbols indicates arrays or matrices

Her x = t, is a column vector representing the temperaturesn in the five notes. uis a two dimensional vector representing the ambient temperaturesti og tu. Theyvector is an ”output” vector – here the 5 temperatures1.

Using the example above, with the five notes, we get:

(20)

(21)

4.Matlab/Simulink

In Simulink a special block is available for this kind of problem. Figure 3 show the model and the content of the State-Space block.

The script file running the model is shown in chapter 5.

Figure 3a: Simulink model /
Figure 3b: State-Space-block

Running the model gives this plot:

Figure 4: Temperature distribution in the wall (see data in the script file)

5.Matlab code

MICTR1Heat transfer – Dynamic models1

02.03.2011/SGt

% Thick wall

% 28.02.2011/SGt

clear

clc

close all

N = 5; % Number of notes

A = 1; % [m2] Area

alfa_i = 10; % [W/(m2K)] Coef. of heat transfer, internal wall

alfa_o = 10; % [W/(m2K)] Coef. of heat transfer, outer wall

L = 0.10; % [m] Thickness of wall

k = 0.20; % [W/(m*K)] Thermal conductivity

rho = 500; % [kg/m3] Density

c_p = 2000; % [J/(kg*K] Spec heat capacity

t_i = 20; % [°C] Temperture, fluid inside

t_o = 0; % [°C] Temperature, fluid outside

t_init = 0; % [°C] Initial temperature

% Some derivatives:

Biot_i = alfa_i/(k/L) % [-] Biot number, internal

Biot_o = alfa_o/(k/L) % [-] Biot number, external

R_k = (L/(N-1))/(k*A) % [K/W] Thermal resistance, inner wall

R_i = 1/(alfa_i*A) % [K/W] Thermal resistance, inside wall

R_o = 1/(alfa_o*A) % [K/W] Thermal resistance, out side wall

M = A*L*rho % [kg] Mass of wall

C_i = c_p*M/(N-1) % [J/K] Thermal capacity, wall

U = 1/(1/alfa_i + L/k + 1/alfa_o) % [W/(m2K)] - Coef. of heat transfer

% Steady-state conditions:

t_wi_2 = t_i - U/alfa_i*(t_i-t_o)

t_wu_2 = t_o + U/alfa_o*(t_i-t_o)

PHI_i_1 = alfa_i * A * (t_i -t_init)

PHI_o_1 = alfa_o * A * (t_init-t_o)

PHI_i_2 = alfa_i * A * (t_i-t_wi_2)

PHI_o_2 = alfa_o * A * (t_wu_2-t_o)

% Matricerne i state space modellen beregnes

% Matrice A

for i = 1 : N;

for j = 1 : N;

A_w(i,j) = 0;

if i == j;

A_w(i,j) = -2/(C_i*R_k);

elseif j==i-1;

A_w(i,j) = 1/(C_i*R_k);

elseif j==i+1;

A_w(i,j) = 1/(C_i*R_k);

end

end

end

A_w(1,1) = -(1+R_k/R_i)*2/(C_i*R_k);

A_w(1,2) = 1*2/(C_i*R_k);

A_w(N,N-1) = 1*2/(C_i*R_k);

A_w(N,N) = -(1+R_k/R_o)*2/(C_i*R_k);

% Matrice B

for i = 1 : N;

for j = 1:2;

B_w(i,j) = 0;

end

end

B_w(1,1) = 2/(C_i*R_i);

B_w (N,2) = 2/(C_i*R_o);

% Matrice C

for i = 1 : N;

for j = 1 : N;

C_w(i,j) = 0;

if i == j;

C_w(i,j) = 1;

end

end

end

% Matrice D

for i = 1 : N;

for j = 1 : 2;

D_w(i,j) = 0;

end

end

A_w

B_w

C_w

D_w

simtime=3600*1;

[tau,t]=sim('sim_wall_3',[0 simtime]);

lw = 2; % Linewidth

figure(1)

plot(tau/60,t)

grid

AXIS([0 simtime/60 0 inf])

ylabel('t [°C]')

title(['N = ' num2str(N) ' Biot_i = ' num2str(Biot_i)...

' Biot_o = ' num2str(Biot_o)]);

xlabel('Time, tau [min]')

figure(2)

PHI_i=alfa_i*A*(t_i-t(:,1));

SUBPLOT(3,1,1)

plot(tau/60,PHI_i)

grid

AXIS([0 simtime/60 -inf inf])

ylabel('PHI_i [W]')

PHI_o=alfa_o*A*(t(:,N)-t_o);

SUBPLOT(3,1,2)

plot(tau/60,PHI_o)

grid

AXIS([0 simtime/60 0 inf])

ylabel('PHI_o [W]')

dif=PHI_o-PHI_i;

SUBPLOT(3,1,3)

plot(tau/60,dif)

grid

% AXIS([0 simtime/60 0 inf])

ylabel('PHI_o - PHI_i [W]')

xlabel('time [min]')

figure(3)

x=0:L/(N-1):L

n_pl = 5;

n_tau = length(tau);

n(1)=1;

for i = 2:n_pl

n(i) = round(n_tau*((i-1)/(n_pl-1)));

end

for i = 1:n_pl

hold on

h1=plot(x,t(n(i),:),'k-');

set(h1,'linewidth',lw);

hold off

end

hold on

title('Temperature vs. distance and time')

xlabel('x [m]')

ylabel('t [°C]');

grid

dif_t = 1;

AXIS([0 L min(min(t_i,t_o),t_init)-dif_t max(max(t_i,t_o),t_init)+dif_t])

hold on

t_ix = ones(length(x),1)*t_i;

h1= plot(x,t_ix,'r-');

set(h1,'linewidth',lw);

hold off

hold on

t_ox = ones(length(x),1)*t_o;

h1= plot(x,t_ox,'b-');

set(h1,'linewidth',lw);

hold off