% Author: Layla S. Mayboudi
% 2021-05-27
% This code creates a stationary thermal PDE model for an axisymmetric pipe. 
% Geometry, temperature distribution, and temperature gradients diagrams
% are plotted.
% A cost comparison is eventually made between the energy waste from the
% copper, PEX, and PEX-covered copper pipes.

% Note: 
% The codes are to be run at the following order, for their inputs are
% required for the consequent steps:
% 1- Run the cu pipe (thermal_model_pipe_steady_Model_Cu.m).
% 2- Run the PEX pipe (thermal_model_pipe_steady_Model_PEX.m).
% 3- Run the PEX-covered cu pipe (thermal_model_pipe_steady_Model_Cu_PEX.m).

close all; clf, clc; %clear all;
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')
  
  
thermalModelS = createpde('thermal','steadystate-axisymmetric'); % create a model object: thermal

insulation_th = 0.5*25.4/1000; %m
length_P = 10/1000; %m

radius_P = 1.125*25.4/2000; offset = 0.995*25.4/2000;  radius_E = insulation_th+radius_P;     %pipe diameter 1" %m
geom_P = [3;4;offset;offset;radius_P;radius_P; 0;length_P;length_P;0]; % create geometry: rectangle solid
geom_I = [3;4;radius_P;radius_P;radius_E;radius_E; 0;length_P;length_P;0]; % create geometry: rectangle solid

gd = [geom_P,geom_I];
ns = char('geom_P','geom_I');
ns = ns';
sf = 'geom_P+geom_I';
dl = decsg(gd,sf,ns);
geometryFromEdges(thermalModelS,dl); % associate geometry and thermal model object
figure(1);
pdegplot(thermalModelS,'EdgeLabels','on','FaceLabels','on');
ylim([0-1/1000 length_P+1/1000]); xlim([0 radius_E+offset+1/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 = 392; % 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(thermalModelS,'ThermalConductivity',k1,'Face',1) 

% Insulating material properties - transient thermal analysis (fiberglass, FG)
k2 = 0.04; % copper thermal conductivity (W/mK) 
k_FG = k2;
cp2 = 700; % copper heatcapacity (J/kgK) 
rho2 = 150 % copper density (kg/m^3) 
q1 = 0; %heat generation (W/m^3)
tp2=thermalProperties(thermalModelS,'ThermalConductivity',k2,'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
thermalBC(thermalModelS,'Edge',2,'ConvectionCoefficient',hc,'AmbientTemperature',T_ambient); % convection heat transfer
thermalBC(thermalModelS,'Edge',1,'Temperature',T_water); % constant temperature

% Plot geometry
figure(2);
grid_size = (radius_E-offset)/20; % identify grid size
mesh1 = generateMesh(thermalModelS,'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
thermalModelS.SolverOptions.ReportStatistics = 'on';
results = solve(thermalModelS); % identify PDESolverOptions Properties
T1_FG = results.Temperature; % identify data to process and plot
Tx1_FG = results.XGradients; % identify data gradient along r-coordinate
Ty1_FG = results.YGradients; % identify data gradient along z-coordinate

% Contour plot thermal model steady temperature versus rz
figure(5);
pdeplot(thermalModelS,'XYData',T1_FG(:,end),'colormap','hot','contour','on','levels',12); % 'colormap','hot' ,'Contour','on'
ylim([0 length_P]); xlim([offset radius_E+2*grid_size]);
title({'Temperature Contours (degC)';''},'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(thermalModelS,'XYData',Tx1_FG(:,end),'colorMap','hot','contour','on','levels',10); 
ylim([0 length_P]); xlim([offset radius_E+2*grid_size]);
title({'Temperature Gradient (dT/dr, degC/m)';''},'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(thermalModelS,'XYData',Ty1_FG(:,end),'colormap','hot','contour','on','levels',10); 
ylim([0 length_P]); xlim([offset radius_E+2*grid_size]);
title({'Temperature Gradient (dT/dz, degC/m)';''},'FontWeight','normal')
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14;
 
%**********************************************************
%Investigate the node density

point_FG = thermalModelS.Mesh.Nodes; % recall node locations from thermal model object
eps1=eps;
nodesIndex_FG = find(point_FG);
xy00_FG=point_FG; r00_FG = xy00_FG(1,:); z00_ZG = xy00_FG(2,:);
length_sorted = sort(z00_ZG); radius_sorted = sort(r00_FG); 

figure(8); plot(length_sorted); 
title({'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);
title({'Radius versus Node Number';''},'FontWeight','normal')
grid on; grid minor; xlabel('Node number'); ylabel('r (m)');
ax = gca; ax.FontSize = 14;

nodesTopEnd_h_FG = find(z00_ZG == length_P); % identify top nodes (at z = length)
nodesCenter_h_FG = find(and(0.5*length_P - 0.5*grid_size <= z00_ZG,z00_ZG <= 0.5*length_P + 0.5*grid_size)); % identify middle nodes
nodesExterior_v_O_FG = find(and(offset - eps <= r00_FG, r00_FG <= offset + eps)); % identify exterior nodes *****
nodesExterior_v_P_FG = find(and(radius_P - eps <= r00_FG, r00_FG <= radius_P + eps)); % identify exterior nodes *****
nodesExterior_v_E_FG = find(and(radius_E - eps <= r00_FG, r00_FG <= radius_E + eps)); % identify exterior nodes *****
xy1_FG=point_FG(:,nodesTopEnd_h_FG); x1_FG=xy1_FG(1,:); y1_FG=xy1_FG(2,:);
xy2_FG=point_FG(:,nodesCenter_h_FG); x2_FG=xy2_FG(1,:); y2_FG=xy2_FG(2,:); % central nodes (middle of the pipe) alogn the length (Z = length/2)
xy3_O_FG=point_FG(:,nodesExterior_v_O_FG); x3_O_FG=xy3_O_FG(1,:); y3_O_FG=xy3_O_FG(2,:); % exterior nodes (exterior surface of the pipe) (r = Radius_P)
xy3_P_FG=point_FG(:,nodesExterior_v_P_FG); x3_P_FG=xy3_P_FG(1,:); y3_P_FG=xy3_P_FG(2,:); % exterior nodes (exterior surface of the pipe) (r = Radius_P)
xy3_E_FG=point_FG(:,nodesExterior_v_E_FG); x3_E_FG=xy3_E_FG(1,:); y3_E_FG=xy3_E_FG(2,:); % exterior nodes (exterior surface of the insulation) (r = Radius_E)
nodesInterior_1_FG = find(and(-eps + offset <= x2_FG, x2_FG <= + eps + offset)); % identify nodes at the fluid interface
nodesInterior_2_P_FG = find(and(radius_P - eps <= x2_FG, x2_FG <= radius_P + eps)); % identify exterior nodes for the pipe *****
nodesInterior_2_E_FG = find(and(radius_E - eps <= x2_FG, x2_FG <= radius_E + eps)); % identify exterior nodes for the insulation *****
xy5_FG=xy2_FG(:,nodesInterior_1_FG); x5_FG=xy5_FG(1,:); y5_FG=xy5_FG(2,:); % identify nodes at the fluid interface
xy6_P=xy2_FG(:,nodesInterior_2_P_FG); x6_P=xy6_P(1,:); y6_P=xy6_P(2,:); % exterior nodes (exterior surface of the pipe) (r = Radius_P, z = length/2)
xy6_E=xy2_FG(:,nodesInterior_2_E_FG); x6_E=xy6_E(1,:); y6_E=xy6_E(2,:); % exterior nodes (exterior surface of the insulation) (r = Radius_E, z = length/2)

% Plot nodes in the middle of the pipe (z = length/2) at the offset, pipe
% internal boundary, insulation external boundary
figure(10);
hold all;
% plot(x3_O,y3_O,'-','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]); % radial nides at pipe''s middle
% plot(x3_P,y3_P,'d','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]); % radial nides at pipe''s middle
% plot(x3_E,y3_E,'+','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]); % radial nides at pipe''s middle
plot(x5_FG,y5_FG,'s','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]); % radial nides at pipe''s middle
plot(x6_P,y6_P,'d','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]); % radial nides at pipe''s middle
plot(x6_E,y6_E,'+','MarkerSize',7,'MarkerFaceColor',[0.5 0.5 0.5]); % radial nides at pipe''s middle
ylim([0 length_P]); xlim([offset-2*grid_size radius_E+2*grid_size]);
title({'Radial Nodes at Pipe''s Middle';''},'FontWeight','normal');
grid on; grid minor; xlabel ('r (m)'); ylabel ('z (m)'); 
ax = gca; ax.FontSize = 14; legend;

% Plot radial-central nodes
figure(11);
hold all;
plot(x2_FG,y2_FG,'o','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]); % radial nides at pipe''s middle
ylim([0 length_P]); xlim([offset-2*grid_size radius_E+2*grid_size]);
title({'Radial Nodes at Pipe''s Middle';''},'FontWeight','normal');
grid on; grid minor; xlabel ('r (m)'); ylabel ('z (m)'); 
ax = gca; ax.FontSize = 14;

% Plot axial-exterior nodes
figure(12);
hold all;
plot(x3_O_FG,y3_O_FG,'s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
plot(x3_P_FG,y3_P_FG,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
plot(x3_E_FG,y3_E_FG,'+','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
ylim([0 length_P]); xlim([offset-2*grid_size radius_E+2*grid_size]);
title({'Axial Nodes at Pipe''s Exterior';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14; legend;

% Plot axial-central nodes tempeature
figure(13); 
hold all;
Tintrp2_FG = interpolateTemperature(results,x2_FG,y2_FG);
p2 = plot(x2_FG,Tintrp2_FG,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Temperature at Pipe''s Middle';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('T (degC)');
ax = gca; ax.FontSize = 14;

figure(14)
hold all;
[Tx2_FG,Ty2_FG] = evaluateTemperatureGradient(results,x2_FG,y2_FG);
p3 = plot(x2_FG,Tx2_FG,'x','MarkerSize',4);    
title({'Temperature Gradient';''},...
    'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('dT/dr (degC/m)');
ax = gca; ax.FontSize = 14;

% Plot axial-central nodes heat flux
figure(15); 
hold all;
[qx2_FG,qy2_FG] = evaluateHeatFlux(results,x2_FG,y2_FG); % heat flux at pipe''s exterior surface
p4 = plot(x2_FG,qx2_FG,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Radial Heat Flux at Pipe''s Middle';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('qr (W/m^2)');
ax = gca; ax.FontSize = 14;

% Plot axial-central heat rate at pipe''s exterior surface
figure(16); 
hold all;
Qn_exterior_P_FG = evaluateHeatRate(results,'Edge',3); % heat rate at pipe''s exterior surface
Qn_exterior_E_FG = evaluateHeatRate(results,'Edge',2); % heat rate at pipe''s exterior surface
p5_P = plot(x3_P_FG,Qn_exterior_P_FG,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p5_E = plot(x3_E_FG,Qn_exterior_E_FG,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Heat Rate at Pipe''s Interior and Exterior Edges';''},'FontWeight','normal');
xlim([offset-2*grid_size radius_E+2*grid_size]);
grid on; grid minor; xlabel('r (m)'); ylabel('qr (W/m^2)');
ax = gca; ax.FontSize = 14;

% heat_flux_P = qx2(end) % 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(17);
[qx,qy] = evaluateHeatFlux(results);
pdeplot(thermalModelS,'XYData',T1_FG(:,end),'Contour','on','FlowData',[qx,qy],'ColorMap','hot')
grid on; grid minor; xlabel('r (m)'); ylabel('qr (W/m^2)');
title({'Temperature Contours and Heat Flux Vector Field';''},'FontWeight','normal');
ax = gca; ax.FontSize = 14;

% Plot axial-central nodes tempeature
figure(18); 
hold all;
Tintrp3_O_FG = interpolateTemperature(results,x3_O_FG,y3_O_FG);
Tintrp3_P_FG = interpolateTemperature(results,x3_P_FG,y3_P_FG);
Tintrp3_E_FG = interpolateTemperature(results,x3_E_FG,y3_E_FG);
p6_P = plot(x3_P_FG,Tintrp3_P_FG,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p6_E = plot(x3_E_FG,Tintrp3_E_FG,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Radial Temperature at Pipe''s Interior and Exterior Edges';''},'FontWeight','normal');
xlim([offset-2*grid_size radius_E+2*grid_size]);
grid on; grid minor; xlabel('r (m)'); ylabel('T (degC)');
ax = gca; ax.FontSize = 14;

figure(19)
hold all;
[Tx3_O_FG,Ty3_O_FG] = evaluateTemperatureGradient(results,x3_O_FG,y3_O_FG);
[Tx3_P_FG,Ty3_P_FG] = evaluateTemperatureGradient(results,x3_P_FG,y3_P_FG);
[Tx3_E_FG,Ty3_E_FG] = evaluateTemperatureGradient(results,x3_E_FG,y3_E_FG);
p32_P = plot(x3_O_FG,-k1*Tx3_O_FG,'>','MarkerSize',4);    
p32_E = plot(x3_P_FG,-k2*Tx3_P_FG,'x','MarkerSize',4);    
p32_E = plot(x3_E_FG,-k2*Tx3_E_FG,'o','MarkerSize',4);    
title({'Radial Heat Flux at Pipe''s Interior and Exterior Edges';''},...
    'FontWeight','normal');
xlim([offset-2*grid_size radius_E+2*grid_size]);
grid on; grid minor; xlabel('r (m)'); ylabel('-k*dT/dr (W/m^2)');
ax = gca; ax.FontSize = 14; legend;

% Plot radial heat rate
figure(20); 
hold all;
[qx3_O,qy3_O] = evaluateHeatFlux(results,x3_O_FG,y3_O_FG); % heat flux at pipe''s exterior surface
[qx3_P_FG,qy3_P_FG] = evaluateHeatFlux(results,x3_P_FG,y3_P_FG); % heat flux at pipe''s exterior surface
[qx3_E_FG,qy3_E_FG] = evaluateHeatFlux(results,x3_E_FG,y3_E_FG); % heat flux at pipe''s exterior surface
p42_O = plot(x3_O_FG,qx3_O,'>','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p42_P = plot(x3_P_FG,qx3_P_FG,'x','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p42_E = plot(x3_E_FG,qx3_E_FG,'o','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Radial Heat Rate at Pipe''s Interior and Exterior Edges';''},...
    'FontWeight','normal');
xlim([offset-2*grid_size radius_E+2*grid_size]);
grid on; grid minor; xlabel('r (m)'); ylabel('qr (W/m)');
ax = gca; ax.FontSize = 14; legend;

%*************************
% Plot central heat flux vectors (W/m^2)
figure(21)
hold all; 
[qx2_FG,qy2_FG] = evaluateHeatFlux(results,x2_FG,y2_FG);
quiver(x2_FG(:),y2_FG(:),qx2_FG(:),qy2_FG(:),'AutoScale','on','color','blue'); %'MarkerSize',2,'LineWidth',2
ylim([0 length_P]); xlim([offset-2*grid_size radius_E+2*grid_size]);
title({'Radial Heat Flux at Pipe''s Middle (W/m^2)';''},'FontWeight','normal');
grid on; grid minor; xlabel 'r (m)'; ylabel 'z (m)';
ax = gca; ax.FontSize = 14;
 
% Plot exterior heat flux vectors (W/m^2)
figure(22)
hold all; 
[qx3_O_FG,qy3_O_FG] = evaluateHeatFlux(results,x3_P_FG,y3_P_FG);
[qx3_P_FG,qy3_P_FG] = evaluateHeatFlux(results,x3_P_FG,y3_P_FG);
[qx3_E_FG,qy3_E_FG] = evaluateHeatFlux(results,x3_E_FG,y3_E_FG);
quiver(x3_O_FG(:),y3_O_FG(:),qx3_O_FG(:),qy3_O_FG(:),'AutoScale','on','color','blue'); %'MarkerSize',2,'LineWidth',2
quiver(x3_P_FG(:),y3_P_FG(:),qx3_P_FG(:),qy3_P_FG(:),'AutoScale','on','color','blue'); %'MarkerSize',2,'LineWidth',2
quiver(x3_E_FG(:),y3_E_FG(:),qx3_E_FG(:),qy3_E_FG(:),'AutoScale','on','color','blue'); %'MarkerSize',2,'LineWidth',2
ylim([0 length_P]); xlim([offset-2*grid_size radius_E+2*grid_size]);
title({'Axial Heat Flux at Pipe''s Exterior Surface (W/m^2)';''},'FontWeight','normal');
grid on; grid minor; xlabel 'r (m)'; ylabel 'z (m)';
ax = gca; ax.FontSize = 14; legend('Interior edge','Exterior edge')

pipe_area_P_FG = (2*pi*radius_P)*length_P % m^2
pipe_area_E_FG = (2*pi*radius_E)*length_P % m^2

flux_P_FG = -k1*Tx3_P_FG(end) % W
flux_E_FG = -k2*Tx3_E_FG(end) % W
flux_W_per_length_P_FG = flux_P_FG*pipe_area_P_FG/length_P % W
flux_W_per_length_E_FG = flux_E_FG*pipe_area_E_FG/length_P % W
flux_BTU_per_length_P_FG = flux_W_per_length_P_FG / (0.293071*3.28084) % BTU/hr/m
flux_BTU_per_length_E_FG = flux_W_per_length_E_FG / (0.293071*3.28084) % BTU/hr/m

%%% Plot axial-central nodes tempeature - Cu
figure(23); 
hold all;
Tintrp2_FG = interpolateTemperature(results,x2_FG,y2_FG);
p2 = plot(x2_FG,Tintrp2_FG,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Radial Temperature at Pipe''s Middle';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('T (degC)');
ax = gca; ax.FontSize = 14;

figure(24)
hold all;
[Tx3_O_FG,Ty3_O_FG] = evaluateTemperatureGradient(results,x3_P_FG,y3_P_FG);
[Tx3_P_FG,Ty3_P_FG] = evaluateTemperatureGradient(results,x3_P_FG,y3_P_FG);
[Tx3_E_FG,Ty3_E_FG] = evaluateTemperatureGradient(results,x3_E_FG,y3_E_FG);
p3_O= plot(x3_O_FG,-k1*Tx3_O_FG,'x','MarkerSize',4);    
p3_P = plot(x3_P_FG,-k1*Tx3_P_FG,'x','MarkerSize',4);    
p3_E = plot(x3_E_FG,-k2*Tx3_E_FG,'x','MarkerSize',4);    
title({'Radial Heat Flux (dT/dr)';''},...
    'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('dT/dr (degC/m)');
ax = gca; ax.FontSize = 14; legend;

%********************
% Compare the PEX, Cu, and Insulated pipes

% figure(25)
% tt1 = plot(x2_PEX,Tintrp2_PEX,'black *',x2_Cu,Tintrp2_Cu,'black o',x2_FG,Tintrp2_FG,'black +','MarkerSize',5)
% title({'Radial Temperature at Pipe''s Middle';''},'FontWeight','normal');
% legend('PEX','Cu: no insulation','Cu: insulation');
% grid on; grid minor; xlabel('r (m)'); ylabel('T (degC)');
% ax = gca; ax.FontSize = 14;
% 
% figure(26)
% tt2 = plot(x2_PEX,-k_PEX*Tx2_PEX,'black *',x2_Cu,-k_Cu*Tx2_Cu,'black o',x3_P_FG,-k_Cu*Tx3_P,'black p',x3_E_FG,-k_FG*Tx3_E,'black s');
% title({'Radial Heat Flux at Pipe''s Middle (dT/dr)';''},'FontWeight','normal');
% legend('PEX','Cu: no insulation','Cu: insulation','Insulation');
% grid on; grid minor; xlabel('r (m)'); ylabel('-k*dT/dr (W/m^2)');
% ax = gca; ax.FontSize = 14;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plot -central nodes tempeature - PEX and Cu
%  Note: To compare PEX and Cu results, the Cu file is to be run first; otherwise,
% comment the Cu sections below.

figure(27); 
hold all;
p21_PEX = plot(x2_PEX,Tintrp2_PEX,'>','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p22_Cu = plot(x2_Cu,Tintrp2_Cu,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p23_FG = plot(x2_FG,Tintrp2_FG,'s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Temperature at Pipe''s Midplane';''},'FontWeight','normal');
%ylim([75 81])
grid on; grid minor; xlabel('r (m)'); ylabel('T (degC)');
ax = gca; ax.FontSize = 14; legend('PEX','Cu','FG');

figure(28)
hold all;
p24_PEX = plot(x2_PEX,-k_PEX*Tx2_PEX,'>','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p25_Cu = plot(x2_Cu,-k_Cu*Tx2_Cu,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p26_FG = plot(x2_FG,-k_FG*Tx2_FG,'s','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Heat Flux';''},...
    'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('-k*dT/dr (W/m^2)');
ax = gca; ax.FontSize = 14; legend('PEX','Cu','FG');

pipe_area_PEX = (2*pi*radius_P)*length_P % m^2
flux_P_PEX = -k_PEX*Tx3_PEX(end) % flux at the pipe's exterior surface (W/m^2 )
flux_W_per_length_PEX = flux_P_PEX*pipe_area_PEX/length_P % heat flux at the pipe's exterior surface per length (W/m)

pipe_area_Cu = (2*pi*radius_P)*length_P % m^2
flux_P_Cu = -k_Cu*Tx3_Cu(end) % heat flux at the pipe's exterior surface (W/m^2 )
flux_W_per_length_Cu = flux_P_Cu*pipe_area_Cu/length_P % heat flux at the pipe's exterior surface per length (W/m)

pipe_area_FG = (2*pi*radius_P)*length_P % m^2
flux_P_FG = -k_FG*Tx3_E_FG(end) % heat flux at the pipe's exterior surface (W/m^2 )
flux_W_per_length_FG = flux_P_FG*pipe_area_FG/length_P % heat flux at the pipe's exterior surface per length (W/m)

% 1 W/ m = 1.04 BTU/h.ft
% flux_BTU_per_length_PEX = flux_W_per_length_PEX / (0.293071*3.28084) % BTU/hr/m
% flux_BTU_per_length_Cu = flux_W_per_length_Cu / (0.293071*3.28084) % BTU/hr/m

Length_t = 25; % m
month = 12;
Time_t = month * 30; %d ays
time_t = Time_t * 24 * 3600; % s

flux_W_PEX = flux_W_per_length_PEX * Length_t % W
flux_W_Cu = flux_W_per_length_Cu * Length_t % W
flux_W_P_FG = flux_W_per_length_P_FG * Length_t % W
flux_W_E_FG = flux_W_per_length_E_FG * Length_t % W

energy_PEX = flux_W_PEX * time_t / 1e6 % J
energy_Cu = flux_W_Cu * time_t / 1e6 % J
%energy_P_FG = flux_W_P_FG * time_t / 1e6 % J
energy_E_FG = flux_W_E_FG * time_t / 1e6 % J

% "Energy." Science & Technology. Gale Research Inc., 1993.
% Standardized results (37 MJ/m^3)

conv = 37; %MJ/m^3

Cubicm_PEX = energy_PEX / conv % m^3
Cubicm_Cu = energy_Cu /conv % m^3
%Cubicm_P_FG = energy_P_FG /conv % m^3
Cubicm_E_FG = energy_E_FG /conv % m^3

% July 1, 2021 natural gas rates Ontario
% EPCOR Natural Gas Limited Partnership (Aylmer)
% The rate for EPCOR includes storage and transportation charges.
% 13.3057 ¢/m3

cost = 13.3057/100; % $/m3

Dollar_PEX = Cubicm_PEX * cost % m^3
Dollar_Cu = Cubicm_Cu * cost % m^3
%Dollar_P = Cubicm_P * cost % m^3
Dollar_E = Cubicm_E_FG * cost % m^3

figure(29)
X = categorical({'PEX','Cu: no insulation','Cu: insulation'});
X = reordercats(X,{'PEX','Cu: no insulation','Cu: insulation'});
Y = [Dollar_PEX Dollar_Cu Dollar_E];
bar(X,Y,'b')
title({'Cost Comparison';''},'FontWeight','normal');
grid on; grid minor; xlabel('Cost ($)'); ylabel('Condition');
