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 | def calculate_dif_eq_int_values(t, sv, control_variables, operating_inputs, parameters):
"""This functions calculates intermediate values for the calculation of the differential equations
Parameters
----------
t : float
Time (s).
sv : dict
Variables calculated by the solver. They correspond to the fuel cell internal states.
sv is a contraction of solver_variables for enhanced readability.
control_variables : dict
Variables controlled by the user.
operating_inputs : dict
Operating inputs of the fuel cell.
parameters : dict
Parameters of the fuel cell model.
Returns
-------
Mext : float
Molar mass of the ambient air outside the stack (kg/mol).
M_H2_N2_in : float
Molar mass of the inlet gas at the anode side (H2/N2 mixture) (kg/mol).
i_n : float
Internal current density (A/m²).
Masm : float
Molar mass of all the gas species in the anode supply manifold (kg/mol).
Maem : float
Molar mass of all the gas species in the anode external manifold (kg/mol).
Mcsm : float
Molar mass of all the gas species in the cathode supply manifold (kg/mol).
Mcem : float
Molar mass of all the gas species in the cathode external manifold (kg/mol).
rho_Cp0 : dict
Volumetric heat capacity of each component in the stack (J.m-3.K-1).
"""
# Extraction of the variables
C_v_acl, C_v_ccl = sv['C_v_acl'], sv['C_v_ccl']
s_acl, s_ccl = sv['s_acl'], sv['s_ccl']
lambda_acl, lambda_mem, lambda_ccl = sv['lambda_acl'], sv['lambda_mem'], sv['lambda_ccl']
C_H2_acl, C_O2_ccl = sv['C_H2_acl'], sv['C_O2_ccl']
T_acl, T_mem, T_ccl = sv['T_acl'], sv['T_mem'], sv['T_ccl']
Pasm, Paem, Pcsm, Pcem = sv.get('Pasm', None), sv.get('Paem', None), sv.get('Pcsm', None), sv.get('Pcem', None)
Phi_asm, Phi_aem = sv.get('Phi_asm', None), sv.get('Phi_aem', None)
Phi_csm, Phi_cem = sv.get('Phi_csm', None), sv.get('Phi_cem', None)
# Extraction of the operating inputs and the parameters
T_des, y_H2_in = operating_inputs['T_des'], operating_inputs['y_H2_in']
Lgc, nb_channel_in_gc, Lm = parameters['Lgc'], parameters['nb_channel_in_gc'], parameters['Lm']
Hmem, Hacl, Hccl = parameters['Hmem'], parameters['Hacl'], parameters['Hccl']
epsilon_gdl, epsilon_cl = parameters['epsilon_gdl'], parameters['epsilon_cl']
epsilon_mpl, kappa_co, epsilon_mc = parameters['epsilon_mpl'], parameters['kappa_co'], parameters['epsilon_mc']
nb_gc, nb_gdl, nb_mpl = parameters['nb_gc'], parameters['nb_gdl'], parameters['nb_mpl']
t_purge, type_auxiliary, type_purge = parameters['t_purge'], parameters['type_auxiliary'], parameters['type_purge']
# Calculation of intermediate values
C_N2_a_mean = (sum(sv[f'C_N2_agc_{i}'] for i in range(1, nb_gc + 1)) / nb_gc)
C_N2_c_mean = (sum(sv[f'C_N2_cgc_{i}'] for i in range(1, nb_gc + 1)) / nb_gc)
# Physical quantities outside the stack
# Molar masses
M = {}
M['ext'] = Phi_ext * Psat(Text) / Pext * M_H2O + \
y_O2_ext * (1 - Phi_ext * Psat(Text) / Pext) * M_O2 + \
(1 - y_O2_ext) * (1 - Phi_ext * Psat(Text) / Pext) * M_N2
M['H2_N2_in'] = y_H2_in * M_H2 + (1 - y_H2_in) * M_N2
# Physical quantities inside the stack
# Pressures
P = {}
for i in range(1, nb_gc + 1):
P[f'agc_{i}'] = (sv[f'C_v_agc_{i}'] + sv[f'C_H2_agc_{i}'] + sv[f'C_N2_agc_{i}']) * R * sv[f'T_agc_{i}']
P[f'cgc_{i}'] = (sv[f'C_v_cgc_{i}'] + sv[f'C_O2_cgc_{i}'] + sv[f'C_N2_cgc_{i}']) * R * sv[f'T_cgc_{i}']
# Total concentration
C_tot = {}
for i in range(1, nb_gc + 1):
C_tot[f'agc_{i}'] = sv[f'C_v_agc_{i}'] + sv[f'C_H2_agc_{i}'] + sv[f'C_N2_agc_{i}']
C_tot[f'cgc_{i}'] = sv[f'C_v_cgc_{i}'] + sv[f'C_O2_cgc_{i}'] + sv[f'C_N2_cgc_{i}']
# Humidities
Phi = {}
for i in range(1, nb_gc + 1):
Phi[f'cgc_{i}'] = sv[f'C_v_cgc_{i}'] / C_v_sat(sv[f'T_cgc_{i}'])
# H2/O2 ratio in the dry anode/cathode gas mixture (H2/N2 or O2/N2) at the GC
y_O2 = {}
for i in range(1, nb_gc + 1):
y_O2[f'cgc_{i}'] = sv[f'C_O2_cgc_{i}'] / (sv[f'C_O2_cgc_{i}'] + sv[f'C_N2_cgc_{i}'])
# Molar masses
for i in range(1, nb_gc + 1):
M[f'agc_{i}'] = sv[f'C_v_agc_{i}'] * R * T_des / P[f'agc_{i}'] * M_H2O + \
sv[f'C_H2_agc_{i}'] * R * T_des / P[f'agc_{i}'] * M_H2 + \
sv[f'C_N2_agc_{i}'] * R * T_des / P[f'agc_{i}'] * M_N2
M[f'cgc_{i}'] = sv[f'C_v_cgc_{i}'] * R * T_des / P[f'cgc_{i}'] * M_H2O + \
sv[f'C_O2_cgc_{i}'] * R * T_des / P[f'cgc_{i}'] * M_O2 + \
sv[f'C_N2_cgc_{i}'] * R * T_des / P[f'cgc_{i}'] * M_N2
# Density of the gas mixture.
rho = {}
for i in range(1, nb_gc + 1):
rho[f'agc_{i}'] = P[f'agc_{i}'] / (R * sv[f'T_agc_{i}']) * M[f'agc_{i}']
for i in range(1, nb_gc + 1):
rho[f'cgc_{i}'] = P[f'cgc_{i}'] / (R * sv[f'T_cgc_{i}']) * M[f'cgc_{i}']
# Vapor ratio over the gas mixture.
x_H2O_v = {}
for i in range(1, nb_gc + 1):
x_H2O_v[f'agc_{i}'] = sv[f'C_v_agc_{i}'] / (sv[f'C_v_agc_{i}'] + sv[f'C_H2_agc_{i}'] + sv[f'C_N2_agc_{i}'])
for i in range(1, nb_gc + 1):
x_H2O_v[f'cgc_{i}'] = sv[f'C_v_cgc_{i}'] / (sv[f'C_v_cgc_{i}'] + sv[f'C_O2_cgc_{i}'] + sv[f'C_N2_cgc_{i}'])
# Dynamic viscosity of the gas mixture.
mu_gaz = {}
for i in range(1, nb_gc + 1):
mu_gaz[f'agc_{i}'] = mu_mixture_gases(['H2O_v', 'H2'], [x_H2O_v[f'agc_{i}'], 1 - x_H2O_v[f'agc_{i}']],
sv[f'T_agc_{i}'])
for i in range(1, nb_gc + 1):
mu_gaz[f'cgc_{i}'] = mu_mixture_gases(['H2O_v', 'O2', 'N2'],
[x_H2O_v[f'cgc_{i}'], y_O2[f'cgc_{i}'] * (1 - x_H2O_v[f'cgc_{i}']),
(1 - y_O2[f'cgc_{i}']) * (1 - x_H2O_v[f'cgc_{i}'])],
sv[f'T_cgc_{i}'])
# Volumetric heat capacity (J.m-3.K-1)
rho_Cp0 = {
**{f'agdl_{i}': calculate_rho_Cp0('agdl', sv[f'T_agdl_{i}'], C_v=sv[f'C_v_agdl_{i}'],
s=sv[f's_agdl_{i}'], C_H2=sv[f'C_H2_agdl_{i}'], C_N2=C_N2_a_mean, epsilon=epsilon_gdl)
for i in range(1, nb_gdl + 1)},
**{f'ampl_{i}': calculate_rho_Cp0('ampl', sv[f'T_ampl_{i}'], C_v=sv[f'C_v_ampl_{i}'],
s=sv[f's_ampl_{i}'], C_H2=sv[f'C_H2_ampl_{i}'], C_N2=C_N2_a_mean, epsilon=epsilon_mpl)
for i in range(1, nb_mpl + 1)},
'acl': calculate_rho_Cp0('acl', T_acl, C_v=C_v_acl, s=s_acl, lambdaa=lambda_acl, C_N2=C_N2_a_mean, C_H2=C_H2_acl,
epsilon=epsilon_cl, epsilon_mc=epsilon_mc),
'mem': calculate_rho_Cp0('mem', T_mem, lambdaa=lambda_mem),
'ccl': calculate_rho_Cp0('ccl', T_ccl, C_v=C_v_ccl, s=s_ccl, lambdaa=lambda_ccl, C_O2=C_O2_ccl, C_N2=C_N2_c_mean,
epsilon=epsilon_cl, epsilon_mc=epsilon_mc),
**{f'cmpl_{i}': calculate_rho_Cp0('cmpl', sv[f'T_cmpl_{i}'], C_v=sv[f'C_v_cmpl_{i}'],
s=sv[f's_cmpl_{i}'], C_O2=sv[f'C_O2_cmpl_{i}'], C_N2=C_N2_c_mean, epsilon=epsilon_mpl)
for i in range(1, nb_mpl + 1)},
**{f'cgdl_{i}': calculate_rho_Cp0('cgdl', sv[f'T_cgdl_{i}'], C_v=sv[f'C_v_cgdl_{i}'],
s=sv[f's_cgdl_{i}'], C_O2=sv[f'C_O2_cgdl_{i}'], C_N2=C_N2_c_mean, epsilon=epsilon_gdl)
for i in range(1, nb_gdl + 1)}
}
# The crossover current density i_n
T_acl_mem_ccl = average([T_acl, T_mem, T_ccl],
weights=[Hacl / (Hacl + Hmem + Hccl), Hmem / (Hacl + Hmem + Hccl),
Hccl / (Hacl + Hmem + Hccl)])
i_H2 = 2 * F * R * T_acl_mem_ccl / Hmem * C_H2_acl * k_H2(lambda_mem, T_mem, kappa_co)
i_O2 = 4 * F * R * T_acl_mem_ccl / Hmem * C_O2_ccl * k_O2(lambda_mem, T_mem, kappa_co)
i_n = i_H2 + i_O2
# Physical quantities inside the auxiliary system
if parameters["type_auxiliary"] == "forced-convective_cathode_with_anodic_recirculation" or \
parameters["type_auxiliary"] == "forced-convective_cathode_with_flow-through_anode":
pass
# # Purge
# if type_purge == "no_purge":
# k_purge = 0
# elif type_purge == "constant_purge":
# k_purge = 1
# elif type_purge == "periodic_purge":
# purge_time, delta_purge = t_purge
# if (t - int(t / (purge_time + delta_purge)) * (purge_time + delta_purge)) <= purge_time:
# k_purge = 1
# else:
# k_purge = 0
# else:
# raise ValueError("The type_purge variable should be correctly referenced.")
#
# # H2/O2 ratio in the dry anode/cathode gas mixture (H2/N2 or O2/N2) at the EM
# y_H2_aem = (Paem - Phi_aem * Psat(T_des) - C_N2_a * R * T_des) / (Paem - Phi_aem * Psat(T_des))
# y_O2_cem = (Pcem - Phi_cem * Psat(T_cgc) - C_N2_c * R * T_cgc) / (Pcem - Phi_cem * Psat(T_cgc))
#
# # Molar masses at the anode side
# if parameters["type_auxiliary"] == "forced-convective_cathode_with_anodic_recirculation":
# M['asm'] = Phi_asm * Psat(T_des) / Pasm * M_H2O + \
# (1 - Phi_asm * Psat(T_des) / Pasm) * M_H2
# M['aem'] = Phi_aem * Psat(T_des) / Paem * M_H2O + \
# (1 - Phi_aem * Psat(T_des) / Paem) * M_H2
# else: #parameters["type_auxiliary"] == "forced-convective_cathode_with_flow-through_anode":
# M['asm'] = Phi_asm * Psat(T_des) / Pasm * M_H2O + \
# y_H2_in * (1 - Phi_asm * Psat(T_des) / Pasm) * M_H2 + \
# (1 - y_H2_in) * (1 - Phi_asm * Psat(T_des) / Pasm) * M_N2
# M['aem'] = Phi_aem * Psat(T_des) / Paem * M_H2O + \
# y_H2_aem * (1 - Phi_aem * Psat(T_des) / Paem) * M_H2 + \
# (1 - y_H2_aem) * (1 - Phi_aem * Psat(T_des) / Paem) * M_N2
# # Molar masses at the cathode side
# M['csm'] = Phi_csm * Psat(T_des) / Pcsm * M_H2O + \
# y_O2_ext * (1 - Phi_csm * Psat(T_des) / Pcsm) * M_O2 + \
# (1 - y_O2_ext) * (1 - Phi_csm * Psat(T_des) / Pcsm) * M_N2
# M['cem'] = Phi_cem * Psat(T_des) / Pcem * M_H2O + \
# y_O2_cem * (1 - Phi_cem * Psat(T_des) / Pcem) * M_O2 + \
# (1 - y_O2_cem) * (1 - Phi_cem * Psat(T_des) / Pcem) * M_N2
#
# # Density/concentration of the gas mixture.
# C_tot_a_in = Pasm_in / (R * T_des)
# rho_asm = Pasm / (R * T_des) * Masm
# rho_agc = P[f'agc_{i}'] / (R * sv[f'T_agc_{i}']) * Magc
# rho_aem = Paem / (R * T_des) * Maem
# if type_auxiliary == "forced-convective_cathode_with_anodic_recirculation":
# rho_asm_in_re = Pasm_in_re / (R * T_des) * Masm_in_re
# rho_aem_out_re = Paem_out_re / (R * T_des) * Maem_out_re
# else:
# rho_asm_in_re, rho_aem_out_re = None, None
# rho_a_ext = Pext / (R * T_des) * Maem_out
# C_tot_a_ext = Pext / (R * T_des) # Boundary condition: at the exit, pressure and temperature are fixed. So, the total concentration is fixed.
# C_tot_c_in = Pcsm_in / (R * T_des)
# rho_csm = Pcsm / (R * T_des) * Mcsm
# rho_cgc = P[f'cgc_{i}'] / (R * sv[f'T_cgc_{i}']) * Mcgc
# rho_cem = Pcem / (R * T_cgc) * Mcem
# rho_c_ext = Pext / (R * T_des) * Mcem_out
# C_tot_c_ext = Pext * Mcem_out / (R * T_des) # Boundary condition: at the exit, pressure and temperature are fixed. So, the total concentration is fixed.
#
# # Vapor ratio over the gas mixture.
# x_H2O_v_asm = Phi_asm * Psat(T_des) / Pasm
# x_H2O_v_agc = C_v_agc / (C_v_agc + C_H2_agc + C_N2_a)
# x_H2O_v_aem = Phi_aem * Psat(T_des) / Paem
# x_H2O_v_a_ext = Phi_a_ext * Psat(T_des) / Pext
# x_H2O_v_csm = Phi_csm * Psat(T_des) / Pcsm
# x_H2O_v_cgc = C_v_cgc / (C_v_cgc + C_O2_cgc + C_N2_c)
# x_H2O_v_cem = Phi_cem * Psat(T_des) / Pcem
# x_H2O_v_c_ext = Phi_c_ext * Psat(T_des) / Pext
#
# # Molar fraction of H2 in the dry gas mixture (H2/N2)
# y_H2_agc = C_H2_agc / (C_H2_agc + C_N2_a)
# y_O2_cgc = C_O2_cgc / (C_O2_cgc + C_N2_c)
#
# # Dynamic viscosity of the gas mixture at the anode side.
# if type_auxiliary == "forced-convective_cathode_with_anodic_recirculation":
# mu_gaz_asm = mu_mixture_gases(['H2O_v', 'H2'], [x_H2O_v_asm, 1 - x_H2O_v_asm], T_des)
# mu_gaz_agc = mu_mixture_gases(['H2O_v', 'H2'], [x_H2O_v_agc, 1 - x_H2O_v_agc], T_agc)
# mu_gaz_aem = mu_mixture_gases(['H2O_v', 'H2'], [x_H2O_v_aem, 1 - x_H2O_v_aem], T_des)
# mu_gaz_a_ext = mu_mixture_gases(['H2O_v', 'H2'], [x_H2O_v_a_ext, 1 - x_H2O_v_a_ext], T_des)
# else: # type_auxiliary == "forced-convective_cathode_with_flow-through_anode"
# mu_gaz_asm = mu_mixture_gases(['H2O_v', 'H2', 'N2'],
# [x_H2O_v_asm, y_H2_in * (1 - x_H2O_v_asm), (1 - y_H2_in) * (1 - x_H2O_v_asm)],
# T_des)
# mu_gaz_agc = mu_mixture_gases(['H2O_v', 'H2', 'N2'],
# [x_H2O_v_agc, y_H2_agc * (1 - x_H2O_v_agc),
# (1 - y_H2_agc) * (1 - x_H2O_v_agc)], T_agc)
# mu_gaz_aem = mu_mixture_gases(['H2O_v', 'H2', 'N2'],
# [x_H2O_v_aem, y_H2_aem * (1 - x_H2O_v_aem),
# (1 - y_H2_aem) * (1 - x_H2O_v_aem)], T_des)
# mu_gaz_a_ext = mu_mixture_gases(['H2O_v', 'H2', 'N2'],
# [x_H2O_v_a_ext, y_H2_aem_out * (1 - x_H2O_v_a_ext), (1 - y_H2_aem_out) * (1 - x_H2O_v_a_ext)],
# T_des)
# # Dynamic viscosity of the gas mixture at the cathode side.
# mu_gaz_csm = mu_mixture_gases(['H2O_v', 'O2', 'N2'],
# [x_H2O_v_csm, y_O2_ext * (1 - x_H2O_v_csm), (1 - y_O2_ext) * (1 - x_H2O_v_csm)],
# T_des)
# mu_gaz_cgc = mu_mixture_gases(['H2O_v', 'O2', 'N2'],
# [x_H2O_v_cgc, y_O2_cgc * (1 - x_H2O_v_cgc), (1 - y_O2_cgc) * (1 - x_H2O_v_cgc)],
# T_cgc)
# mu_gaz_cem = mu_mixture_gases(['H2O_v', 'O2', 'N2'],
# [x_H2O_v_cem, y_O2_cem * (1 - x_H2O_v_cem), (1 - y_O2_cem) * (1 - x_H2O_v_cem)],
# T_des)
# mu_gas_c_ext = mu_mixture_gases(['H2O_v', 'O2', 'N2'],
# [x_H2O_v_c_ext, y_O2_cem_out * (1 - x_H2O_v_c_ext),
# (1 - y_O2_cem_out) * (1 - x_H2O_v_c_ext)],
# T_des)
#
# # Boundary velocities
# if type_auxiliary == "forced-convective_cathode_with_anodic_recirculation":
# v_re = Ware / rho_aem_out_re / A_T_a
# else: # type_auxiliary == "forced-convective_cathode_with_flow-through_anode"
# v_re = None
else: # parameters["type_auxiliary"] == "no_auxiliary"
# Set to None the variables not used when "no_auxiliary" system is considered
v_re, Lman_to_endplate, Lman_to_man_gc, k_purge = [None] * 4
return {'rho_Cp0': rho_Cp0, 'v_re': v_re, 'k_purge': k_purge, 'rho': rho, 'C_tot': C_tot, 'mu_gaz': mu_gaz, 'P': P,
'i_n': i_n}
|