%% Housekeeping
clear;
close all;
clc;
format short
%% Given Data:

del_T = 10;
Phi = linspace(0,1);                                                       % Volume fraction of metal
Phif = 0.75;                                                               % Volume fraction of fields metal

% Fields metal:

k_f = 18.75;                                                               % Thermal conductivity of fields metal in W/mK
h_f = 25400;                                                               % Latent Heat of fields matal in J/kg
rho_f = 7880;                                                              % Density of fields metal in Kg/m^3
cp_f = 325;                                                                % Specific heat capacity of fields metal in J/KgK
Cp_f = cp_f*rho_f;                                                         % Specific heat capacity of fields metal in J/m^3K
H_f = h_f*rho_f;                                                           % Latent heat of fields metal in J/m^3

% Paraffin:

k_p = 0.2;                                                                 % Thermal conductivity of paraffin in W/mK
h_p = 243500;                                                              % Latent Heat of paraffin in J/kg
rho_p = 774;                                                               % Density of paraffin in Kg/m^3
cp_p = 2160;                                                               % Specific heat capacity of paraffin in J/KgK
Cp_p = cp_p*rho_p;                                                         % Specific heat capacity of paraffin in J/m^3K
H_p = h_p*rho_p;                                                           % Latent heat of paraffin in J/m^3

% Cu:

k_Cu = 385;                                                                % Thermal conductivity of copper in W/mK
rho_Cu = 8960;                                                             % Density of copper in Kg/m^3
cp_Cu= 390;                                                                % Specific heat capacity of copper in J/KgK
Cp_Cu = cp_Cu*rho_Cu;                                                      % Specific heat capacity of copper in J/m^3K

%Al:

k_Al = 205;                                                                % Thermal conductivity of Aluminum in W/mK
rho_Al = 2700;                                                             % Density of Aluminum in kg/m^3
cp_Al = 900;                                                               % Specific heat capacity of Aluminum in J/kgK
Cp_Al = cp_Al*rho_Al;                                                      % Specific heat capacity of Aluminum in J/m^3K

% GCF:

k_GCF = 600;                                                               % Thermal conductivity of copper in W/mK
rho_GCF = 2266;                                                            % Density of copper in Kg/m^3
cp_GCF = 720;                                                              % Specific heat capacity of GCF in J/kgK
Cp_GCF = cp_GCF*rho_GCF;                                                   % Specific heat capacity of GCF in J/m^3K

%Paraffin + Fields metal:

h_pcm = h_f.*Phif + (h_p.*(1-Phif));                                       % Effective Latent heat of pcm in J/kg
cp_pcm = cp_f.*Phif + (cp_p.*(1-Phif));                                    % Effective Specific heat capacity of pcm in J/kgK
rho_pcm = rho_f.*Phif + (rho_p.*(1-Phif));                                 % Effective density of pcm in Kg/m^3
H_pcm = h_pcm.*rho_pcm;                                                    % Effective Latent heat of pcm in J/m^3
Cp_pcm = cp_pcm.*rho_pcm;                                                  % Effective Specific heat capacity of pcm in J/m^3K

%% Setting up: 

% Thermal conductivity:
k1 = k_f.*Phif + (k_p.*(1-Phif));                                          % Effective Thermal conductivity of pcm in W/mK
k2 = k_Cu.*Phi + (k1.*(1-Phi));                                            % Effective Thermal conductivity of Cu+pcm(k1) W/mK
k3 = k_Cu.*Phi + (k_p.*(1-Phi));                                           % Effective Thermal conductivity of Cu+paraffin W/mK
k4 = k_Al.*Phi + (k_p.*(1-Phi));                                           % Effective Thermal conductivity of Al+paraffin in W/mK
k5 = k_GCF.*Phi + (k_p.*(1-Phi));                                          % Effective Thermal conductivity of GCF+paraffin in W/mK
k6 = k_GCF.*Phi + (k1.*(1-Phi));                                           % Effective Thermal conductivity of GCF+pcm in W/mK
k7 = k_Al.*Phi + (k1.*(1-Phi));                                            % Effective Thermal conductivity of Al+pcm in W/mK

% Energy density:
E_eff_1 = Phi*0 + ((H_pcm + (Cp_pcm*del_T)))*0.000001;            % Effective Energy density of pcm in J/cm^3
E_eff_2 = ((((Cp_Cu*del_T)).*Phi) + ((H_pcm + (Cp_pcm*del_T)).*(1-Phi)))*0.000001;  % Effective Energy density of Cu+pcm in J/cm^3
E_eff_3 = ((((Cp_Cu*del_T)).*Phi) + ((H_p + (Cp_p*del_T)).*(1-Phi)))*0.000001;      % Effective Energy density of Cu+paraffin in J/cm^3
E_eff_4 = ((((Cp_Al*del_T)).*Phi) + ((H_p + (Cp_p*del_T)).*(1-Phi)))*0.000001;      % Effective Energy density of Al+paraffin in J/cm^3
E_eff_5 = ((((Cp_GCF*del_T)).*Phi) + ((H_p + (Cp_p*del_T)).*(1-Phi)))*0.000001;     % Effective Energy density of GCF+paraffin in J/cm^3
E_eff_6 = ((((Cp_GCF*del_T)).*Phi) + ((H_pcm + (Cp_pcm*del_T)).*(1-Phi)))*0.000001;   % Effective Energy density of GCF+pcm in J/cm^3
E_eff_7 = ((((Cp_Al*del_T)).*Phi) + ((H_pcm + (Cp_pcm*del_T)).*(1-Phi)))*0.000001;    % Effective Energy density of Al+pcm in J/cm^3

X1 = times(k1,E_eff_1/0.000001);
X2 = times(k2,E_eff_2/0.000001);
X3 = times(k3,E_eff_3/0.000001);
X4 = times(k4,E_eff_4/0.000001);
X5 = times(k5,E_eff_5/0.000001);
X6 = times(k6,E_eff_6/0.000001);
X7 = times(k7,E_eff_7/0.000001);

Eta_eff1 = sqrt(X1);                                                         % Figure of merit of pcm in (Wsec^0.5)/Km^2
Eta_eff2 = sqrt(X2);                                                         % Figure of merit of Cu+pcm in (Wsec^0.5)/Km^2
Eta_eff3 = sqrt(X3);                                                         % Figure of merit of Cu+paraffin in (Wsec^0.5)/Km^2
Eta_eff4 = sqrt(X4);                                                         % Figure of merit of Al+paraffin in (Wsec^0.5)/Km^2
Eta_eff5 = sqrt(X5);                                                         % Figure of merit of GCF+paraffin in (Wsec^0.5)/Km^2
Eta_eff6 = sqrt(X6);                                                         % Figure of merit of GCF+pcm in (Wsec^0.5)/Km^2
Eta_eff7 = sqrt(X7);                                                         % Figure of merit of Al+pcm in (Wsec^0.5)/Km^2

%% Plot

colors = [60/255 133/255 244/255; 225/255 68/255 55/255; 244/255 160/255 0; 15/255 157/255 88/255; 144/255 103/255 167/255; 128/255 133/255 133/255];
set(0, 'DefaultLineLineWidth', 1.5);
set(groot, 'defaultAxesColorOrder', colors);
 
figure(1); clf; grid on;
hold on;
y = [Eta_eff1; Eta_eff2; Eta_eff3; Eta_eff4; Eta_eff5; Eta_eff6; Eta_eff7];

%plot(Phi,y)
plot(Phi(1,:),y(1,:),'k')
hold on;
plot(Phi(1,:),y(2,:))
hold on
plot(Phi(1,:),y(3,:),'--')
hold on
plot(Phi(1,:),y(4,:),'--')
hold on
plot(Phi(1,:),y(5,:),'--')
hold on
plot(Phi(1,:),y(6,:))
hold on
plot(Phi(1,:),y(7,:))

title('Volume fraction of FM = 0.75')
xlabel('Volume fraction of metal');
ylabel('Figure of Merit [(Wsec^0.5)/Km^2]');
ltag = {'FM+P','Cu+FM+P','Cu+P','Al+P','GCF+P','GCF+FM+P','Al+FM+P'};
legend(ltag,'Location','Northwest');
hold off;