% Based on Equation (4) from Forgez et al. 2010
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Thermal inputs
Rout = 6.651562529;
cp = 978.2969842;
Rin = 1.408390014;
m = 0.095;
CpRinRout = m*cp*(Rin+Rout);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the length of time of the experiment
tfin = length(texp);
% Calculate the open circuit voltage at each time, this is based on the
% state of charge of the cell
% Initially predefine a vector
OCVatt = zeros(tfin,1);
% Go through each time
for i = 1:tfin
% Find the location of the state of charge
[~, idx] = min(abs(SoC-(100*SoCexp(i))));
% Use this location to find the open circuit voltage
OCVatt(i) = OCV(idx);
end
% Calculate the average value for the ambient temperature
TaAverage = mean(Taexp);
% Calculate the heat generation at each time
%Qdot = abs(Iexp.*(Vexp-OCVatt));
Qdot = abs(Iexp.*(Vexp-OCVatt));;
% Predefine a vector for the surface temperature
Ts = zeros(tfin,1);
% The initial temperature is the same as the initial temperature from the
% experiment
Ts(1) = Tsexp(1);
E_heat = 0;
E_heat_exp = 0;
% Go through each time
for i = 2:tfin
% Calculate the surface temperature
Ts(i) = (((TaAverage-Ts(i-1))/CpRinRout + (Qdot(i)*Rout)/CpRinRout) * (texp(i) - texp(i-1))) + Ts(i-1);
E_heat = E_heat + (Ts(i) - TaAverage) / Rout;
E_heat_exp = E_heat_exp + (Tsexp(i) - TaAverage) / Rout;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CREATE PLOTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(0,'DefaultAxesFontName', 'Times New Roman')
set(0,'DefaultAxesFontSize', 12)
figure('units','centimeters','position',[10 10 8.8 5.0])
hold on
subplot(2,1,1);
p1 = plot(texp/60,Ts,'b',texp/60,Tsexp,'g');
xlabel('Time (mins)')
ylabel(sprintf('Temperature (%cC)', char(176)))
legend('Calculated Temperature','Measured Temperature','Location','northwest')
subplot(2,1,2);
p3 = plot(texp/60,abs(Tsexp-Ts));
%p3 = plot(texp/60,100*abs((Tsexp-Ts)./(Tsexp+273)));
p3.Color = [0 0 1.0000];
p3.LineStyle = '-';
p3.LineWidth = 1;
p3.Marker = 'none';
xlabel('Time (mins)')
ylabel(sprintf('Absolute Error (%cC)', char(176)))
%ylabel(sprintf('Percentage Error (%)', char(176)))
%legend('Calculated Temperature','Measured Temperature','northwest')