-
Notifications
You must be signed in to change notification settings - Fork 13
/
thermal_sim.py
155 lines (119 loc) · 5.02 KB
/
thermal_sim.py
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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
4-step thermal simulator
Peter Attia
Last modified February 23, 2019
INPUTS:
-C1, C2, C3: C rate parameters (C4 is calculated)
OUTPUT: Arbitrary 'degradation' value
"""
import math
import numpy as np
def sim(C1, C2, C3):
# DETERMINE C4
C4 = 0.2/(1/6 - (0.2/C1 + 0.2/C2 + 0.2/C3))
# CONSTANT PARAMETERS
R = 0.009 # [=] m, radius of 18650 cylindrical cell
L = 0.065 # [=] m, length of 18650 cylindrical cell
R_int = 0.017 # [=] Ohms, internal resistance
Tinit = 30 # [=] deg C, initial temperature
Tinf = 30 # [=] deg C, environmental temperature
V = math.pi * R**2 * L # [=] 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 # [=] 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
# DEGRADATION PARAMETERS. ESTIMATED FROM SMITH, DAHN ET AL 2011
Ea = 0.122 # [=] eV, activation energy
kb = 8.617E-5 # [=] eV/K, Boltzmann's constant
# CALCULATED PARAMETERS
I1 = C1*1.1 # [=] A, C1 discharge current
I2 = C2*1.1 # [=] A, C2 discharge current
I3 = C3*1.1 # [=] A, C1 discharge current
I4 = C4*1.1 # [=] A, C2 discharge current
Qn = 1.1*20/100 # [=] Ah, capacity filled during one step
t1 = Qn / I1 * 3600 # [=] s, time of C1
t2 = Qn / I2 * 3600 # [=] s, time of C2
t3 = Qn / I3 * 3600 # [=] s, time of C3
t4 = Qn / I4 * 3600 # [=] s, time of C4
#total_time = t1 + t2 + t3 + t4; # [=] s, total time for both C1 and C2 steps
alpha = k/(rho*Cp) # [=] m^2/s, thermal diffusivity
power1 = I1**2 * R_int # [=] W, power generated by cell during C1
e_gen1 = power1 / V # [=] W/m^3, volumetric heat generation term
power2 = I2**2 * R_int # [=] W, power generated by cell during C2
e_gen2 = power2 / V # [=] W/m^3, volumetric heat generation term
power3 = I3**2 * R_int # [=] W, power generated by cell during C3
e_gen3 = power3 / V # [=] W/m^3, volumetric heat generation term
power4 = I4**2 * R_int # [=] W, power generated by cell during C4
e_gen4 = power4 / V # [=] W/m^3, volumetric heat generation term
""""
--------------------------BEGIN SIMULATION--------------------------
"""
## Set up problem
dr = 0.001 # [=] m, length step
dt = 20 # [=] s, time step. <2s for minimal error, 5s is reasonable
N = int(np.round(R/dr + 1))
r = np.arange(-dr,R+2*dr,dr)
T = Tinit*np.ones((N+2,1)) # initialize domain to be IC
#T[N+1,0] = T[N,0]-dr*k*h*(T[N,0]-Tinf) # Env. boundary condition
# PRE-INITIALIZE VARIABLES
time = 0
deg_rates = 0
def fin_el(Tin, e_gen):
# FINITE ELEMENT SIMULATION
mat = np.zeros((N+2,N+2))
rhs = np.zeros((N+2,1))
# Internal domain
for i in np.arange(2,N+1):
mat[i,i] = (r[i] * dr**2 * dt) * (1/(alpha*dt)+2/(dr**2))
mat[i,i+1] = (r[i] * dr**2 * dt) * (-1/(2*r[i]*dr)-1/(dr**2))
mat[i,i-1] = (r[i] * dr**2 * dt) * (1/(2*r[i]*dr)-1/(dr**2))
rhs[i] = (r[i] * dr**2 * dt) * ( Tin[i,0]/(alpha*dt)+e_gen/k )
# Boundary condition at r = 0
mat[0,0] = 1
mat[0,2] = -1
# Cartesian singularity at r = 0
mat[1,1] = (dr ** 2 * dt) * (4/(dr**2))
mat[1,2] = (dr ** 2 * dt) * ((1/dt - 4/(dr**2)))
rhs[1] = (dr ** 2 * dt) * (Tin[1,0]/(dt) + e_gen/k)
# Boundary condition at r = R
mat[N+1,N] = -h
mat[N+1,N+1] = -k/(2*dr)
mat[N+1,N-1] = k/(2*dr)
rhs[N+1] = -h*Tinf
return np.linalg.lstsq(mat,rhs)[0] # T(r)
## C1 step
while time < t1:
T = fin_el(T, e_gen1)
time = time + dt
for Tidx in T:
deg_rates = deg_rates + math.exp(-Ea/(kb*Tidx))
## C2 step
while time < t1 + t2:
T = fin_el(T, e_gen2)
time = time + dt
for Tidx in T:
deg_rates = deg_rates + math.exp(-Ea/(kb*Tidx))
## C3 step
while time < t1 + t2 + t3:
T = fin_el(T, e_gen3)
time = time + dt
for Tidx in T:
deg_rates = deg_rates + math.exp(-Ea/(kb*Tidx))
## C4 step
while time < t1 + t2 + t3 + t4:
T = fin_el(T, e_gen4)
time = time + dt
for Tidx in T:
deg_rates = deg_rates + math.exp(-Ea/(kb*Tidx))
# arbitrary factor to give lifetimes on the order of 100s of cycles
return deg_rates*1e15