Coverage for fxdgm/Helper_DoubleLayerCapacity.py: 98%
87 statements
« prev ^ index » next coverage.py v7.6.12, created at 2025-02-13 13:26 +0000
1# Jan Habscheid
2# Jan.Habscheid@rwth-aachen.de
4import numpy as np
5from scipy.optimize import fixed_point, fsolve
8# Helper functions
9def Phi_pot_center(Phi_pot:np.ndarray) -> np.ndarray:
10 '''
11 Returns vector with the center of the electric potential values.
13 Parameters
14 ----------
15 Phi_pot : np.ndarray
16 Input vector with electric potential values.
18 Returns
19 -------
20 np.ndarray
21 Vector with the center of the electric potential values.
22 '''
23 return (Phi_pot[1:] + Phi_pot[:-1]) / 2
25def dx(Phi_pot:np.ndarray) -> float:
26 '''
27 Returns the difference between the first two electric potential values.
29 Assumes that the electric potential values are equally spaced.
31 Parameters
32 ----------
33 Phi_pot : np.ndarray
34 Input vector with electric potential values.
36 Returns
37 -------
38 float
39 Difference between the first two electric potential values.
40 '''
41 return Phi_pot[1] - Phi_pot[0]
43def C_dl(Q_DL:np.ndarray, Phi_pot:np.ndarray) -> np.ndarray:
44 '''
45 Double Layer Capacity
47 Parameters
48 ----------
49 Q_DL : np.ndarray
50 Charge of the system
51 Phi_pot : np.ndarray
52 Electric potential values
54 Returns
55 -------
56 np.ndarray
57 Double Layer Capacity
58 '''
59 return (Q_DL[1:] - Q_DL[:-1]) / dx(Phi_pot)
61def n(p:np.ndarray, K:str|float) -> np.ndarray:
62 '''
63 Calculates the total number density
65 Parameters
66 ----------
67 p : np.ndarray
68 Pressure
69 K : str | float
70 Bulk modulus, use 'incompressible' for an incompressible mixture
72 Returns
73 -------
74 np.ndarray
75 Total number density
76 '''
77 if K == 'incompressible':
78 return np.ones_like(p)
79 n_Dimensionless = (p-1) / K + 1
80 return n_Dimensionless
82def Q_num_(y_A:np.ndarray, y_C:np.ndarray, n:np.ndarray, x:np.ndarray, z_A:float=-1.0, z_C:float=1.0) -> float:
83 '''
84 Calculates the charge of the system
86 Q = ∫_Ω n^F dΩ
88 Parameters
89 ----------
90 y_A : np.ndarray
91 Anion fraction
92 y_C : np.ndarray
93 Cation fraction
94 n : np.ndarray
95 Total number density
96 x : np.ndarray
97 Spatial discretization
98 z_A : float, optional
99 Charge number of anions, by default -1.0
100 z_C : float, optional
101 Charge number of cations, by default 1.0
103 Returns
104 -------
105 float
106 Charge of the system
107 '''
108 nF_dimensionless = (z_C * y_C + z_A * y_A) * n
109 nF_int = -np.trapz(nF_dimensionless, x)
110 return nF_int
112def Q_num_dim(y_A:np.ndarray, y_C:np.ndarray, n:np.ndarray, x:np.ndarray, z_A:float, z_C:float, nR_m:float, e0:float, LR:float) -> float:
113 '''
114 Calculates the charge of the system
116 Q = ∫_Ω n^F dΩ
118 Parameters
119 ----------
120 y_A : np.ndarray
121 Anion fraction
122 y_C : np.ndarray
123 Cation fraction
124 n : np.ndarray
125 Total number density
126 x : np.ndarray
127 Spatial discretization
128 z_A : float, optional
129 Charge number of anions, by default -1.0
130 z_C : float, optional
131 Charge number of cations, by default 1.0
132 nR_m : float
133 Reference number density in 1/m^3
134 e0 : float
135 Dielectric constant
136 LR : float
137 Reference length in m
139 Returns
140 -------
141 float
142 Charge of the system
143 '''
144 Q_DL = Q_num_(y_A, y_C, n, x, z_A, z_C)
145 Q_DL *= nR_m * e0 * LR
146 Q_DL *= 1e+6
147 Q_DL *= 1/(1e+4)
148 return Q_DL
150def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, solvation:float) -> float:
151 '''
152 Calculates charge of the system using the analytical method in dimensionless units
154 Q = λ^2(∂ₓ φᴸ − ∂ₓ φᴿ)
156 Parameters
157 ----------
158 y_A_R : float
159 Value of anion fraction at the right boundary
160 y_C_R : float
161 Value of cation fraction at the right boundary
162 y_N_R : float
163 Value of neutral fraction at the right boundary
164 z_A : float
165 Charge number of anions
166 z_C : float
167 Charge number of cations
168 phi_L : float
169 Electric potential at the left boundary [V]
170 phi_R : float
171 Electric potential at the right boundary
172 p_R : float
173 Pressure at the right boundary
174 K : str | float
175 Bulk modulus, use 'incompressible' for an incompressible mixture
176 Lambda2 : float
177 Dimensionless parameter
178 a2 : float
179 Dimensionless parameter
180 solvation : float
181 Solvation number
183 Returns
184 -------
185 float
186 Charge of the system in dimensionless units
187 '''
188 if solvation != 0:
189 raise ValueError('Solvation number must be 0 for the analytical method')
190 z_N = 0
191 E_p = p_R
192 if K == 'incompressible':
193 # Assume E_p = p_R = 0
194 D_A = y_A_R / (np.exp(-a2*p_R-z_A*phi_R))
195 D_C = y_C_R / (np.exp(-a2*p_R-z_C*phi_R))
196 D_N = y_N_R / (np.exp(-a2*p_R))
198 # Analytical method
199 E_p = p_R
200 CLambda_L = np.log(D_A * np.exp(-z_A * phi_L - E_p) + D_C * np.exp(-z_C * phi_L - E_p) + D_N * np.exp(-z_N * phi_L - E_p))
202 Lambda = np.sqrt(Lambda2)
203 Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2) * (np.sqrt(CLambda_L))
204 return Q_DL
205 else:
206 C_A = y_A_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_A*phi_R))
207 C_C = y_C_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_C*phi_R))
208 C_N = y_N_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_N*phi_R))
210 CLambda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N
212 Left = np.sqrt((1/CLambda_L)**(-1/(a2*K)) + 1 - K - E_p)
214 Lambda = np.sqrt(Lambda2)
215 a = np.sqrt(a2)
217 Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * a * np.sqrt(2) * (Left)# - Right)
218 return Q_DL
219 raise ValueError('Invalid input for K')
221def Q_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, solvation:float) -> float:
222 '''
223 Calculates charge of the system using the analytical method in dimensionless units
225 Q = λ^2(∂ₓ φᴸ − ∂ₓ φᴿ)
227 Parameters
228 ----------
229 y_A_R : float
230 Value of anion fraction at the right boundary
231 y_C_R : float
232 Value of cation fraction at the right boundary
233 y_N_R : float
234 Value of neutral fraction at the right boundary
235 z_A : float
236 Charge number of anions
237 z_C : float
238 Charge number of cations
239 phi_L : float
240 Value of electric potential at the left boundary (dimensionless units)
241 phi_R : float
242 Value of electric potential at the right
243 p_R : float
244 Value of pressure at the right boundary
245 K : str | float
246 Bulk modulus, use 'incompressible' for an incompress
247 Lambda2 : float
248 Dimensionless parameter
249 a2 : float
250 Dimensionless parameter
251 nR_m : float
252 Reference number density in 1/m^3
253 e0 : float
254 Dielectric constant
255 LR : float
256 Reference length in m
257 solvation : float
258 Solvation number
260 Returns
261 -------
262 float
263 Charge of the system in µAs/cm³
264 '''
265 Q_DL = Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, phi_L, phi_R, p_R, K, Lambda2, a2, solvation)
266 Q_DL *= nR_m * e0 * LR
267 Q_DL *= 1e+6
268 Q_DL *= 1/(1e+4)
269 return Q_DL
273def C_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, solvation:float) -> float:
274 '''
275 Calculates the double layer capacity of the system using the analytical method in dimensionless units
277 C_dl = ∂Q/∂ φᴸ
279 Parameters
280 ----------
281 y_A_R : float
282 Value of anion fraction at the right boundary
283 y_C_R : float
284 Value of cation fraction at the right boundary
285 y_N_R : float
286 Value of neutral fraction at the right boundary
287 z_A : float
288 Charge number of anions
289 z_C : float
290 Charge number of cations
291 phi_L : float
292 Electric potential at the left boundary [V]
293 phi_R : float
294 Electric potential at the right boundary
295 p_R : float
296 Pressure at the right boundary
297 K : str | float
298 Bulk modulus, use 'incompressible' for an incompressible mixture
299 Lambda2 : float
300 Dimensionless parameter
301 a2 : float
302 Dimensionless parameter
303 solvation : float
304 Solvation number
306 Returns
307 -------
308 float
309 Double layer capacity of the system in dimensionless units
310 '''
311 if solvation != 0:
312 raise ValueError('Solvation number must be 0 for the analytical method')
313 z_N = 0
314 E_p = p_R
315 if K == 'incompressible':
316 # Assume E_p = p_R = 0
317 D_A = y_A_R / (np.exp(-a2*p_R-z_A*phi_R))
318 D_C = y_C_R / (np.exp(-a2*p_R-z_C*phi_R))
319 D_N = y_N_R / (np.exp(-a2*p_R))
321 CLambda_L = np.log(D_A * np.exp(-z_A * phi_L - E_p * a2) + D_C * np.exp(-z_C * phi_L - E_p * a2) + D_N * np.exp(-z_N * phi_L - E_p * a2))
323 xi_CLambda_L = D_A * np.exp(-z_A * phi_L - E_p * a2) + D_C * np.exp(-z_C * phi_L - E_p * a2) + D_N * np.exp(-z_N * phi_L - E_p * a2)
324 dxi_CLambda_L = - D_A * z_A * np.exp(-z_A * phi_L - E_p * a2) - D_C * z_C * np.exp(-z_C * phi_L - E_p * a2) - D_N * z_N * np.exp(-z_N * phi_L - E_p * a2)
326 Lambda = np.sqrt(Lambda2)
328 PreTerm = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2)
329 DerivedTerm = dxi_CLambda_L / (2 * xi_CLambda_L * np.sqrt(np.log(xi_CLambda_L)))
331 return PreTerm * DerivedTerm
332 else:
333 C_A = y_A_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_A*phi_R))
334 C_C = y_C_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_C*phi_R))
335 C_N = y_N_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_N*phi_R))
337 z = 1 - K - E_p
338 y = 1/(a2*K)
340 CLambda = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
341 dClambda = -C_A * z_A * np.exp(-z_A*phi_L) - C_C * z_C * np.exp(-z_C*phi_L) - C_N * z_N * np.exp(-z_N*phi_L)
343 Lambda = np.sqrt(Lambda2)
344 a = np.sqrt(a2)
346 DerivedTerm = y * CLambda**(y-1) * dClambda / (2 * np.sqrt(CLambda**y + z))
348 PreTerm = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * a * np.sqrt(2)
350 return PreTerm * DerivedTerm
351 raise ValueError('Invalid input for K')
355def C_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, k:float, T:float, solvation:float) -> float:
356 '''
357 Calculates the double layer capacity of the system using the analytical method in dimensionless units
359 C_dl = ∂Q/∂ φᴸ
361 Parameters
362 ----------
363 y_A_R : float
364 Value of anion fraction at the right boundary
365 y_C_R : float
366 Value of cation fraction at the right boundary
367 y_N_R : float
368 Value of neutral fraction at the right boundary
369 z_A : float
370 Charge number of anions
371 z_C : float
372 Charge number of cations
373 phi_L : float
374 Value of electric potential at the left boundary (dimensionless units)
375 phi_R : float
376 Value of electric potential at the right
377 p_R : float
378 Value of pressure at the right boundary
379 K : str | float
380 Bulk modulus, use 'incompressible' for an incompress
381 Lambda2 : float
382 Dimensionless parameter
383 a2 : float
384 Dimensionless parameter
385 nR_m : float
386 Reference number density in 1/m^3
387 e0 : float
388 Dielectric constant
389 LR : float
390 Reference length in m
391 k : float
392 Boltzmann constant
393 T : float
394 Temperature
395 solvation : float
396 Solvation number
398 Returns
399 -------
400 float
401 Double layer capacity of the system in µAs/cm²
402 '''
403 C_DL = C_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, phi_L, phi_R, p_R, K, Lambda2, a2, solvation)
404 C_DL *= nR_m * e0 * LR
405 # ! ToDo
406 C_DL *= 1e+6
407 C_DL *= 1/(1e+4)
408 C_DL *= 1 /(k*T/e0)
409 return C_DL