-
Notifications
You must be signed in to change notification settings - Fork 1
/
li_s_battery_model.py
711 lines (564 loc) · 27.6 KB
/
li_s_battery_model.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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
# -*- coding: utf-8 -*-
"""
Created on Thu May 9 12:30:28 2019
@author: dkorff
This module contains the set up and running of the simulation for the Li-S
model.
"""
import numpy as np
import pandas as pd
import time
import importlib
import cantera as ct
from matplotlib import pyplot as plt
import matplotlib
from assimulo.solvers import IDA
from assimulo.problem import Implicit_Problem
from assimulo.exception import TerminateSimulation
from li_s_battery_inputs import inputs
from li_s_battery_init import anode as an
from li_s_battery_init import sep
from li_s_battery_init import cathode as cat
from li_s_battery_init import sol_init
from li_s_battery_post import label_columns
from li_s_battery_post import tag_strings
from li_s_battery_post import plot_sim
from li_s_battery_post import plot_meanPS
def main():
res_class = eval(inputs.test_type)
# plt.close('all')
t_count = time.time()
SV_0 = sol_init.SV_0
SV_dot_0 = np.zeros_like(SV_0)
t_0 = 0.
t_f = 3600./inputs.C_rate #63006.69900049 93633
algvar = sol_init.algvar
atol = np.ones_like(SV_0)*1e-6
atol[cat.ptr_vec['eps_S8']] = 1e-15
atol[cat.ptr_vec['eps_Li2S']] = 1e-15
atol[cat.ptr_vec['rho_k_el']] = 1e-16 # 1e-16 for Bessler
rtol = 1e-6; sim_output = 50
rtol_ch = 1e-6
atol_ch = np.ones_like(SV_0)*1e-6
atol_ch[cat.ptr_vec['eps_S8']] = 1e-15
atol_ch[cat.ptr_vec['eps_Li2S']] = 1e-15
atol_ch[cat.ptr_vec['rho_k_el']] = 1e-30
rate_tag = str(inputs.C_rate)+"C"
# if 'cascade' in inputs.ctifile:
# ncols = 2 + inputs.flag_req
# else:
ncols = 1 + inputs.flag_req
fig, axes = plt.subplots(sharey="row", figsize=(9,12), nrows=3, ncols = (ncols)*inputs.n_cycles, num=1)
plt.subplots_adjust(wspace = 0.15, hspace = 0.4)
fig.text(0.35, 0.85, rate_tag, fontsize=20, bbox=dict(facecolor='white', alpha = 0.5))
# Set up user function to build figures based on inputs
"----------Equilibration----------"
print('\nEquilibrating...')
# Set external current to 0 for equilibration
cat.set_i_ext(0)
# cat.nucleation_flag = 1
# Create problem object
bat_eq = res_class(res_class.res_fun, SV_0, SV_dot_0, t_0)
bat_eq.external_event_detection = True
bat_eq.algvar = algvar
# Create simulation object
sim_eq = IDA(bat_eq)
sim_eq.atol = atol
sim_eq.rtol = rtol
sim_eq.verbosity = sim_output
sim_eq.make_consistent('IDA_YA_YDP_INIT')
t_eq, SV_eq, SV_dot_eq = sim_eq.simulate(t_f)
# Put solution into pandas dataframe with labeled columns
SV_eq_df = label_columns(t_eq, SV_eq, an.npoints, sep.npoints, cat.npoints)
# Obtain tag strings for dataframe columns
tags = tag_strings(SV_eq_df)
# plot_sim(tags, SV_eq_df, 'Equilibrating', 0, fig, axes)
t_equilibrate = time.time() - t_count
print("Equilibration time = ", t_equilibrate, '\n')
print('Done equilibrating\n')
for cycle_num in np.arange(0, inputs.n_cycles):
"------------Discharging-------------"
print('Discharging...')
# cat.np_L = inputs.np_Li2S_init
cat.nucleation_flag = np.zeros((inputs.npoints_cathode, 1))
# New initial conditions from previous simulation
if cycle_num == 0:
# SV_0 = SV_0
SV_0 = SV_eq[-1, :]
# SV_dot_0 = SV_dot_0
SV_dot_0 = SV_dot_eq[-1, :]
else:
SV_0 = SV_0_cycle
SV_dot_0 = SV_dot_0_cycle
# Set external current
cat.set_i_ext(cat.i_ext_amp)
# Update problem instance initial conditions
bat_dch = res_class(res_class.res_fun, SV_0, SV_dot_0, t_0)
bat_dch.external_event_detection = True
bat_dch.algvar = algvar
# Re-initialize simulation object
sim_dch = IDA(bat_dch)
sim_dch.atol = atol
sim_dch.rtol = rtol
# sim_dch.maxh = 57
sim_dch.verbosity = sim_output
sim_dch.make_consistent('IDA_YA_YDP_INIT')
t_dch, SV_dch, SV_dot_dch = sim_dch.simulate(131865.32)
SV_dch_df = label_columns(t_dch, SV_dch, an.npoints, sep.npoints, cat.npoints)
# Obtain tag strings for dataframe columns
tags = tag_strings(SV_dch_df)
plot_sim(tags, SV_dch_df, 'Discharging', 0+2*cycle_num, fig, axes)
# plot_meanPS(SV_dch_df, tags, 'Discharging')
print('Done Discharging\n')
"--------Re-equilibration---------"
if inputs.flag_req == 1:
print('Re-equilibrating...')
# New initial conditions from previous simulation
SV_0 = SV_dch[-1, :]
SV_dot_0 = SV_dot_dch[-1, :]
# Set external current
cat.set_i_ext(0)
# Update problem instance initial conditions
bat_req = res_class(res_class.res_fun, SV_0, SV_dot_0, t_0)
bat_req.external_event_detection = True
bat_req.algvar = algvar
# Re-initialize simulation object
sim_req = IDA(bat_req)
sim_req.atol = atol
sim_req.rtol = rtol
sim_req.verbosity = sim_output
sim_req.make_consistent('IDA_YA_YDP_INIT')
t_req, SV_req, SV_dot_req = sim_req.simulate(t_f)
SV_req_df = label_columns(t_req, SV_req, an.npoints, sep.npoints, cat.npoints)
# plot_sim(tags, SV_req_df, 'Re-Equilibrating', 1, fig, axes)
print('Done re-equilibrating\n')
else:
SV_req = SV_dch
SV_dot_req = SV_dot_dch
"-----------Charging-----------"
if 'dog' in inputs.ctifile:
print('Charging...')
SV_0 = SV_req[-1, :]
SV_dot_0 = SV_dot_req[-1, :]
SV_0[cat.ptr_vec['eps_S8']] = cat.eps_cutoff
cat.set_i_ext(-cat.i_ext_amp)
# Update problem instance initial conditions
bat_ch = res_class(res_class.res_fun, SV_0, SV_dot_0, t_0)
bat_ch.external_event_detection = True
bat_ch.algvar = algvar
# Re-initialize simulation object
sim_ch = IDA(bat_ch)
sim_ch.atol = atol_ch
sim_ch.rtol = rtol_ch
if cycle_num > 0:
sim_ch.maxh = 0.1
sim_ch.verbosity = sim_output
sim_ch.make_consistent('IDA_YA_YDP_INIT')
t_ch, SV_ch, SV_dot_ch = sim_ch.simulate(t_f)
SV_ch_df = label_columns(t_ch, SV_ch, an.npoints, sep.npoints, cat.npoints)
plot_sim(tags, SV_ch_df, 'Charging', 1+inputs.flag_req+2*cycle_num, fig, axes)
plot_meanPS(SV_ch_df, tags, 'Charging')
SV_ch_df = 0
#
# print('Max S_8(e) concentration = ', max(SV_ch[:, 6]))
# SV_0_cycle = SV_ch[-1, :]
# SV_dot_0_cycle = SV_ch[-1, :]
print('Done Charging\n')
exp_data_01C = pd.read_csv(r'0.1C Data.csv', header=None)
exp_data_05C = pd.read_csv(r'0.5C Data.csv', header=None)
exp_data_1C = pd.read_csv(r'1C Data.csv', header=None)
Bessler = pd.read_csv(r'Bessler Dennis Data.csv', header=None)
SV_copy = SV_dch_df.copy()
SV_copy.loc[:, 'Time'] *= -cat.i_ext_amp*inputs.A_cat/3600/(cat.m_S_tot_0)
"Set up your figure"
fig = plt.figure(3)
ax = fig.add_axes([0.2,0.2,0.6,0.75])
fig.set_size_inches((10.,5.0))
# ax2 = ax.twinx()
"Formatting for the figure:"
fs = 20 #font size for plots
lw = 3.0 #line width for plots
# font = plt.matplotlib.font_manager.FontProperties(family='Times New Roman',size=fs-1)
for tick in ax.xaxis.get_major_ticks():
tick.label1.set_fontsize(fs)
tick.label1.set_fontname('Times New Roman')
for tick in ax.yaxis.get_major_ticks():
tick.label1.set_fontsize(fs)
tick.label1.set_fontname('Times New Roman')
color = matplotlib.cm.plasma(inputs.C_rate)
p1, = plt.plot(SV_copy.loc[:, 'Time'], SV_copy.loc[:, 'Phi_ed1'], 'k-', linewidth=lw)
# p2, = plt.plot(exp_data_01C.iloc[:,0], exp_data_01C.iloc[:,1], 'ro')
# p3, = plt.plot(exp_data_05C.iloc[:,0], exp_data_05C.iloc[:,1], 'co')
# p4, = plt.plot(exp_data_1C.iloc[:,0], exp_data_1C.iloc[:,1], 'ko')
p5, = plt.plot(Bessler.iloc[:,0], Bessler.iloc[:,1], 'mo', ms=4)
plt.xlim((0, 1700))
plt.xticks([0, 400, 800, 1200, 1600])
plt.ylim((1.8, 2.6))
plt.ylabel(r'Cell Voltage $[\mathrm{V}]$', fontstyle='normal', fontname='Times new Roman', fontsize=fs+2, labelpad=5.0)
plt.xlabel(r'Capacity $[\mathrm{Ah} \hspace{0.5} \mathrm{kg}^{-1}_{\mathrm{sulfur}}]$', fontstyle='normal', fontname='Times new Roman', fontsize=fs+2, labelpad=5.0)
# plt.legend(["Discharge", "Charge"])
file_name_dch = 'dch'+str(inputs.C_rate)+"C_"+inputs.mech+'.csv'
SV_dch = SV_dch_df.copy()
SV_dch.loc[:, 'Time'] *= -cat.i_ext_amp*inputs.A_cat/3600/(cat.m_S_0 + cat.m_S_el)
SV_dch.to_csv(file_name_dch, index=False, header=True)
t_elapsed = time.time() - t_count
print('t_cpu=', t_elapsed, '\n')
return SV_eq_df, SV_dch_df, SV_ch_df, tags
"============================================================================="
"===========RESIDUAL CLASSES AND HELPER FUNCTIONS BEYOND THIS POINT==========="
"============================================================================="
from li_s_battery_init import sulfur_obj as sulfur
from li_s_battery_init import Li2S_obj as Li2S
from li_s_battery_init import carbon_obj as carbon
from li_s_battery_init import sulfur_el_s as S_el_s
from li_s_battery_init import Li2S_el_s as L_el_s
from li_s_battery_init import carbon_el_s as C_el_s
from li_s_battery_init import Li2S_tpb
from li_s_battery_init import lithium_obj as lithium
from li_s_battery_init import lithium_el_s as lithium_s
from li_s_battery_init import conductor_obj as conductor
from li_s_battery_init import elyte_obj as elyte
from li_s_battery_functions import set_state, set_geom, set_rxn
from li_s_battery_functions import set_state_sep
from li_s_battery_functions import dst
from math import pi, exp, tanh
class cc_cycling(Implicit_Problem):
def res_fun(t, SV, SV_dot):
res = np.zeros_like(SV)
ptr = cat.ptr; F = ct.faraday; R = ct.gas_constant; T = inputs.T
"""=============================CATHODE============================="""
"""CC BOUNDARY"""
j = 0; offset = cat.offsets[int(j)]
i_ext = cat.get_i_ext()
s2 = set_state(SV, offset, cat.ptr)
# Set electronic current and ionic current boundary conditions
i_el_p = i_ext
i_io_p = 0
N_io_p = 0
"""=============================CATHODE============================="""
"""INTERIOR NODES"""
for j in np.arange(1, cat.npoints):
# Set previous outlet fluxes to new inlet fluxes
i_el_m = i_el_p
i_io_m = i_io_p
N_io_m = N_io_p
s1 = dict(s2)
# Update offset to NEXT node
offset = cat.offsets[int(j)]
s2 = set_state(SV, offset, cat.ptr)
# Set variables to CURRENT NODE value
offset = cat.offsets[int(j-1)]
eps_S8 = max(SV[offset + ptr['eps_S8']], cat.eps_cutoff)
eps_Li2S = max(SV[offset + ptr['eps_Li2S']], cat.eps_cutoff)
eps_el = 1 - cat.eps_C_0 - eps_S8 - eps_Li2S
# Set states for THIS node
carbon.electric_potential = s1['phi_ed']
elyte.electric_potential = s1['phi_el']
conductor.electric_potential = s1['phi_ed']
elyte.X = s1['X_k']
D_el = cat.D_el*eps_el**(cat.bruggeman)
# Current node plus face boundary fluxes
i_el_p = cat.sigma_eff*(s1['phi_ed'] - s2['phi_ed'])*cat.dyInv
N_io_p, i_io_p = dst(s1, s2, D_el, cat.dy, cat.dy)
sdot_C = C_el_s.get_net_production_rates(elyte)
mult = tanh(eps_S8/cat.eps_dropoff)
sdot_S8 = S_el_s.get_creation_rates(sulfur) - mult*S_el_s.get_destruction_rates(sulfur)
sdot_S = S_el_s.get_net_production_rates(elyte)
mult = tanh(eps_Li2S/cat.eps_dropoff)
sdot_Li2S = L_el_s.get_creation_rates(Li2S) - mult*(L_el_s.get_destruction_rates(Li2S))
sdot_L = L_el_s.get_net_production_rates(elyte)
sdot_tpb = Li2S_tpb.get_creation_rates(Li2S) - mult*(Li2S_tpb.get_destruction_rates(Li2S))
sdot_tpb_el = mult*Li2S_tpb.get_creation_rates(elyte) - Li2S_tpb.get_destruction_rates(elyte)
# default_np_L = 1e5
# nucleation_coeff = 8e12
# # Calculate new particle radii based on new volume fractions
# C_Li2S2 = SV[offset + ptr['rho_k_el'][-1]]
# nucl_mult = tanh(C_Li2S2/cat.nucl_thresh)
# m_S_el = inputs.A_cat*eps_el*cat.dy*np.dot(cat.W_S_k, SV[offset + ptr['rho_k_el']])
# m_S = sulfur.density_mass*eps_S8*cat.dy + m_S_el
# local_C_rate = -i_el_m/1675/m_S
# if C_Li2S2 > 1e-15 and cat.nucleation_flag[j-1] == 0:
# np_L = nucleation_coeff*exp(1.7953*local_C_rate)
# cat.np_L[j-1] = np_L
# cat.nucleation_flag[j-1] = 1
# print(np_L, local_C_rate, -i_el_m)
# elif cat.nucleation_flag[j-1] == 1:
# np_L = nucl_mult*cat.np_L[j-1]
# else:
# np_L = default_np_L
np_S = inputs.np_S8_init
np_L = inputs.np_Li2S_init
A_S = cat.A_S_0*(eps_S8/cat.eps_S_0)**1.5
# A_S = 2*pi*np_S*(3*eps_S8/2/np_S/pi)**(2/3)
A_L = cat.A_L_0*(eps_Li2S/cat.eps_L_0)**1.5
# A_L = 2*pi*np_L*(3*eps_Li2S/2/np_L/pi)**(2/3)
r_S = 3*eps_S8/A_S
r_L = 3*eps_Li2S/A_L
tpb_len = 3*eps_Li2S/(r_L**2)
A_C = cat.A_C_0 #- (pi*np_S*r_S**2) - (pi*np_L*r_L**2)
# print(A_C, SV[offset + ptr['phi_ed']])
R_C = sdot_C*A_C
if eps_S8 < 1e-5:
R_S = 0*sdot_S*A_S
sw = 0
else:
R_S = sdot_S*A_S
sw = 1
if eps_Li2S < 1e-5 and i_ext == 0:
R_L = 0*sdot_L*A_L + 0*sdot_tpb_el*tpb_len
sw2 = 0
else:
R_L = sdot_L*A_L + sdot_tpb_el*tpb_len
sw2 = 1
i_C = (C_el_s.get_net_production_rates(conductor)*A_C +
(Li2S_tpb.get_creation_rates(conductor)*mult -
Li2S_tpb.get_destruction_rates(conductor))*tpb_len)
i_Far = (i_C)*F/cat.dyInv
# if eps_S8 < 1e-3:
# print(-i_Far + i_el_m - i_el_p)
# Net rate of formation
R_net = R_C + R_S + R_L
R_net[cat.ptr['iFar']] += (-i_Far + i_el_m - i_el_p)/cat.dy/F
"""Calculate change in Sulfur"""
res[offset + ptr['eps_S8']] = (SV_dot[offset + ptr['eps_S8']]
- sw*sulfur.volume_mole*sdot_S8*A_S)
"""Calculate change in Li2S"""
res[offset + ptr['eps_Li2S']] = (SV_dot[offset + ptr['eps_Li2S']]
- sw2*Li2S.volume_mole*(sdot_Li2S*A_L
+ sdot_tpb*tpb_len))
# print(sdot_Li2S, A_L, sdot_tpb, tpb_len, res[offset + ptr['eps_Li2S']], 'node =', j)
"""Calculate change in electrolyte"""
res[offset + ptr['rho_k_el']] = (SV_dot[offset + ptr['rho_k_el']] -
(R_net + (N_io_m - N_io_p)*cat.dyInv)/eps_el
+ SV[offset + ptr['rho_k_el']]*(- SV_dot[offset + ptr['eps_S8']]
- SV_dot[offset + ptr['eps_Li2S']])/eps_el)
# print(R_net, i_ext)
"""Calculate change in delta-phi double layer"""
res[offset + ptr['phi_dl']] = (SV_dot[offset + ptr['phi_dl']] -
(-i_Far + i_el_m - i_el_p)*cat.dyInv/cat.C_dl/A_C)
# print(i_Far, i_ext, eps_S8)
"""Algebraic expression for charge neutrality in all phases"""
res[offset + ptr['phi_ed']] = i_el_m - i_el_p + i_io_m - i_io_p
"""Calculate change in S8 nucleation sites"""
res[offset + ptr['np_S8']] = SV[offset + ptr['np_S8']] - np_S
"""Calculate change in Li2S nucleation sites"""
res[offset + ptr['np_Li2S']] = SV[offset + ptr['np_Li2S']] - np_L
"""=============================CATHODE============================="""
"""SEPARATOR BOUNDARY"""
i_el_m = i_el_p
i_io_m = i_io_p
N_io_m = N_io_p
s1 = dict(s2)
# Shift forward to NEXT node, first separator node (j=0)
j = 0; offset = sep.offsets[int(j)]
s2 = set_state_sep(SV, offset, sep.ptr)
# Shift back to THIS node, set THIS node outlet conditions
j = cat.npoints-1; offset = cat.offsets[int(j)]
# Set variables to CURRENT NODE value
# geom = set_geom(SV, offset, cat.ptr)
eps_S8 = max(SV[offset + ptr['eps_S8']], cat.eps_cutoff)
eps_Li2S = max(SV[offset + ptr['eps_Li2S']], cat.eps_cutoff)
eps_el = 1 - cat.eps_C_0 - eps_S8 - eps_Li2S
carbon.electric_potential = s1['phi_ed']
elyte.electric_potential = s1['phi_el']
conductor.electric_potential = s1['phi_ed']
elyte.X = s1['X_k']
# Set outlet boundary conditions for THIS node
i_el_p = 0
D_el = cat.D_el*eps_el**(cat.bruggeman)
N_io_p, i_io_p = dst(s1, s2, D_el, cat.dy, sep.dy)
sdot_C = C_el_s.get_net_production_rates(elyte)
mult = tanh(eps_S8/cat.eps_dropoff)
sdot_S8 = S_el_s.get_creation_rates(sulfur) - mult*S_el_s.get_destruction_rates(sulfur)
sdot_S = S_el_s.get_net_production_rates(elyte)
mult = tanh(eps_Li2S/cat.eps_dropoff)
sdot_Li2S = L_el_s.get_creation_rates(Li2S) - mult*(L_el_s.get_destruction_rates(Li2S))
sdot_L = L_el_s.get_net_production_rates(elyte)
sdot_tpb = Li2S_tpb.get_creation_rates(Li2S) - mult*(Li2S_tpb.get_destruction_rates(Li2S))
sdot_tpb_el = mult*Li2S_tpb.get_creation_rates(elyte) - Li2S_tpb.get_destruction_rates(elyte)
# print(L_el_s.delta_gibbs, i_ext, t)
# Calculate new particle radii based on new volume fractions
# C_Li2S2 = SV[offset + ptr['rho_k_el'][-1]]
# nucl_mult = tanh(C_Li2S2/cat.nucl_thresh)
# m_S_el = inputs.A_cat*eps_el*cat.dy*np.dot(cat.W_S_k, SV[offset + ptr['rho_k_el']])
# m_S = sulfur.density_mass*eps_S8*cat.dy + m_S_el
# local_C_rate = -i_el_m/1675/m_S
# if C_Li2S2 > 1e-15 and cat.nucleation_flag[-1] == 0:
# np_L = nucleation_coeff*exp(1.7953*local_C_rate)
# cat.np_L[-1] = np_L
# cat.nucleation_flag[-1] = 1
# print(np_L, t, '\n')
# elif cat.nucleation_flag[-1] == 1:
# np_L = nucl_mult*cat.np_L[-1]
# else:
# np_L = default_np_L
np_S = inputs.np_S8_init
np_L = inputs.np_Li2S_init #SV[offset + ptr['np_Li2S']]
A_S = cat.A_S_0*(eps_S8/cat.eps_S_0)**1.5
# A_S = 2*pi*np_S*(3*eps_S8/2/np_S/pi)**(2/3)
A_L = cat.A_L_0*(eps_Li2S/cat.eps_L_0)**1.5
# A_L = 2*pi*np_L*(3*eps_Li2S/2/np_L/pi)**(2/3)
r_S = 3*eps_S8/A_S
r_L = 3*eps_Li2S/A_L
tpb_len = 3*eps_Li2S/(r_L**2)
A_C = cat.A_C_0 #- (pi*np_S*r_S**2) - (pi*np_L*r_L**2)
# print(A_C, SV[offset + ptr['phi_ed']], '\n\n')
# cat.A_C_vec = np.append(cat.A_C_vec, A_C)
R_C = sdot_C*A_C
if eps_S8 < 1e-5:
R_S = 0*sdot_S*A_S
sw = 0
else:
R_S = sdot_S*A_S
sw = 1
if eps_Li2S < 1e-5 and i_ext == 0:
R_L = 0*sdot_L*A_L + 0*sdot_tpb_el*tpb_len
sw2 = 0
else:
R_L = sdot_L*A_L + sdot_tpb_el*tpb_len
sw2 = 1
# print(eps_S8, '\n', R_S, '\n')
i_C = (C_el_s.get_net_production_rates(conductor)*A_C +
(Li2S_tpb.get_creation_rates(conductor)*mult -
Li2S_tpb.get_destruction_rates(conductor))*tpb_len)
i_Far = (i_C)*F/cat.dyInv
# if eps_S8 < 1e-3:
# print(-i_Far + i_el_m - i_el_p,
# -t*cat.i_ext_amp*inputs.A_cat/3600/(cat.m_S_0 + cat.m_S_el), '\n\n')
# Net rate of formation
R_net = R_C + R_S + R_L
R_net[cat.ptr['iFar']] += (-i_Far + i_el_m - i_el_p)/cat.dy/F
"""Calculate change in Sulfur"""
res[offset + ptr['eps_S8']] = (SV_dot[offset + ptr['eps_S8']]
- sw*sulfur.volume_mole*sdot_S8*A_S)
"""Calculate change in Li2S"""
res[offset + ptr['eps_Li2S']] = (SV_dot[offset + ptr['eps_Li2S']]
- sw2*Li2S.volume_mole*(sdot_Li2S*A_L
+ sdot_tpb*tpb_len))
# print(sdot_Li2S, A_L, sdot_tpb, tpb_len, res[offset + ptr['eps_Li2S']], 'node =', j+1, '\n\n')
"""Calculate change in electrolyte"""
res[offset + ptr['rho_k_el']] = (SV_dot[offset + ptr['rho_k_el']] -
(R_net + (N_io_m - N_io_p)*cat.dyInv)/eps_el
+ SV[offset + ptr['rho_k_el']]*(- SV_dot[offset + ptr['eps_S8']]
- SV_dot[offset + ptr['eps_Li2S']])/eps_el)
# print(R_net, i_ext, eps_S8, '\n\n')
"""Calculate change in delta-phi double layer"""
res[offset + ptr['phi_dl']] = (SV_dot[offset + ptr['phi_dl']] -
(-i_Far + i_el_m - i_el_p)*cat.dyInv/cat.C_dl/A_C)
# print(i_Far, i_ext, eps_S8, '\n\n')
"""Algebraic expression for charge neutrality in all phases"""
res[offset + ptr['phi_ed']] = i_el_m - i_el_p + i_io_m - i_io_p
# print(i_el_m, i_el_p, i_io_m, i_io_p, res[offset+ptr['phi_ed']], '\n')
"""Calculate change in S8 nucleation sites"""
res[offset + ptr['np_S8']] = SV[offset + ptr['np_S8']] - np_S
"""Calculate change in Li2S nucleation sites"""
res[offset + ptr['np_Li2S']] = SV[offset + ptr['np_Li2S']] - np_L
# print(np_L, SV[offset + ptr['np_Li2S']], i_ext, offset, '\n')
"""============================SEPARATOR============================"""
"""INTERIOR NODES"""
for j in np.arange(1, sep.npoints):
i_io_m = i_io_p
N_io_m = N_io_p
s1 = dict(s2)
# Set NEXT node conditions
offset = sep.offsets[int(j)]
s2 = set_state_sep(SV, offset, sep.ptr)
# Shift back to THIS node
offset = sep.offsets[int(j-1)]
D_el = sep.D_el
# THIS node plus face boundary conditions
N_io_p, i_io_p = dst(s1, s2, D_el, sep.dy, sep.dy)
res[offset + sep.ptr['rho_k_el']] = (SV_dot[offset + sep.ptr['rho_k_el']]
- (N_io_m - N_io_p)*sep.dyInv/sep.epsilon_el)
res[offset + sep.ptr['phi']] = i_io_m - i_io_p
"""============================SEPARATOR============================"""
"""ANODE BOUNDARY"""
i_io_m = i_io_p
N_io_m = N_io_p
s1 = dict(s2)
# Shift forward to NEXT node
j = 0; offset = an.offsets[int(j)]
s2 = set_state(SV, offset, an.ptr)
# Shift back to THIS node
offset = sep.offsets[-1]
D_el = sep.D_el
# Current node plus face boundary conditions
N_io_p, i_io_p = dst(s1, s2, D_el, sep.dy, an.dy)
res[offset + sep.ptr['rho_k_el']] = (SV_dot[offset + sep.ptr['rho_k_el']]
- (N_io_m - N_io_p)*sep.dyInv/sep.epsilon_el)
res[offset + sep.ptr['phi']] = i_io_m - i_io_p
"""==============================ANODE=============================="""
"""INTERIOR NODES"""
i_io_m = i_io_p
N_io_m = N_io_p
i_el_m = 0
s1 = dict(s2)
j = 0
offset = an.offsets[int(j)]
i_el_p = i_ext
i_io_p = 0
N_io_p = 0
elyte.X = s1['X_k']
elyte.electric_potential = s1['phi_el']
lithium.electric_potential = s1['phi_ed']
conductor.electric_potential = s1['phi_ed']
sdot_Li = lithium_s.get_net_production_rates(elyte)
sdot_Far = lithium_s.get_net_production_rates(conductor)
R_net = sdot_Li*an.A_Li
i_Far = sdot_Far*an.A_Li*F*an.dy
res[offset + an.ptr['rho_k_el']] = (SV_dot[offset + an.ptr['rho_k_el']]
- (R_net + (N_io_m - N_io_p)*an.dyInv)/an.eps_el)
res[offset + an.ptr['phi_dl']] = (SV_dot[offset + an.ptr['phi_dl']]
- (-i_Far + i_el_m - i_el_p)*an.dyInv/an.C_dl/an.A_Li)
res[offset + an.ptr['phi_ed']] = SV[offset + an.ptr['phi_ed']]
"""==============================ANODE=============================="""
"""CC BOUNDARY"""
# print(i_ext, t, '\n\n')
return res
"========================================================================="
def state_events(self, t, y, yd, sw):
event1 = np.zeros([cat.npoints])
event2 = np.zeros([cat.npoints])
event1 = 1 - y[cat.ptr_vec['eps_S8']]
# event2 = y[cat.ptr_vec['eps_S8']]
event3 = np.zeros([cat.npoints])
event4 = np.zeros([cat.npoints])
event3 = 1 - y[cat.ptr_vec['eps_Li2S']]
# event4 = y[cat.ptr_vec['eps_Li2S']]
event5 = np.zeros([cat.npoints])
# event5 = 2.8 - y[cat.ptr_vec['phi_ed']]
event6 = np.zeros([cat.npoints])
event6 = y[cat.ptr_vec['phi_ed']] - 1.8
event7 = np.zeros([cat.npoints*elyte.n_species])
event7 = y[cat.ptr_vec['rho_k_el']]
events = np.concatenate((event1, event2, event3, event4, event5, event6,
event7))
return events
"========================================================================="
def handle_event(self, solver, event_info):
state_info = event_info[0]
offset = cat.npoints
if any(state_info[0:offset]):
print('Sulfur volume fraction over 1')
raise TerminateSimulation
if any(state_info[1*offset:2*offset]):
print('WARNING: Sulfur volume fraction below 0')
raise TerminateSimulation
elif any(state_info[2*offset:3*offset]):
print('Li2S volume fraction over 1')
raise TerminateSimulation
elif any(state_info[3*offset:4*offset]):
print('Li2S volume fraction below 0')
raise TerminateSimulation
elif any(state_info[4*offset:5*offset]):
print('Cell voltage hit 2.8')
raise TerminateSimulation
elif any(state_info[5*offset:6*offset]):
print('Cell voltage hit 1.8')
raise TerminateSimulation
elif any(state_info):
print('Stop condition')
raise TerminateSimulation
if __name__ == "__main__":
SV_eq, SV_dch, SV_ch, tags = main()
"============================================================================="