% Conduction through three layered planar wall with convection.
% Ref: A.J. Chapman "Fundamentals of Heat Transfer", Macmillan, 1987
% Example 2.2, pages 46-47. 
%.......................................................................
% Outside (_o) convection at 35C to facing brick (_f) to brick (_b)
% to wall plaster (_p) to inside (_i) convection at 22C.  Note the 
% essential BC are hidden in the convection air temperatures (T_o, T_i).

%   convection    x_1       x_2       x_3       x_4   convection
%  [air, h_o, T_o] *---(1)---*---(2)---*---(3)---* [air, h_i, T_i]
%                 T_1       T_2       T_3       T_4
%                (P1)  (L2)      (L2)      (L2) (P1)            
% Connectivity list:    e     i      j
%                       1     1      2   % L2
%                       2     2      3   % L2
%                       3     3      4   % L2
%                       4     1      0   % P1
%                       5     4      0   % P1
%.......................................................................

%            Set given constants
L_f = 0.1  ; L_b = 0.15 ; L_p = 0.0125 ;  % conduction lengths
K_f = 1.32 ; K_b = 0.69 ; K_p = 0.48 ;    % thermal conductivities
T_o = 35   ; T_i = 22 ;                   % air temperatures (out, in)
h_o = 30   ; h_i = 8  ;                   % convection coefficients

%         Manually assemble 4 by 4 square matrix terms
%DOF: 1               2                3                  4
S=[(K_f/L_f+h_o), -K_f/L_f,           0,                 0,        ; ...
   -K_f/L_f,      (K_f/L_f+ K_b/L_b), -K_b/L_b,          0,        ; ...
    0,            -K_b/L_b,           (K_b/L_b+K_p/L_p), -K_p/L_p, ; ...
    0,             0,                 -K_p/L_p,          (K_p/L_p+h_i)]

C = [h_o*T_o; 0;  0;  h_i*T_i]  % Assemble 4 by 1 source terms
                                % Here EBCs are automatic (why)
T = S \ C                       % Solve for four temperatures (C)
 
q_f = -K_f * (T(2) - T(1)) / L_f % heat flux through wall (W/m^2)

%            Plot the temperature results
X = [0, L_f, (L_f+L_b), (L_f+L_b+L_p)] ;   % merge lengths
clf ; hold on ; grid on                    % clear & keep
dX = (X(4) - X(1)) / 10 ;
axis([(X(1)-dX) (X(4)+dX) T_i T_o])
plot (X, T, 'ko') ; plot (X(1), T_o, 'r<') % nodal and surface
legend ('Nodal values', 'Air surface')
plot (X, T, 'k-') ; plot (X(4), T_i, 'r<') % elements points
title ('Temperature Through a Three Material Wall')
xlabel ('Position (m)') ; ylabel ('Temperature (C)')

% S = [ 43.2000  -13.2000         0         0
%      -13.2000   17.8000   -4.6000         0
%             0   -4.6000   43.0000  -38.4000
%             0         0  -38.4000   46.4000 ]
% C = [ 1050 0 0 176]
% T = [ 34.0925 32.0301 26.1119 25.4030 ]
% q_f = 27.2238