-
Notifications
You must be signed in to change notification settings - Fork 2
/
battery.m
127 lines (101 loc) · 3.5 KB
/
battery.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
function [] = battery(C_rate)
%% Thermal model of heating an 18650 battery during use
% Peter Attia, November 8, 2017
% Reference: http://cecs.wright.edu/~sthomas/htchapter02.pdf
% Equation 2-26
% This simulation requires the PDE toolbox.
% For help with the PDE solver ('pdepe'), see:
% https://www.mathworks.com/help/matlab/ref/pdepe.html
%% Assumptions
% * Long cylinder (R/L = 0.14)
% * Heat transfer through SS can is fast (for can, T_in = T_out)
%% Initial condition
% * T(t=0,r) = T_init
%% Boundary conditions
% (1) dT/dr(t,r=0) = 0 (symmetry at cylinder centerline)
% (2) -k(dT/dr)|(r=R) = h(T(t=0,R) - T_env) (convection)
% NOTE: BC (1) reduces to constant temperature case (T(t,x=R) = T_init)
% if h is large
%%%%%%%%%%%%%%%%%%%%%%%%
%% PARAMETERS
global k e_gen Tinit h alpha
%% RELATIVELY CONSTANT PARAMETERS
cell_capacity = 1.1; % [=] Ah
I = cell_capacity*C_rate; % [=] A, current (discharge)
R_int = 0.017; % [=] Ohms, internal resistance
Tinit = 30; % [=] deg C, initial temperature
V = pi*(0.009)^2*(0.065); % [=] m^3, cell volume
%% ESTIMATED PARAMETERS
% Values for k, Cp, and rho from Drake et al, JPS (2014)
% http://www.uta.edu/faculty/jaina/MTL/pubs/Drake-JPS2014.pdf
k = 0.20; % [=] W/m-K, thermal conductivity
Cp = 1000;%1720; % [=] J/kg-K, specific heat capacity
rho = 2362; % [=] kg/m^3, density
% Value for rho estimated from mass of cell is in good agreement with above
% value
% http://www.a123batteries.com/v/vspfiles/images/pdf/APR18650M1A.pdf
% rho = (39/1000)/V; % [=] kg/m^3, average density
% Value for h taken from Engineering Toolbox: h = 10 for air, 500 for oil
h = 10; % [=] W/m^2-K, heat transfer coefficient
%% CALCULATED PARAMETERS
alpha = k/(rho*Cp);
power = I^2*R_int; % [=] W, power generated by cell
e_gen = power/V; % [=] W/m^3, volumetric heat generation term
total_t = 3600/C_rate; % [=] s, total time
%%%%%%%%%%%%%%%%%%%%%%%%
close all
filename = ['battery_' num2str(C_rate) 'C'];
m = 1; % cylinder geometry
r = linspace(0,0.009,100); % r = 0 to 0.009 m (9 mm)
t = linspace(0,total_t,100); % t = 0 to total_t s
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,r,t);
% Extract the first solution component as T
T = sol(:,:,1);
%% Surface plot
pcolor(r.*1000,t./60,T)
xlabel('Distance from center, r (mm)')
ylabel('Time, t (min)')
hcb = colorbar; title(hcb,'Temperature (^oC)')
title([num2str(C_rate) 'C'])
print(filename,'-dsvg')
%% Animated plot of T(r,t)
figure
SOC = t.*C_rate/3600*100;
y_limits = [29 ceil(max(T(:))+1)];
for i = 1:length(T(:,1))
plot(r.*1000,T(i,:));
title(['t = ' num2str(t(i)./60,3) ' min; SOC = ' num2str(SOC(i),3) '%'],'FontSize',16)
xlabel('Distance from center, r (mm)')
ylabel('Temperature (^oC)')
ylim(y_limits)
% Create GIF
drawnow
frame = getframe(2);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
if i == 1
imwrite(imind,cm,[filename '.gif'],'gif','Loopcount',1);
else
imwrite(imind,cm,[filename '.gif'],'gif','WriteMode','append','DelayTime',0.05);
end
end
% --------------------------------------------------------------------------
function [c,f,s] = pdex1pde(r,t,u,DuDx)
% Main PDE
global alpha k e_gen
c = 1/alpha;
f = DuDx;
s = e_gen/k;
% --------------------------------------------------------------------------
function u0 = pdex1ic(r)
% Initial condition
global Tinit
u0 = Tinit;
% --------------------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
% Boundary conditions
global h k Tinit
pl = 0;
ql = 1;
pr = h*(ur-Tinit);
qr = k;