Skip to content

Commit

Permalink
dynamics base model
Browse files Browse the repository at this point in the history
  • Loading branch information
@Ombrini
Ombrini committed Oct 24, 2023
1 parent a595a5b commit 6a805c0
Show file tree
Hide file tree
Showing 7 changed files with 155 additions and 57 deletions.
52 changes: 46 additions & 6 deletions mpet/config/configuration.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -293,7 +293,7 @@ def read(self, folder=None, filenamebase='input_dict', full=False):
# only update generated distributions
if section == 'system':
for key in ['psd_num', 'psd_len', 'psd_area', 'psd_vol',
'psd_vol_FracVol', 'G']:
'psd_vol_FracVol', 'G','G_2']:
self[key] = d[key]
elif section in ['anode', 'cathode']:
trode = section[0]
Expand Down Expand Up @@ -444,6 +444,7 @@ def _process_config(self, prevDir=None):
# set default 1C_current_density
limtrode = self['limtrode']
theoretical_1C_current = self[limtrode, 'cap'] / 3600. # A/m^2
# print('Theoretical 1C current density: ', theoretical_1C_current)
param = '1C_current_density'
if self[param] is None:
# set to theoretical value
Expand All @@ -469,6 +470,8 @@ def _process_config(self, prevDir=None):
self._distr_part()
# Gibss free energy, must be done after distr_part
self._G()
# particle distributions for second particle type
self._G_2()
# Electrode parameters that depend on invidividual particle
self._indvPart()

Expand Down Expand Up @@ -519,6 +522,10 @@ def _scale_electrode_parameters(self):

if self[trode, 'lambda'] is not None:
self[trode, 'lambda'] = self[trode, 'lambda'] / kT
if self[trode, 'lambda_1'] is not None:
self[trode, 'lambda_1'] = self[trode, 'lambda_1'] / kT
if self[trode, 'lambda_2'] is not None:
self[trode, 'lambda_2'] = self[trode, 'lambda_2'] / kT
if self[trode, 'lambda_m'] is not None:
self[trode, 'lambda_m'] = self[trode, 'lambda_m'] / kT
if self[trode, 'kintra'] is not None:
Expand Down Expand Up @@ -715,15 +722,16 @@ def _distr_part(self):

def _G(self):
"""
Generate Gibbs free energy distribution and store in config.
Generate Intra-particle connectivity distribution and store in config.
"""
self['G'] = {}
for trode in self['trodes']:
Nvol = self['Nvol'][trode]
Npart = self['Npart'][trode]
mean = self['G_mean'][trode]
stddev = self['G_stddev'][trode]
if np.allclose(stddev, 0, atol=1e-12):
stddev = 50*mean
# stddev = self['G_stddev'][trode]
if np.allclose(stddev, 0, atol=1e-17):
G = mean * np.ones((Nvol, Npart))
else:
var = stddev**2
Expand All @@ -735,6 +743,29 @@ def _G(self):
self['G'][trode] = G * constants.k * constants.T_ref * self['t_ref'] \
/ (constants.e * constants.F * self[trode, 'csmax'] * self['psd_vol'][trode])

def _G_2(self):
"""
Generate Intra-particle connectivity of second specie and store in config.
"""
self['G_2'] = {}
for trode in self['trodes']:
Nvol = self['Nvol'][trode]
Npart = self['Npart'][trode]
mean = self['G_mean_2'][trode]
stddev = 50*mean
# stddev = self['G_stddev_2'][trode]
if np.allclose(stddev, 0, atol=1e-17):
G_2 = mean * np.ones((Nvol, Npart))
else:
var = stddev**2
mu = np.log((mean**2) / np.sqrt(var + mean**2))
sigma = np.sqrt(np.log(var / (mean**2) + 1))
G_2 = np.random.lognormal(mu, sigma, size=(Nvol, Npart))

# scale and store
self['G_2'][trode] = G_2 * constants.k * constants.T_ref * self['t_ref'] \
/ (constants.e * constants.F * self[trode, 'csmax'] * self['psd_vol'][trode])

def _indvPart(self):
"""
Generate particle-specific parameter values and store in config.
Expand Down Expand Up @@ -794,8 +825,17 @@ def _indvPart(self):
* self['t_ref'] / plen**2
self[trode, 'indvPart']['E_D'][i, j] = self[trode, 'E_D'] \
/ (constants.k * constants.N_A * constants.T_ref)
self[trode, 'indvPart']['k0'][i, j] = self[trode, 'k0'] \
/ (constants.e * F_s_ref)
if self[trode, 'k0'] is not None:
self[trode, 'indvPart']['k0'][i, j] = self[trode, 'k0'] \
/ (constants.e * F_s_ref)
if self[trode, 'k0_1'] is not None:
F_s_1ref = plen * constants.N_A * self[trode, 'csmax'] * self[trode,"stoich_1"] / self['t_ref']
self[trode, 'indvPart']['k0_1'][i, j] = self[trode, 'k0_1'] \
/ (constants.e * F_s_1ref)
if self[trode, 'k0_2'] is not None:
F_s_2ref = plen * constants.N_A * self[trode, 'csmax'] * (1 - self[trode,"stoich_1"]) / self['t_ref']
self[trode, 'indvPart']['k0_2'][i, j] = self[trode, 'k0_2'] \
/ (constants.e * F_s_2ref)
self[trode, 'indvPart']['E_A'][i, j] = self[trode, 'E_A'] \
/ (constants.k * constants.N_A * constants.T_ref)
self[trode, 'indvPart']['Rfilm'][i, j] = self[trode, 'Rfilm'] \
Expand Down
5 changes: 3 additions & 2 deletions mpet/config/constants.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,13 @@

#: parameter that are defined per electrode with a ``_{electrode}`` suffix
PARAMS_PER_TRODE = ['Nvol', 'Npart', 'mean', 'stddev', 'cs0', 'simBulkCond', 'sigma_s',
'simPartCond', 'G_mean', 'G_stddev', 'L', 'P_L', 'poros', 'BruggExp',
'simPartCond', 'G_mean', 'G_stddev','simPartCond_2', 'G_mean_2', 'G_stddev_2',
'L', 'P_L', 'poros', 'BruggExp',
'specified_psd', 'rhom', 'cp', 'k_h']
#: subset of ``PARAMS_PER_TRODE``` that is defined for the separator as well
PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp', 'k_h', 'cp', 'rhom']
#: parameters that are defined for each particle, and their type
PARAMS_PARTICLE = {'N': int, 'kappa': float, 'kappa1': float, 'kappa2': float, 'beta_s': float,
'D': float, 'D1': float, 'D2': float, 'k0': float,
'D': float, 'D1': float, 'D2': float, 'k0': float, 'k0_1': float, 'k0_2': float,
'Rfilm': float, 'delta_L': float, 'Omega_a': float,
'E_D': float,'E_A': float}
14 changes: 12 additions & 2 deletions mpet/config/schemas.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,13 @@ def tobool(value):
'G_mean_c': Use(float),
'G_stddev_c': Use(float),
'G_mean_a': Use(float),
'G_stddev_a': Use(float)},
'G_stddev_a': Use(float),
Optional('simPartCond_2_c', default=False): Use(tobool),
Optional('simPartCond_2_a', default=False): Use(tobool),
Optional('G_mean_2_c', default=0): Use(float),
Optional('G_stddev_2_c', default=0): Use(float),
Optional('G_mean_2_a', default=0): Use(float),
Optional('G_stddev_2_a', default=0): Use(float)},
'Geometry': {'L_c': Use(float),
'L_a': Use(float),
'L_s': Use(float),
Expand Down Expand Up @@ -175,10 +181,14 @@ def tobool(value):
Optional('cwet', default=None): Use(float)},
'Reactions': {Optional('rxnType_filename', default=None): str,
'rxnType': str,
'k0': Use(float),
Optional('k0', default=None): Use(float),
Optional('k0_1', default=None): Use(float),
Optional('k0_2', default=None): Use(float),
Optional('E_A', default=0.): Use(float),
'alpha': Use(float),
Optional('lambda', default=None): Use(float),
Optional('lambda_1', default=None): Use(float),
Optional('lambda_2', default=None): Use(float),
Optional('lambda_m', default=None): Use(float),
Optional('kintra', default=None): Use(float),
'Rfilm': Use(float)}}
Expand Down
36 changes: 34 additions & 2 deletions mpet/mod_cell.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,31 +290,63 @@ def DeclareEquations(self):
# Simulate the potential drop along the connected
# particles
simPartCond = config['simPartCond'][trode]
simPartCond_2 = config['simPartCond_2'][trode]
if simPartCond_2:
# list the content of config
stoich_1 = config['c',"stoich_1"]
stoich_2 = 1 - stoich_1
tail_mid = 1e-3
for vInd in range(Nvol[trode]):
phi_bulk = self.phi_bulk[trode](vInd)
for pInd in range(Npart[trode]):
G_l = config["G"][trode][vInd,pInd]

# cbar_l = self.particles[trode][vInd,pInd].cbar()
if simPartCond_2:
c2bar_l = self.particles[trode][vInd,pInd].c2bar()
G_l_2 = config["G_2"][trode][vInd,pInd]

# G_tot_l = G_l + (G_l_2-G_l)/(1+np.exp(300*(stoich_1-cbar_l)))
# G_tot_l = G_l + (G_l_2-G_l)*0.5*(np.tanh((cbar_l - stoich_1)/tail_mid) + 1)
G_tot_l = G_l + (G_l_2-G_l)/(1+np.exp(1000*(0-c2bar_l)))

phi_n = self.phi_part[trode](vInd, pInd)
if pInd == 0: # reference bulk phi
phi_l = phi_bulk
else:
phi_l = self.phi_part[trode](vInd, pInd-1)
if pInd == (Npart[trode] - 1): # No particle at end of "chain"
G_r = 0
G_tot_r = 0
phi_r = phi_n
else:
G_r = config["G"][trode][vInd,pInd+1]

if simPartCond_2:
c2bar_r = self.particles[trode][vInd,pInd+1].c2bar()
# cbar_r = self.particles[trode][vInd,pInd+1].cbar()
G_r_2 = config["G_2"][trode][vInd,pInd+1]
# G_tot_r = G_r + (G_r_2-G_r)/(1+np.exp(300*(stoich_1-cbar_r)))
# G_tot_r = G_r + (G_r_2-G_r)*0.5*(np.tanh((cbar_r - stoich_1)/tail_mid) + 1)
G_tot_r = G_r + (G_r_2-G_r)/(1+np.exp(1000*(0-c2bar_r)))

phi_r = self.phi_part[trode](vInd, pInd+1)
# charge conservation equation around this particle
eq = self.CreateEquation(
"phi_ac_trode{trode}vol{vInd}part{pInd}".format(
vInd=vInd, trode=trode, pInd=pInd))
if simPartCond_2:
G_l_tot = G_tot_l
G_r_tot = G_tot_r
else:
G_l_tot = G_l
G_r_tot = G_r
if simPartCond:
# -dcsbar/dt = I_l - I_r
eq.Residual = (
self.particles[trode][vInd,pInd].dcbardt()
+ ((-G_l * (phi_n - phi_l))
- (-G_r * (phi_r - phi_n))))
+ ((-G_l_tot * (phi_n - phi_l))
- (-G_r_tot * (phi_r - phi_n))))
else:
eq.Residual = self.phi_part[trode](vInd, pInd) - phi_bulk

Expand Down
80 changes: 45 additions & 35 deletions mpet/mod_electrodes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,9 +43,6 @@ def __init__(self, config, trode, vInd, pInd,
self.c2 = dae.daeVariable(
"c2", mole_frac_t, self,
"Concentration in 'layer' 2 of active particle", [self.Dmn])
self.interLayerRxn = dae.daeVariable(
"interLayerRxn", dae.no_t, self,
"Rate of reaction between layers", [self.Dmn])
self.cbar = dae.daeVariable(
"cbar", mole_frac_t, self,
"Average concentration in active particle")
Expand All @@ -63,6 +60,9 @@ def __init__(self, config, trode, vInd, pInd,
if self.get_trode_param("type") not in ["ACR2"]:
self.Rxn1 = dae.daeVariable("Rxn1", dae.no_t, self, "Rate of reaction 1")
self.Rxn2 = dae.daeVariable("Rxn2", dae.no_t, self, "Rate of reaction 2")
self.interLayerRxn = dae.daeVariable("interLayerRxn", dae.no_t, self,
"Rate of reaction between layers",
[self.Dmn])
else:
self.Rxn1 = dae.daeVariable("Rxn1", dae.no_t, self, "Rate of reaction 1", [self.Dmn])
self.Rxn2 = dae.daeVariable("Rxn2", dae.no_t, self, "Rate of reaction 2", [self.Dmn])
Expand Down Expand Up @@ -253,14 +253,24 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises):
else:
eta1_eff = eta1 + self.Rxn1()*self.get_trode_param("Rfilm")
eta2_eff = eta2 + self.Rxn2()*self.get_trode_param("Rfilm")
Rxn1 = self.calc_rxn_rate(
eta1_eff, c1_surf, self.c_lyte(), self.get_trode_param("k0"),
self.get_trode_param("E_A"), self.T_lyte(), act1R_surf, act_lyte,
self.get_trode_param("lambda"), self.get_trode_param("alpha"))
Rxn2 = self.calc_rxn_rate(
eta2_eff, c2_surf, self.c_lyte(), self.get_trode_param("k0"),
self.get_trode_param("E_A"), self.T_lyte(), act2R_surf, act_lyte,
self.get_trode_param("lambda"), self.get_trode_param("alpha"))
if self.get_trode_param("lambda_1") is not None:
Rxn1 = self.calc_rxn_rate(
eta1_eff, c1_surf, self.c_lyte(), self.get_trode_param("k0_1"),
self.get_trode_param("E_A"), self.T_lyte(), act1R_surf, act_lyte,
self.get_trode_param("lambda_1"), self.get_trode_param("alpha"))
Rxn2 = self.calc_rxn_rate(
eta2_eff, c2_surf, self.c_lyte(), self.get_trode_param("k0_2"),
self.get_trode_param("E_A"), self.T_lyte(), act2R_surf, act_lyte,
self.get_trode_param("lambda_2"), self.get_trode_param("alpha"))
else:
Rxn1 = self.calc_rxn_rate(
eta1_eff, c1_surf, self.c_lyte(), self.get_trode_param("k0"),
self.get_trode_param("E_A"), self.T_lyte(), act1R_surf, act_lyte,
self.get_trode_param("lambda"), self.get_trode_param("alpha"))
Rxn2 = self.calc_rxn_rate(
eta2_eff, c2_surf, self.c_lyte(), self.get_trode_param("k0"),
self.get_trode_param("E_A"), self.T_lyte(), act2R_surf, act_lyte,
self.get_trode_param("lambda"), self.get_trode_param("alpha"))
if self.get_trode_param("type") in ["ACR2"]:
for i in range(N):
eq1 = self.CreateEquation("Rxn1_{i}".format(i=i))
Expand Down Expand Up @@ -321,30 +331,30 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises):
c1, c2, mu1R, mu2R, self.get_trode_param("D"), Dfunc,
self.get_trode_param("E_D"), Flux1_bc, Flux2_bc, dr, self.T_lyte(),
noise1, noise2)
if self.get_trode_param("shape") == "sphere":
area_vec = 4*np.pi*edges**2
elif self.get_trode_param("shape") == "cylinder":
area_vec = 2*np.pi*edges # per unit height
RHS1 += -np.diff(Flux1_vec * area_vec)
RHS2 += -np.diff(Flux2_vec * area_vec)
# lambda_m = self.get_trode_param("lambda_m")
# gamma_ts = 1/((1-c1)*(1-c2))
# kinterlayer = self.get_trode_param("kintra")
# delta_mu = mu2R - mu1R - np.log(c2/c1)
# interLayerRxn = kinterlayer/gamma_ts\
# * (c1*np.exp(-(((lambda_m+delta_mu)**2)/(4*lambda_m)))
# - c2*np.exp(-(((lambda_m-delta_mu)**2)/(4*lambda_m))))
# # interLayerRxn = kinterlayer/gamma_ts\
# # * (- c2*np.exp(-(((lambda_m-delta_mu)**2)/(4*lambda_m))))
# # interLayerRxn = (kinterlayer * (1 - c1) * (1 - c2) * (act1R - act2R))
# RxnTerm1 = -interLayerRxn
# RxnTerm2 = interLayerRxn
# for i in range(N):
# eq = self.CreateEquation("interLayerRxn{i}".format(i=i))
# eq.Residual = self.interLayerRxn(i) - interLayerRxn[i]

# RHS1 += RxnTerm1
# RHS2 += RxnTerm2
if self.get_trode_param("shape") == "sphere":
area_vec = 4*np.pi*edges**2
elif self.get_trode_param("shape") == "cylinder":
area_vec = 2*np.pi*edges # per unit height
RHS1 += -np.diff(Flux1_vec * area_vec)
RHS2 += -np.diff(Flux2_vec * area_vec)
# lambda_m = self.get_trode_param("lambda_m")
# gamma_ts = 1/((1-c1)*(1-c2))
# kinterlayer = self.get_trode_param("kintra")
# delta_mu = mu2R - mu1R - np.log(c2/c1)
# interLayerRxn = kinterlayer/gamma_ts\
# * (c1*np.exp(-(((lambda_m+delta_mu)**2)/(4*lambda_m)))
# - c2*np.exp(-(((lambda_m-delta_mu)**2)/(4*lambda_m))))
# # interLayerRxn = kinterlayer/gamma_ts\
# # * (- c2*np.exp(-(((lambda_m-delta_mu)**2)/(4*lambda_m))))
# # interLayerRxn = (kinterlayer * (1 - c1) * (1 - c2) * (act1R - act2R))
# RxnTerm1 = -interLayerRxn
# RxnTerm2 = interLayerRxn
# for i in range(N):
# eq = self.CreateEquation("interLayerRxn{i}".format(i=i))
# eq.Residual = self.interLayerRxn(i) - interLayerRxn[i]

# RHS1 += RxnTerm1
# RHS2 += RxnTerm2

dc1dt_vec = np.empty(N, dtype=object)
dc2dt_vec = np.empty(N, dtype=object)
Expand Down
6 changes: 3 additions & 3 deletions mpet/plot/plot_data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -699,8 +699,8 @@ def animate(tind):
# Set up colors.
# Define if you want smooth or discrete color changes
# Option: "smooth" or "discrete"
color_changes = "discrete"
# color_changes = "smooth"
# color_changes = "discrete"
color_changes = "smooth"
# Discrete color changes:
if color_changes == "discrete":
# Make a discrete colormap that goes from green to yellow
Expand All @@ -721,7 +721,7 @@ def animate(tind):
# Smooth colormap changes:
if color_changes == "smooth":
# generated with colormap.org
cmaps = np.load("colormaps_custom.npz")
cmaps = np.load(r"C:\Users\pierfrancescoo\Documents\Phase-field\mpet-LFMP\mpet\mpet\plot\colormaps_custom.npz")
cmap_data = cmaps["GnYlRd_3"]
cmap = mpl.colors.ListedColormap(cmap_data/255.)

Expand Down
Loading

0 comments on commit 6a805c0

Please sign in to comment.