% Author: Layla S. Mayboudi
% 2021-05-27
% This code creates a transient thermal PDE model for an externally-finned pipe.
% It plots the geometry and heat flux flow diagrams.
% Selected nodes and their temperatures are also displayed.

% Note:
% A transient model used ThemalConductivity, MassDensity and SpecificHeat. 

clear all; 
close all; clf, clc;
set(groot,'DefaultAxesColorOrder',[0 0 0],...
    'DefaultSurfaceEdgeColor',[0 0 0],...
    'DefaultLineLineWidth',1,...
    'DefaultAxesLineStyleOrder','-|--|:|-.|+|o|*|.|x|s|d|^|v|>|<|p|h');      %may be added to set the line style and color.
ax.FontSize = 20; ax.FontName = 'Cambria'; ax.FontSmoothing = 'on'; ax.FontWeight = 'normal'; %'bold'; 

'DefaultAxesLineStyleOrder','-|--|:|-.|+|o|*|.|x|s|d|^|v|>|<|p|h'      %may be added to set the line style and color.
set(gcf,'defaultSurfaceEdgeColor','b')
  
thermalModelT = createpde('thermal','transient-axisymmetric'); % create a model object: thermal

%radius_P = 0.875*25.4/2000; offset = 0.681*25.4/2000;         %pipe diameter 3/4" %m
radius_P = 1.125*25.4/2000; 
offset = 0.995*25.4/2000;  %radius_E = insulation_th+radius_P; %pipe diameter 1" %m
radius_O = offset;
thickness = radius_P - offset;

%insulation_th = 0.5*25.4/1000; %m
length_P = 25/1000; %m
length_I = 1/2000; %m
height_I = offset*1.5; %m
radius_E = height_I+radius_P;

geom_P = [3;4;offset;offset;radius_P;radius_P; 0;length_P;length_P;0]; % create geometry: rectangle solid

% create pattern of rectangles (fins)
n_pattern = 11; % number of rectangles in pattern
dz1 = (length_P-n_pattern*length_I)/(n_pattern-1)
geom_I2 = zeros(length(geom_P),n_pattern);
dz = dz1;
sf = 'geom_P';
for i = 1:n_pattern
    space = dz*(i-1);
    side_1 = space + (i-1)*length_I;
    geom_I2(:,i) = [3;4;radius_P;radius_P;radius_E;radius_E; side_1;length_I+side_1;length_I+side_1;side_1]; % create geometry: rectangle solid
    geom_names(i)=string(['geom_I2_' num2str(i)]);
    sf=[sf '+' geom_names{i}];

end
gd = [geom_P,geom_I2];
ns=char('geom_P',geom_names{1:n_pattern});
ns = ns';
dl = decsg(gd,sf,ns);
geometryFromEdges(thermalModelT,dl); % associate geometry and thermal model object
figure(1);
%pdegplot(thermalModelT);
pdegplot(thermalModelT,'EdgeLabels','on','FaceLabels','on');
ylim([0-1/1000 length_P+1/1000]); xlim([0 radius_E+offset+3/1000]);
title({'Geometry with edge and face labels displayed';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14;

% Pipe material properties - transient thermal analysis (copper, Cu)
k1 = 400; % copper thermal conductivity (W/mK)
k_Cu = k1;
cp1 = 385; % copper heatcapacity (J/kgK)
rho1 = 8960 % copper density (kg/m^3)
q1 = 0; %heat generation (W/m^3)
tp1=thermalProperties(thermalModelT,'ThermalConductivity',k1,'MassDensity',rho1,'SpecificHeat',cp1,'Face',1) 

% Fin material properties - transient thermal analysis (aluminium, Al)
% http://www.matweb.com/search/datasheet.aspx?bassnum=AMEAL00&ckck=1
k2 = 238; % copper thermal conductivity (W/mK) 
k_Al = k2;
cp2 = 900; % copper heatcapacity (J/kgK) 
rho2 = 2700 % copper density (kg/m^3) 
q2 = 0; %heat generation (W/m^3)
tp2=thermalProperties(thermalModelT,'ThermalConductivity',k2,'MassDensity',rho2,'SpecificHeat',cp2,'Face',[2:12]) 
internal_heat = internalHeatSource(thermalModelT,q1,'Face',2);

% Boundary conditions
T_ambient = 25;  % degC = 75 degF
T_water = 80;   % degC = 180 degF   %or T_water = 60;
hc = 10; % degC 180 degF % convection heat transfer coefficient W/m^2K
q = 0;
thermalBC(thermalModelT,'Edge',[1:57],'ConvectionCoefficient',hc,...
    'AmbientTemperature',T_ambient); % convection heat transfer
thermalBC(thermalModelT,'Edge',[22 23 45 46],'HeatFlux',q); % convection heat transfer
thermalBC(thermalModelT,'Edge',1,'Temperature',T_water); % constant temperature

% Initial conditions
t_f = 0.1; % hours
tfinal = t_f * 60 *60; % set final time (s) for 1 hour = 3600 s
time_step = 0.25; % time step (s)
tlist = 0.1:time_step:tfinal; % set time array
time_l = length(tlist);
tt = round(tlist(time_l)/60,1) % time in minutes
ttt = round(tlist(time_l)/60,1);
t2 = tlist/60; % time in minutes
T0 = T_ambient;
ic = thermalIC(thermalModelT,T0);

bc = findThermalBC(thermalModelT.BoundaryConditions,'Edge',[1 2])
ic = findThermalIC(thermalModelT.InitialConditions,'Face',[1 2])

% Plot geometry
figure(2);
grid_size = length_P/100;
mesh1 = generateMesh(thermalModelT,'Hmax',grid_size); % generate mesh *****
pdeplot(mesh1); % plot PDE thermal model object with node ID's displayed
ylim([-1/1000 length_P+1/1000]); xlim([0 radius_E+2/1000]);
title({'Triangular elements';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14;

% Plot node numbers
figure(3);
pdeplot(mesh1,'NodeLabels','on'); % plot PDE thermal model object with node ID's displayed
ylim([0 2*grid_size]); xlim([offset offset+2*grid_size]); % set axis conditions
title('','FontWeight','normal');
title({'Node ID''s';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14;

% Plot element numbers
figure(4);
pdeplot(mesh1,'ElementLabels','on'); % plot PDE thermal model object with node ID's displayed
ylim([0 2*grid_size]); xlim([offset offset+2*grid_size]); % set axis conditions
title('','FontWeight','normal');
title({'Elements ID''s';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14;

% Results
thermalModelT.SolverOptions.ReportStatistics = 'on'
resultT = solve(thermalModelT,tlist) % identify PDESolverOptions Properties
T1 = resultT.Temperature; % identify data to process and plot *****
Tx1 = resultT.XGradients; % identify data gradient along x-coordinate *****
Ty1 = resultT.YGradients; % identify data gradient along y-coordinate *****
T_st1 = resultT.SolutionTimes; % identify data gradient along y-coordinate *****

% Contour plot thermal model steady temperature versus rz
figure(5);
pdeplot(thermalModelT,'XYData',T1(:,time_l),'colormap','hot','contour','on','levels',10); % 'colormap','hot' ,'Contour','on'
ylim([0 length_P]); xlim([offset-2*grid_size radius_E+2*grid_size]);
title({'Temperature contours (degC)'; ['(t = ',...
    num2str(tt),' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14;

% Contour plot temperature gradient (dT/dr) versus rz
figure(6);
pdeplot(thermalModelT,'XYData',Tx1(:,time_l),'colorMap','hot','contour','on','levels',10); 
ylim([0 length_P]); xlim([offset-2*grid_size radius_E+2*grid_size]);
title({'Temperature gradient contours (dT/dr, degC/m)'; ['(t = ',...
    num2str(tt),' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14;

% Contour plot temperature gradient (dT/dz) versus rz
figure(7);
pdeplot(thermalModelT,'XYData',Ty1(:,time_l),'colormap','hot','contour','on','levels',10); 
ylim([0 length_P]); xlim([offset-2*grid_size radius_E+2*grid_size]);
title({'Temperature gradient contours (dT/dz, degC/m)'; ['(t = ',...
    num2str(tt),' min)'];''},'FontWeight','normal')
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14;
 
%**********************************************************
%Investigate the node density

point = thermalModelT.Mesh.Nodes; % recall node locations from thermal model object
nodesIndex = find(point);
xy00=point; r00 = xy00(1,:); z00 = xy00(2,:);
radius_sorted = sort(r00); index_sort = find(radius_sorted);
length_sorted = z00(index_sort);

% figure(8);
% test_Plot = plot(radius_sorted,length_sorted,'o')
% title({'Sorted length versus node number';''},'FontWeight','normal')
% grid on; grid minor; xlabel('Node number'); ylabel('z (m)');
% ax = gca; ax.FontSize = 14;

figure(9); plot(radius_sorted,'o');
title({'Sorted radius versus node number';''},'FontWeight','normal')
grid on; grid minor; xlabel('Node number'); ylabel('r (m)');
ax = gca; ax.FontSize = 14;

% figure(10); plot(length_sorted,'o');
% title({'Sorted length versus node number';''},'FontWeight','normal')
% grid on; grid minor; xlabel('Node number'); ylabel('r (m)');
% ax = gca; ax.FontSize = 14;

nodesTop_l_h = find(and(length_P - eps <= z00,z00 <= length_P + eps)); % identify middle nodes
xy1=point(:,nodesTop_l_h); x1=xy1(1,:); y1=xy1(2,:); % top nodes (at z = length)

nodesCenter_h = find(and(0.5*length_P - eps <= z00,z00 <= 0.5*length_P + eps)); % identify middle nodes
xy2=point(:,nodesCenter_h); x2=xy2(1,:); y2=xy2(2,:); % central nodes (middle of the pipe) alogn the length (Z = length/2)

nodesExterior_v_O = find(and(offset - eps <= r00, r00 <= offset + eps)); % identify exterior nodes of the pipe *****
nodesExterior_v_P = find(and(radius_P - eps <= r00, r00 <= radius_P + eps)); % identify exterior nodes of the pipe *****
nodesExterior_v_E = find(and(radius_E - eps <= r00, r00 <= radius_E + eps)); % identify exterior nodes of the heat sink *****
xy3_O=point(:,nodesExterior_v_O); x3_O=xy3_O(1,:); y3_O=xy3_O(2,:); % interior nodes (interior surface of the pipe) (r = Offset)
xy3_P=point(:,nodesExterior_v_P); x3_P=xy3_P(1,:); y3_P=xy3_P(2,:); % exterior nodes (exterior surface of the pipe) (r = Radius_P)
xy3_E=point(:,nodesExterior_v_E); x3_E=xy3_E(1,:); y3_E=xy3_E(2,:); % exterior nodes (exterior surface of the heat sink) (r = Radius_E)

nodesInterior_1 = find(and(-eps + offset <= x2, x2 <= + eps + offset)); % nodes at the fluid interface

nodesInterior_2_O = find(and(-eps + radius_O <= x2, x2 <= + eps + radius_O)); % nodes at the fluid interface
nodesInterior_2_P = find(and(-eps + radius_P <= x2, x2 <= + eps + radius_P)); % nodes at the fluid interface
nodesInterior_2_E = find(and(-eps + radius_E <= x2, x2 <= + eps + radius_E)); % nodes at the fluid interface
xy4_O=xy2(:,nodesInterior_2_O); x4_O=xy4_O(1,:); y4_O=xy4_O(2,:); % exterior nodes (exterior surface of the pipe) (r = Radius_P, z = length/2)
xy4_P=xy2(:,nodesInterior_2_P); x4_P=xy4_P(1,:); y4_P=xy4_P(2,:); % exterior nodes (exterior surface of the pipe) (r = Radius_P, z = length/2)
xy4_E=xy2(:,nodesInterior_2_E); x4_E=xy4_E(1,:); y4_E=xy4_E(2,:); % exterior nodes (exterior surface of the heat sink) (r = Radius_E, z = length/2)

Tintrp1 = interpolateTemperature(resultT,x1,y1,1:time_l);
Tintrp2 = interpolateTemperature(resultT,x2,y2,1:time_l);
Tintrp3_O = interpolateTemperature(resultT,x3_O,y3_O,1:time_l);
Tintrp3_P = interpolateTemperature(resultT,x3_P,y3_P,1:time_l);
Tintrp3_E = interpolateTemperature(resultT,x3_E,y3_E,1:time_l);
Tintrp4_O = interpolateTemperature(resultT,x4_O,y4_O,1:time_l);
Tintrp4_P = interpolateTemperature(resultT,x4_P,y4_P,1:time_l);
Tintrp4_E = interpolateTemperature(resultT,x4_E,y4_E,1:time_l);

[Tx1,Ty1] = evaluateTemperatureGradient(resultT,x1,y1,1:time_l);
[Tx2,Ty2] = evaluateTemperatureGradient(resultT,x2,y2,1:time_l);
[Tx3_O,Ty3_O] = evaluateTemperatureGradient(resultT,x3_O,y3_O,1:time_l);
[Tx3_P,Ty3_P] = evaluateTemperatureGradient(resultT,x3_P,y3_P,1:time_l);
[Tx3_E,Ty3_E] = evaluateTemperatureGradient(resultT,x3_E,y3_E,1:time_l);
[Tx4_O,Ty4_O] = evaluateTemperatureGradient(resultT,x4_O,y4_O,1:time_l);
[Tx4_P,Ty4_P] = evaluateTemperatureGradient(resultT,x4_P,y4_P,1:time_l);
[Tx4_E,Ty4_E] = evaluateTemperatureGradient(resultT,x4_E,y4_E,1:time_l);

[qx,qy] = evaluateHeatFlux(resultT);
[qx1,qy1] = evaluateHeatFlux(resultT,x1,y1,1:time_l);
[qx2,qy2] = evaluateHeatFlux(resultT,x2,y2,1:time_l);
[qx3_O,qy3_O] = evaluateHeatFlux(resultT,x3_O,y3_O,1:time_l);
[qx3_P,qy3_P] = evaluateHeatFlux(resultT,x3_P,y3_P,1:time_l);
[qx3_E,qy3_E] = evaluateHeatFlux(resultT,x3_E,y3_E,1:time_l);
[qx4_O,qy4_O] = evaluateHeatFlux(resultT,x4_O,y4_O,1:time_l);
[qx4_P,qy4_P] = evaluateHeatFlux(resultT,x4_P,y4_P,1:time_l);
[qx4_E,qy4_E] = evaluateHeatFlux(resultT,x4_E,y4_E,1:time_l);

% Plot nodes in the middle of the pipe (z = length/2) at the offset, pipe
% x1: top nodes (at z = length)
% x2: central nodes (middle of the pipe) alogn the length (z = length/2)
% x3_O: interior nodes (interior surface of the pipe) (r = Offset)
% x3_P: exterior nodes (exterior surface of the pipe) (r = Radius_P)
% x3_E: exterior nodes (exterior surface of the heat sink) (r = Radius_E)
% x4_O; interior nodes (interior surface of the pipe) (r = Offset, z = length/2)
% x4_P; exterior nodes (exterior surface of the pipe) (r = Radius_P, z = length/2)
% x4_E: exterior nodes (exterior surface of the heat sink) (r = Radius_E, z = length/2)

% internal boundary, insulation external boundary
figure(10);
hold all;
plot(x1,y1,'s','MarkerSize',3,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',2:1:length(x1));
plot(x2,y2,'+','MarkerSize',3,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',2:1:length(x2));
title({'Radial nodes at pipe''s top surface and midplane';''},'FontWeight','normal');
plot(x3_O,y3_O,'s','MarkerSize',3,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:length(x3_O));
plot(x3_P,y3_P,'d','MarkerSize',3,'MarkerFaceColor',[1 1 1],'MarkerIndices',1:2:length(x3_P));
plot(x3_E,y3_E,'+','MarkerSize',3,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:1:length(x3_E));
legend('Top','Midplane','Interior','Interface','Exterior');
%title({'Axial nodes at pipe''s interior, interface,';'and exterior surfaces';''},'FontWeight','normal');
%plot(x4_O,y4_O,'s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:1:length(x4_O));
%plot(x4_P,y4_P,'d','MarkerSize',5,'MarkerFaceColor',[1 1 1],'MarkerIndices',1:1:length(x4_P));
%plot(x4_E,y4_E,'+','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:1:length(x4_E));
%legend('Interior','Interface','Exterior');
%title({'Nodes at pipe''s interior, interface,';'and exterior surfaces';''},'FontWeight','normal');
%title({'Radial nodes at pipe''s middle';''},'FontWeight','normal');
grid on; grid minor; xlabel ('r (m)'); ylabel ('z (m)'); 
ylim([0 length_P]); xlim([offset-2*grid_size radius_E+2*grid_size]);
ax = gca; ax.FontSize = 14;

% Plot axial-central nodes tempeature
figure(11); 
hold all;
p1 = plot(x1,Tintrp1(:,time_l),'s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:length(x1));
p2 = plot(x2,Tintrp2(:,time_l),'+','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:length(x2));
title({'Temperature at pipe''s top surface and midplane'; ['(t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('r (m)'); ylabel('T (degC)');
ax = gca; ax.FontSize = 14; legend('Top','Midplane'); 

figure(12);
hold all;
p3 = plot(x1,Tx1(:,time_l),'s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:length(x1));
p4 = plot(x2,Tx2(:,time_l),'+','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:length(x2));
title({'Temperature gradient at pipe''s top surface and midplane'; ['(t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('r (m)'); ylabel('dT/dr (degC/m)');
ax = gca; ax.FontSize = 14; legend('Top','Midplane','location','east');

figure(13);
hold all;
copper_x1_P = find(x1 <= radius_P + 0.5*eps); % identify middle nodes
aluminum_x1_E = find(and(radius_P - 0.5*eps <= x1,x1 <= radius_E + 0.5*eps)); % identify middle nodes
[Tx1_C,Ty1_C] = evaluateTemperatureGradient(resultT,sort(x1(copper_x1_P)),y1(copper_x1_P),1:time_l);
[Tx1_Al,Ty1_Al] = evaluateTemperatureGradient(resultT,sort(x1(aluminum_x1_E)),y1(aluminum_x1_E),1:time_l);
p3 = plot(sort(x1(copper_x1_P)),-k1*Tx1_C(:,time_l),'-s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p4 = plot(sort(x1(aluminum_x1_E)),-k2*Tx1_Al(:,time_l),'-+','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Heat flux (-k*dT/dr) at pipe''s top surface and midplane'; ['(t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('r (m)'); ylabel('-k*dT/dr (W/m^2)');
ax = gca; ax.FontSize = 14; legend('Top','Midplane');

% Plot axial-central nodes heat flux
figure(14); 
hold all;
p5 = plot(x1,qx1(:,time_l),'s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:length(x1));
p6 = plot(x2,qx2(:,time_l),'+','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:length(x2));
title({'Transient heat rate at pipe''s top surface and midplane'; ['(t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('r (m)'); ylabel('q (W)');
ax = gca; ax.FontSize = 14; legend('Top','Midplane');

% Plot axial-central heat rate at pipe''s exterior surface
figure(15); 
hold all;
mm = 50;
mn = (1:mm);
list = tlist(mn);
Qn_fin_horizontal_1 = evaluateHeatRate(resultT,'Edge',12); % heat rate at pipe''s exterior surface of pipe in contact with the heat sink
Qn_fin_vertical_1 = evaluateHeatRate(resultT,'Edge',52); % heat rate at pipe''s exterior surface of pipe in contact with the heat sink
Qn_fin_horizontal_2 = evaluateHeatRate(resultT,'Edge',11); % heat rate at pipe''s exterior surface of pipe in contact with the heat sink
Qn_fin_vertical_2 = evaluateHeatRate(resultT,'Edge',34); % heat rate at pipe''s exterior surface of pipe in contact with the heat sink
Qn_fin_horizontal_3 = evaluateHeatRate(resultT,'Edge',35); % heat rate at pipe''s exterior surface of pipe in contact with the heat sink
p5_horizontal_1 = plot(list,Qn_fin_horizontal_1(mn),'s','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:1:mm);
p5_vertical_1 = plot(list,Qn_fin_vertical_1(mn),'d','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:1:mm);
p5_horizontal_2 = plot(list,Qn_fin_horizontal_2(mn),'+','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:1:mm);
p5_Contact = plot(list,Qn_fin_vertical_2(mn),'>','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:1:mm);
p5_horizontal_3 = plot(list,Qn_fin_horizontal_3(mn),'<','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:1:mm);
title({'Trasient heat rate per edge length at pipe''s edges';''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('t (s)'); ylabel('q (W/m)');
%xlim([0 tlist]); 
ax = gca; ax.FontSize = 14; legend('Edge 12','Edge 52','Edge 11','Edge 34','Edge 35');

heat_flux_P = qx2(time_l) % W/m^2
heat_flux_per_length = heat_flux_P * 2*pi()*radius_E % W/m
% 1 W/ m = 1.04 BTU/h.ft

figure(16);
pdeplot(thermalModelT,'XYData',T1(:,time_l),'Contour','on','FlowData',[qx,qy],'ColorMap','hot')
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
title({'Temperature contours (degC) and heat flux vector field'; ['(t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
ax = gca; ax.FontSize = 14;

% Plot axial-central nodes tempeature
figure(17);
hold all;
p6_O = plot(y3_O,Tintrp3_O(:,time_l),'s','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:mm);
p6_P = plot(y3_P,Tintrp3_P(:,time_l),'d','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:mm);
p6_E = plot(y3_E,Tintrp3_E(:,time_l),'+','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:mm);
title({'Temperature at pipe''s interior, interface,';['and exterior surfaces (t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
xlim([0 length_P])
grid on; grid minor; xlabel('z (m)'); ylabel('T (degC)');
ax = gca; ax.FontSize = 14; legend('Interior','Interface','Exterior','location','best');

% Plot axial-central nodes tempeature
figure(18);
hold all;
p6_O = plot(y3_O,Tx3_O(:,time_l),'s','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:mm);
p6_P = plot(y3_P,Tx3_P(:,time_l),'d','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:mm);
p6_E = plot(y3_E,Tx3_E(:,time_l),'+','MarkerSize',5,...
    'MarkerFaceColor',[0.5 0.5 0.5],'MarkerIndices',1:2:mm);
title({'Temperature gradient at pipe''s interior, interface,';['and exterior surfaces (t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('z (m)'); ylabel('dT/dr (degC/m)');
ax = gca; ax.FontSize = 14; legend('Interior','Interface','Exterior');

figure(19)
hold all;
p32_O = plot(y3_O(:),-k1*Tx3_O(:,time_l),'s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);  
p32_P = plot(y3_P(:),-k1*Tx3_P(:,time_l),'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p32_E = plot(y3_E(:),-k2*Tx3_E(:,time_l),'+','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Heat flux at pipe''s interior, interface,';['and exterior surfaces (t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('z (m)'); ylabel('-k*dT/dr (W/m^2)');
ax = gca; ax.FontSize = 14; legend('Interior','Interface','Exterior');

figure(20)
hold all;
p32_O = plot(x4_O(:),qx4_O(:,time_l),'s','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]);  
p32_P = plot(x4_P(:),qx4_P(:,time_l),'d','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]);
p32_E = plot(x4_E(:),qx4_E(:,time_l),'>','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Heat rate at pipe''s interior, interface,';['and exterior surfaces (t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('r (m)'); ylabel('q (W)');
ax = gca; ax.FontSize = 14; legend('Interior','Interface','Exterior','location','best');

figure(21)
hold all;
p32_O = plot(x4_O(:),Tx4_O(:,time_l),'s','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]);  
p32_P = plot(x4_P(:),Tx4_P(:,time_l),'d','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]);
p32_E = plot(x4_E(:),Tx4_E(:,time_l),'>','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Temperature gradient at pipe''s interior, interface,';['and exterior surfaces (t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('r (m)'); ylabel('dT/dr (degC/m)');
ax = gca; ax.FontSize = 14; legend('Interior','Interface','Exterior','location','best');

figure(22)
hold all;
p32_O = plot(x4_O(:),Tintrp4_O(:,time_l),'s','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]);  
p32_P = plot(x4_P(:),Tintrp4_P(:,time_l),'d','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]);
p32_E = plot(x4_E(:),Tintrp4_E(:,time_l),'>','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Temperature at pipe''s interior, interface,';['and exterior surfaces (t = ',...
    num2str(ttt),' min)'];''},'FontWeight','normal') %time_l
grid on; grid minor; xlabel('r (m)'); ylabel('T (degC)');
ax = gca; ax.FontSize = 14; legend('Interior','Interface','Exterior','location','best');

figure(23); 
hold all;
p41 = plot(tlist,Tintrp4_O(:,:),'-','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5]);
p42 = plot(tlist,Tintrp4_P(:,:),'-d','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5]);
p43 = plot(tlist,Tintrp4_E(:,:),'-+','MarkerSize',4,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Transient temperature at pipe''s interior, interface and';'exterior surfaces';''},'FontWeight','normal');
grid on; grid minor; xlabel('t (s)'); ylabel('T (degC)');
xlim([0 15]); ax = gca; ax.FontSize = 14; 
legend('Interior','Interface','Exterior','location','East');

figure(24); 
hold all;
p41 = plot(tlist,Tx4_O(:,:),'-s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p42 = plot(tlist,Tx4_P(:,:),'-d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p43 = plot(tlist,Tx4_E(:,:),'-+','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Transient temperature gradient at pipe''s interface and';'exterior surfaces';''},'FontWeight','normal');
grid on; grid minor; xlabel('t (s)'); ylabel('dT/dr (degC/m)');
xlim([0 15]); ax = gca; ax.FontSize = 14; 
legend('Interface','Exterior','location','East');

