% 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
length_P = 10/1000; %m
radius_P = 1.125*25.4/2000; offset = 0.875*25.4/2000;       %pipe diameter 1" %m
thickness = radius_P - offset;
grid_size = thickness/10 % identify grid size
geom1 = [3;4;offset;offset;radius_P;radius_P; 0;length_P;length_P;0]; % create geometry: rectangle solid

gd = [geom1];
sf = 'geom1';
ns = char('geom1');
ns = ns';
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_P+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 - steady-state thermal analysis
k_PEX = 0.46; % PEX thermal conductivity (W/mK) 
k1 = k_PEX; % to be replaced with any other material
q1 = 0; %heat generation (W/m^3)
tp1=thermalProperties(thermalModelS,'ThermalConductivity',k1) 

% Boundary conditions
T_ambient = 25;  %ambient temperature degC
T_water = 80;   %water temperature degC
hc = 10; % convection heat transfer coefficient W/m^2K
thermalBC(thermalModelS,'Edge',3,'ConvectionCoefficient',hc,'AmbientTemperature',T_ambient); % convection heat transfer
thermalBC(thermalModelS,'Edge',1,'Temperature',T_water); % constant temperature

% Plot geometry
figure(2);

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_P+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_PEX = results.Temperature; % identify data to process and plot
Tx1_PEX = results.XGradients; % identify data gradient along r-coordinate
Ty1_PEX = results.YGradients; % identify data gradient along z-coordinate

% Contour plot thermal model steady temperature versus rz
figure(5);
pdeplot(thermalModelS,'XYData',T1_PEX(:,end),'colormap','hot','contour','on','levels',10); % 'colormap','hot' ,'Contour','on'
ylim([0 length_P]); xlim([offset-2*grid_size radius_P+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_PEX(:,end),'colorMap','hot','contour','on','levels',10); 
ylim([0 length_P]); xlim([offset-2*grid_size radius_P+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_PEX(:,end),'colormap','hot','contour','on','levels',10); 
ylim([0 length_P]); xlim([offset-2*grid_size radius_P+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_PEX = thermalModelS.Mesh.Nodes; % recall node locations from thermal model object
eps1=eps;
nodesIndex_PEX = find(point_PEX);
xy00_PEX=point_PEX; r00_PEX = xy00_PEX(1,:); z00_PEX = xy00_PEX(2,:);
length_sorted = sort(z00_PEX); radius_sorted = sort(r00_PEX); 

nodesTopEnd_h_PEX = find(z00_PEX == length_P); % identify top nodes
nodesCenter_h_PEX = find(and(0.5*length_P - 0.5*grid_size <= z00_PEX,z00_PEX <= 0.5*length_P + 0.5*grid_size)); % identify midplane nodes
nodesExterior_v_PEX = find(and(radius_P - eps <= r00_PEX, r00_PEX <= radius_P + eps)); % identify exterior nodes
xy1_PEX=point_PEX(:,nodesTopEnd_h_PEX); x1_PEX=xy1_PEX(1,:); y1_PEX=xy1_PEX(2,:);
xy2_PEX=point_PEX(:,nodesCenter_h_PEX); x2_PEX=xy2_PEX(1,:); y2_PEX=xy2_PEX(2,:); % central nodes (midplane of the pipe)
xy3_PEX=point_PEX(:,nodesExterior_v_PEX); x3_PEX=xy3_PEX(1,:); y3_PEX=xy3_PEX(2,:); % exterior nodes (exterior surface of the pipe)

nodesInterior_1_PEX = find(and(- eps <= y2_PEX, y2_PEX <= + eps)); % identify top nodes
nodesInterior_2_PEX = find(and(radius_P - eps <= y2_PEX, y2_PEX <= radius_P + eps)); % identify exterior nodes
xy5_PEX=xy2_PEX(:,nodesInterior_1_PEX); x5_PEX=xy5_PEX(1,:); y5=xy5_PEX(2,:);
xy6_PEX=xy2_PEX(:,nodesInterior_2_PEX); x6_PEX=xy6_PEX(1,:); y6=xy6_PEX(2,:);

% Plot -central nodes
figure(8);
hold all;
plot(x2_PEX,y2_PEX,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]); %  nides at pipe''s midplane
ylim([0 length_P]); xlim([offset-2*grid_size radius_P+2*grid_size]);
title({'Nodes at Pipe''s Midplane';''},'FontWeight','normal');
grid on; grid minor; xlabel ('r (m)'); ylabel ('z (m)'); 
ax = gca; ax.FontSize = 14;

% Plot -central nodes tempeature
figure(9); 
hold all;
Tintrp2_PEX = interpolateTemperature(results,x2_PEX,y2_PEX);
p2 = plot(x2_PEX,Tintrp2_PEX,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Temperature at Pipe''s Midplane';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('T (degC)');
ylim([79,81])
ax = gca; ax.FontSize = 14;

figure(10)
hold all;
[Tx2_PEX,Ty2_PEX] = evaluateTemperatureGradient(results,x2_PEX,y2_PEX);
p3 = plot(x2_PEX,Tx2_PEX,'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;

figure(10)
hold all;
[Tx2_PEX,Ty2_PEX] = evaluateTemperatureGradient(results,x2_PEX,y2_PEX);
p3 = plot(x2_PEX,-k1*Tx2_PEX,'x','MarkerSize',4);    
title({'Heat Flux (-k*dT/dr)';''},...
    'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('q (W/m^2)');
ax = gca; ax.FontSize = 14;

%**************************************

% Plot-exterior nodes
figure(11);
hold all;
plot(x3_PEX,y3_PEX,'o','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
ylim([0 length_P]); xlim([offset-2*grid_size radius_P+2*grid_size]);
title({'Nodes at Pipe''s Exterior Surface';''},'FontWeight','normal');
grid on; grid minor; xlabel('r (m)'); ylabel('z (m)');
ax = gca; ax.FontSize = 14;

% Plot-central nodes tempeature
figure(12); 
hold all;
Tintrp3_PEX = interpolateTemperature(results,x3_PEX,y3_PEX);
p6 = plot(y3_PEX,Tintrp3_PEX,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
title({'Temperature at Pipe''s Exterior Surface';''},'FontWeight','normal');
grid on; grid minor; xlabel('z (m)'); ylabel('T (degC)');
ylim([75 77]); 
ax = gca; ax.FontSize = 14;

figure(13)
hold all;
[Tx3_PEX,Ty3_PEX] = evaluateTemperatureGradient(results,x3_PEX,y3_PEX);
p32 = plot(y3_PEX,Tx3_PEX,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);    
title({'Temperature Gradient at Pipe''s Exterior Surface';''},...
    'FontWeight','normal');
grid on; grid minor; xlabel('z (m)'); ylabel('dT/dr (degC/m)');
ylim([-1110 -1108]); 
ax = gca; ax.FontSize = 14;

figure(14)
hold all;
[Tx3_PEX,Ty3_PEX] = evaluateTemperatureGradient(results,x3_PEX,y3_PEX);
p33 = plot(y3_PEX,-k1*Tx3_PEX,'d','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);    
title({'Heat Flux at Pipe''s Exterior Surface';''},...
    'FontWeight','normal');
grid on; grid minor; xlabel('z (m)'); ylabel('q (W/m^2)');
ylim([505 515]); 
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(15); 
hold all;
Tintrp2_PEX = interpolateTemperature(results,x2_PEX,y2_PEX);
p2_PEX = plot(x2_PEX,Tintrp2_PEX,'>','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p2_Cu = plot(x2_Cu,Tintrp2_Cu,'d','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');

figure(16)
hold all;
[Tx2_PEX,Ty2_PEX] = evaluateTemperatureGradient(results,x2_PEX,y2_PEX);
p21_PEX = plot(x2_PEX,-k_PEX*Tx2_PEX,'>','MarkerSize',5,'MarkerFaceColor',[0.5 0.5 0.5]);
p22_Cu = plot(x2_Cu,-k_Cu*Tx2_Cu,'d','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');

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)
% 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


