% File: FE_cyl_conv_3_L2.m
%................................................................
%      FE Solution to Convection Outside a Thick Cylinder
%      r_1       r_2       r_3       r_4   convection
%  T_1 *---(1)---*---(2)---*---(3)---* [water, h_w, T_w]
%      T_1       T_2       T_3       T_4
%  Three Linear (L2) conduction line (radial) elements and 
%  one point (P1) convection element.
%................................................................
% Ref: Akin, "FEA w Error Estimators", 2005, section 8.2

%            Set constants
r_1 = 9.375  ; r_2 = 10.717  ; r_3 = 12.058 ; r_4 = 13.4 ; % radii
K = 1.736e-4 ; h_w = 2.89e-4 ; % conductivity, convection coeff
T_w = 45     ; T_1 = 400 ;     % water temp, inside temperature

% Eq (8.6) K/L for cylindrical conduction element (cancelling 2 pi)
S_1 = K * (r_2 + r_1) / (r_2 - r_1) / 2 ;   % conducting layer 1
S_2 = K * (r_3 + r_2) / (r_3 - r_2) / 2 ;   % conducting layer 2
S_3 = K * (r_4 + r_3) / (r_4 - r_3) / 2 ;   % conducting layer 3

% Elements 3 L2 1 P1: Manual assemble 4 by 4 square matrix terms
X = [r_1, r_2, r_3, r_4 ] ;                % merge lengths

% DOF: 1        2            3           4
S = [ S_1     -S_1          0,           0,            ; ...
     -S_1     (S_1 + S_2)   -S_2         0,            ; ...
      0,      -S_2          (S_2 + S_3)  -S_3          ; ...
      0,       0,           -S_3         (S_3 + r_4*h_w) ] ;

% DOF: 1  2   3   4 (convection)
C =   [0; 0;  0;  r_4*h_w*T_w] ;  % Assemble 4 by 1 source terms

% Save reaction data destroyed by following solver trick
Row_1_S = S(1:4) ; Row_1_C = C(1) ;

%  Apply essential BC at node 1, via trick to avoid partitions
C (:) = C (:) - T_1 * S (:, 1) ;  % carry known column to RHS
S (:, 1) = 0 ; S (1, :) = 0 ;     % clear, restore symmetry
S (1, 1) = 1 ; C (1) = T_1  ;     % trick not to change array size

T = S \ C  ;  % Solve system for four temperatures
% Running gives:  T = 400.0000  281.9159  177.8204  84.6259

% Recover reaction (heat flow in) from matrix system
q_in = ( Row_1_S * T - Row_1_C ) * 2 * pi ;  % here 0.964188 BTU/s
 
% Approximate element values = A(r)* q_r(r) = -k*A*dT/dr, A=2pi*r_cg
% q_1 = -K*pi*(r_1+r_2)*(T(2) - T(1)) / (r_2 - r_1) % radial heat flux 
% q_2 = -K*pi*(r_2+r_3)*(T(3) - T(2)) / (r_3 - r_2) % radial heat flux
% q_3 = -K*pi*(r_3+r_4)*(T(4) - T(3)) / (r_4 - r_3) % radial heat flux 

%...............................................................
%                     Add exact solution
r_a = 9.375; r_b = 13.4 ; % radii        (in)
T_a = 400  ; T_w = 45   ; % temperatures (F)
h_w = 2.89e-4           ; % convection   (BTU/s in^2 F)
k   = 1.739e-4         ;  % conductivity (BTU/s in F)

%               Analytic variable temperature
n  = 25               ;    % number of divisions
dr = (r_b - r_a) / n  ;    % plot increment
r = [r_a:dr:r_b]      ;    % radius of interest
Te = zeros (size(r))  ;    % pre-allocate array
Te = T_a - log (r/r_a)*(T_a - T_w) / (k/r_b/h_w + log (r_b/r_a)) ;

fprintf ('Exact surface temperature is__ %g F \n', Te(n+1)) % 84.6439 F
fprintf ('FE Surface temperature is ____ %g F \n', T(4))    % 84.6259 F 

%             Analytic constant heat flow, here 0.964627 BTU/s 
flow = 2*pi*(T_a - T_w) / (log (r_b/r_a) / k + 1 / (r_b * h_w)) ;
fprintf ('Exact heat flow is ___ %g BTU/s \n', flow)
fprintf ('FE heat flow is ______ %g BTU/s \n', q_in)
%...............................................................

%            Plot the temperature results
clf ; hold on ; grid on                    % clear & keep
title  ('Cylinder with internal temperature and external convection')
xlabel ('Radius (inches)') ; ylabel ('Temperature (F)')
plot (X, T, 'ko')                % 4 nodal values         
plot (r,  Te,    'r--')          % analytic radial
plot (X(4), T_w, 'k<')           % surface temperature
legend ('Three FE', 'Analytic', 'Surface')  % legend for 3 plots
plot (X, T, 'k-')                % 3 line elements
[f_ns] = num2str (flow) ;        % convert flow to string and list
text ((r_a+dr/2), ((T_a+T_w)/2.5), ['Exact heat flow = ', f_ns, ' (BTU/s)'])
% end FE_cyl_conv_3_L2.m

%    >> FE_cyl_conv_3_L2
%    Exact surface temperature is ___ 84.6439 F 
%    FE Surface temperature is ______ 84.6259 F 
%    Exact heat flow is ___ 0.964627 BTU/s 
%    FE heat flow is ______ 0.964188 BTU/s