Thermal finite element model

Geometrical considerations

The cross-sectional area for the bar is as follows.

a = pi/4*d^2; % m^2 x-sectional area

Dividing the bar into n sections (“finite elements”) such that we have length of each dx gives the following.

n = 100; % number of chunks
dx = L/n; % m ... length of chunk

Material considerations

The following are the material properties of copper.

cp = 390; % SI ... specific heat capacity
rho = 8920; % SI ... density
ks = 401; % SI ... thermal conductivity

Lumping

From the geometrical and material considerations above, we can develop a lumped thermal resistance R and thermal capacitance c of each cylindrical section of the bar of length dx. From and , these parameters are as follows.

R = dx/(ks*a); % thermal resistance
dV = dx*a; % m^3 ... section volume
dm = rho*dV; % kg ... section mass
c = dm*cp; % section volume

Linear graph model

The linear graph model is shown in with the corresponding normal tree overlayed.

State-space model

The state variables are clearly the temperatures of C_i: T_C1, ⋯, T_Cn. Therefore, the order of the system is n.

The state, input, and output variables are $$\begin{aligned} \bm{x} = \begin{bmatrix} T_{C_1} \cdots T_{C_n} \end{bmatrix}^\top, \quad \bm{u} = \begin{bmatrix} T_S \end{bmatrix},\text{ and} \quad \bm{y} = \bm{x}. \end{aligned}$$

Elemental, continuity, and compatibility equations

Consider the elemental, continuity, and compatibility equations, below, for the first, a middle, and the last elements. The following makes the assumption of homogeneity, which yields R_i = R and C_i = C.

Deriving the state equations for sections 1, i, and n

For each of the first, a representative middle, and the last elements, we can derive the state equation with relatively few substitutions, as follows. $$\begin{align} \dot{T}_{C_1} &= \frac{1}{C} q_{C_1} \tag{elemental} \\ &= \frac{1}{C} (q_{R_1} - q_{R_2}) \tag{continuity} \\ &= \frac{1}{RC} (T_{R_1} - T_{R_2}) \tag{elemental} \\ &= \frac{1}{RC} (T_S - T_{C_1} - T_{C_1} + T_{C_2}) \tag{compatibility} \\ &= \frac{1}{RC} (T_S - 2 T_{C_1} + T_{C_2}).\\ \dot{T}_{C_i} &= \frac{1}{C} q_{C_i} \tag{elemental} \\ &= \frac{1}{C} (q_{R_{i}} - q_{R_{i+1}}) \tag{continuity} \\ &= \frac{1}{RC} (T_{R_i} - T_{R_{i+1}}) \tag{elemental} \\ &= \frac{1}{RC} (T_{C_{i-1}} - 2 T_{C_i} + T_{C_{i+1}}). \tag{compatibility} \\ \dot{T}_{C_n} &= \frac{1}{C} q_{C_n} \tag{elemental} \\ &= \frac{1}{C} q_{R_{n}} \tag{continuity} \\ &= \frac{1}{RC} T_{R_n} \tag{elemental} \\ &= \frac{1}{RC} (T_{C_{n-1}} - T_{C_n}). \tag{compatibility} \end{align}$$

Let τ = RC. The A and B matrices are, then $$\begin{aligned} A &= \begin{bmatrix} -2/\tau & 1/\tau & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 1/\tau & -2/\tau & 1/\tau & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ & & & \ddots & \ddots & \ddots & & & & & \\ \vdots & & & & 1/\tau & -2/\tau & 1/\tau & & & & \vdots \\ & & & & & \ddots & \ddots & \ddots & & & \\ 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 1/\tau & -2/\tau & 1/\tau \\ 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 1/\tau & -1/\tau \end{bmatrix} \nonumber \\ B &= \begin{bmatrix} 1/\tau \\ 0 \\ \vdots \\ 0 \end{bmatrix}_{n\times 1}. \end{aligned}$$

The outputs are the states: y = x. Or, in standard form with identity matrix I, the matrices are: $$\begin{aligned} C = I_{n\times n} \quad\text{and}\quad D = 0_{n\times 1}. \end{aligned}$$

Simulation of a step response

Define the A matrix.

A = zeros(n);
% first row
A(1,1) = -2/(R*c);
A(1,2) = 1/(R*c);
% last row
A(n,n-1) = 1/(R*c);
A(n,n) = -1/(R*c);
% middle rows
for i = 2:(n-1)
  A(i,i-1) = 1/(R*c);
  A(i,i) = -2/(R*c);
  A(i,i+1) = 1/(R*c);
end

Now define B, C, and D.

B = zeros([n,1]);
B(1) = 1/(R*c);
C = eye(n);
D = zeros([n,1]);

Create a state-space model.

sys = ss(A,B,C,D);

Simulate a unit step in the input temperature.

Tmax = 1200; % sec ... final sim time
t = linspace(0,Tmax,100);
y = step(sys,t);

Plot the step response

To prepare for creating a 3D plot, we need to make a grid of points.

x = dx/2:dx:(L-dx/2);
[X,T] = meshgrid(x,t);

Now we’re ready to plot. The result is shown in .

figure
contourf(X,T,y)
shading(gca,'interp')
xlabel('x')
ylabel('time')
zlabel('temp (K)')