From adcb309f80db55b4fadbb2d0bfe606b6dc49975b Mon Sep 17 00:00:00 2001
From: JanHab <Jan.Habscheid@web.de>
Date: Mon, 19 Aug 2024 11:04:15 +0200
Subject: [PATCH] DLKap - changed notation for n
---
.../DoubleLayerCapacity/Compressibility.py | 178 ++++++++++++++++++
examples/DoubleLayerCapacity/Lambda.py | 2 +-
examples/DoubleLayerCapacity/Molarity.py | 2 +-
examples/DoubleLayerCapacity/Solvation.py | 2 +-
.../ValidationAnalyticalMethod.py | 4 +-
src/Helper_DoubleLayerCapacity.py | 2 +-
...Helper_DoubleLayerCapacity.cpython-312.pyc | Bin 8440 -> 8438 bytes
7 files changed, 184 insertions(+), 6 deletions(-)
create mode 100644 examples/DoubleLayerCapacity/Compressibility.py
diff --git a/examples/DoubleLayerCapacity/Compressibility.py b/examples/DoubleLayerCapacity/Compressibility.py
new file mode 100644
index 0000000..a7fb37a
--- /dev/null
+++ b/examples/DoubleLayerCapacity/Compressibility.py
@@ -0,0 +1,178 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
+sys.path.insert(0, src_path)
+
+from Eq02 import solve_System_2eq
+from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
+
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Import plotting library
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Define Parameter
+e0 = 1.602e-19 # [As]
+k = 1.381e-23 # [J/K]
+T = 293.15 # [K]
+epsilon0 = 8.85e-12 #[F/m]
+F = 9.65e+4 # [As/mol]
+NA = 6.022e+23 # [1/mol] - Avogadro constant
+nR_mol = 55
+nR_m = nR_mol * NA * 1/(1e-3)# [1/m^3]
+pR = 1.01325 * 1e+5 # [Pa]
+LR = 20e-8
+chi = 80 # [-]
+
+# Parameter and bcs for the electrolyte
+a2 = (pR)/(nR_m * k * T)
+K = 'incompressible'
+kappa = 0
+Molarity = 0.01
+y_R = Molarity / nR_mol
+z_A, z_C = -1.0, 1.0
+
+K_vec = ['incompressible', 100_000, 50_000, 20_000]
+Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
+a2 = (pR)/(nR_m * k * T)
+
+# Solver settings
+number_cells = 1024
+p_right = 0
+phi_right = 0
+
+# phi^L domain
+Vol_start = 0
+Volt_end = 0.75
+n_Volts = 10
+
+phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
+
+y_A, y_C, y_S, phi, p, x, n_sol = [], [], [], [], [], [], []
+for index, K_ in enumerate(K_vec):
+ print('-------------------------------------------------------------------')
+ print('K: ', str(K_))
+ print('-------------------------------------------------------------------')
+ phi_inner, y_A_inner, y_C_inner, y_S_inner, p_inner,x_inner, n_inner = [], [], [], [], [], [], []
+ for i, phi_bcs in enumerate(phi_left):
+ print('-------------------------------------------------------------------')
+ print('K: ', str(K_))
+ print('phi_bcs: ', phi_bcs)
+ print('-------------------------------------------------------------------')
+ # ToDo: relax_param
+ relax_param = 0.3
+ if phi_bcs > 20:
+ relax_param = 0.05
+ if phi_bcs > 25:
+ relax_param = 0.03
+ y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K_, Lambda2, a2, number_cells, relax_param=relax_param, refinement_style='hard_hard_log', rtol=1e-4, max_iter=10_000, return_type='Vector')
+ y_A_inner.append(y_A_)
+ y_C_inner.append(y_C_)
+ y_S_inner.append(1 - y_A_ - y_C_)
+ phi_inner.append(phi_*(k*T/e0))
+ p_inner.append(p_)
+ n_inner.append(n(p_, K_))
+ x_inner.append(x_)
+ y_A.append(y_A_inner)
+ y_C.append(y_C_inner)
+ y_S.append(y_S_inner)
+ phi.append(phi_inner)
+ p.append(p_inner)
+ x.append(x_inner)
+ n_sol.append(n_inner)
+phi_left *= (k*T)/e0
+
+
+
+# Charge
+Q = []
+for i in range(len(K_vec)):
+ Q_inner = []
+ for j in range(len(phi_left)):
+ # ToDo: nR_m -> n
+ Q_inner.append(Q_num_(y_A[i][j], y_C[i][j], n_sol[i][j], x[i][j]))
+ Q.append(Q_inner)
+Q = np.array(Q)
+
+# Double layer capacity
+dx_ = phi_left[1] - phi_left[0] # [1/V], uniformly distributed
+C_DL = (Q[:,1:] - Q[:,:-1])/dx_ # [µAs/cm³] # ? Right unit?
+C_DL = np.array(C_DL)
+phi_left_center = (phi_left[1:] + phi_left[:-1])/2 # Center points for plotting C_DL
+
+
+
+# Visualizations
+# Charge
+plt.figure()
+[plt.plot(phi_left, Q[i], label=f'K = {K_vec[i]}') for i in range(len(K_vec))]
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('Q $[-]$')
+plt.tight_layout()
+plt.show()
+
+# Double layer capacity
+plt.figure()
+[plt.plot(phi_left_center, C_DL[i], label=f'K = {K_vec[i]}') for i in range(len(K_vec))]
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('$C_{dl} [-]$')
+plt.tight_layout()
+plt.show()
+
+# y_alpha_L
+y_A_np, y_C_np, y_S_np = np.array(y_A), np.array(y_C), np.array(y_S)
+markers = ['--', '-', ':']
+colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple']
+plt.figure()
+for i in range(len(K_vec)):
+ clr = colors[i]
+ plt.plot(phi_left, y_A_np[i,:,0], '--', color=clr)
+ plt.plot(phi_left, y_C_np[i,:,0], '-', color=clr)
+ plt.plot(phi_left, y_S_np[i,:,0], ':', color=clr)
+dummy, = plt.plot(0, 0, color='grey', linestyle='--', label='$y_A$')
+dummy, = plt.plot(0, 0, color='grey', linestyle='-', label='$y_C$')
+dummy, = plt.plot(0, 0, color='grey', linestyle=':', label='$y_S$')
+for index, K_ in enumerate(K_vec):
+ dummy, = plt.plot(0, 0, color=colors[index], linestyle='-', label=f'K = {K_}')
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('$y_\\alpha^L$')
+plt.tight_layout()
+plt.show()
+
+n_sol_np = np.array(n_sol)
+n_A_np, n_C_np, n_S_np = y_A_np * n_sol_np, y_C_np * n_sol_np, y_S_np * n_sol_np
+plt.figure()
+clr = colors[0]
+for i in range(len(K_vec)):
+ clr = colors[i]
+ plt.plot(phi_left, n_A_np[i,:,0], '--', color=clr)
+ plt.plot(phi_left, n_C_np[i,:,0], '-', color=clr)
+ plt.plot(phi_left, n_S_np[i,:,0], ':', color=clr)
+dummy, = plt.plot(0, 0, color='grey', linestyle='--', label='$n_A$')
+dummy, = plt.plot(0, 0, color='grey', linestyle='-', label='$n_C$')
+dummy, = plt.plot(0, 0, color='grey', linestyle=':', label='$n_S$')
+for index, K_ in enumerate(K_vec):
+ dummy, = plt.plot(0, 0, color=colors[index], linestyle='-', label=f'K = {K_}')
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('$n_\\alpha^L [-]$')
+plt.tight_layout()
+plt.show()
diff --git a/examples/DoubleLayerCapacity/Lambda.py b/examples/DoubleLayerCapacity/Lambda.py
index fb1ad02..fdf0c34 100644
--- a/examples/DoubleLayerCapacity/Lambda.py
+++ b/examples/DoubleLayerCapacity/Lambda.py
@@ -12,7 +12,7 @@ src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'sr
sys.path.insert(0, src_path)
from Eq02 import solve_System_2eq
-from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n_num, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
+from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
# Remove the src directory from sys.path after import
diff --git a/examples/DoubleLayerCapacity/Molarity.py b/examples/DoubleLayerCapacity/Molarity.py
index bcb103d..00e231c 100644
--- a/examples/DoubleLayerCapacity/Molarity.py
+++ b/examples/DoubleLayerCapacity/Molarity.py
@@ -12,7 +12,7 @@ src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'sr
sys.path.insert(0, src_path)
from Eq02 import solve_System_2eq
-from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n_num, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
+from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
# Remove the src directory from sys.path after import
diff --git a/examples/DoubleLayerCapacity/Solvation.py b/examples/DoubleLayerCapacity/Solvation.py
index 65c8c2c..70fe587 100644
--- a/examples/DoubleLayerCapacity/Solvation.py
+++ b/examples/DoubleLayerCapacity/Solvation.py
@@ -12,7 +12,7 @@ src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'sr
sys.path.insert(0, src_path)
from Eq02 import solve_System_2eq
-from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n_num, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
+from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
# Remove the src directory from sys.path after import
diff --git a/examples/DoubleLayerCapacity/ValidationAnalyticalMethod.py b/examples/DoubleLayerCapacity/ValidationAnalyticalMethod.py
index 4655f65..c9d626e 100644
--- a/examples/DoubleLayerCapacity/ValidationAnalyticalMethod.py
+++ b/examples/DoubleLayerCapacity/ValidationAnalyticalMethod.py
@@ -12,7 +12,7 @@ src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'sr
sys.path.insert(0, src_path)
from Eq02 import solve_System_2eq
-from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n_num, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
+from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
# Remove the src directory from sys.path after import
@@ -75,7 +75,7 @@ for i, phi_bcs in enumerate(phi_left):
Q_num = []
for j in range(len(phi_left)):
- Q_num.append(Q_num_(y_A_num[j], y_C_num[j], n_num(p_num[j], K), x_num[j]))
+ Q_num.append(Q_num_(y_A_num[j], y_C_num[j], n(p_num[j], K), x_num[j]))
Q_num = np.array(Q_num)
dx_ = phi_left[1] - phi_left[0] # [1/V], Assumption: phi^L is uniformly distributed
diff --git a/src/Helper_DoubleLayerCapacity.py b/src/Helper_DoubleLayerCapacity.py
index 4064f56..50df4ba 100644
--- a/src/Helper_DoubleLayerCapacity.py
+++ b/src/Helper_DoubleLayerCapacity.py
@@ -59,7 +59,7 @@ def C_dl(Q_DL:np.ndarray, Phi_pot:np.ndarray) -> np.ndarray:
'''
return (Q_DL[1:] - Q_DL[:-1]) / dx(Phi_pot)
-def n_num(p:np.ndarray, K:str|float) -> np.ndarray:
+def n(p:np.ndarray, K:str|float) -> np.ndarray:
'''
Calculates the total number density
diff --git a/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc b/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc
index e6acfa710b8e3b774a12597414fa1b95d243aa11..278a45c3840a132f67cf7d53ce418765cdeb8a48 100644
GIT binary patch
delta 230
zcmez2_|1{`G%qg~0}vE)97@};kynC?i7{`oHkUtZkrYsLb1hdNlK_ZSBnu+sK!p5c
zX<iFPmCYf%c^p!zAQ3eXp^i`hVrfjiC@VAhhiERN_2x7&e|AM{pr|H)5nc@v0yZGM
zwjja|MCeZbCu_jyFj-%2E4K<r0Axy$!{qmJvW#w%dF3mmnfkdpxjrxmGfGcLo}oH1
S6U<!TwA^c<*XD)t;fw%&qcX(+
delta 226
zcmez7_`{L+G%qg~0}!<AKajR%Bd-J(8*5&CUTN-RbuNGATa0;|tGW7^1Z06KisV3q
zJcv-3EXixZsJc0bH;+S94J4utA~cW`+!B>xD$<+$S+s!BW^;;|Kf8<#P*#(_2&Zzv
zB3qDtI}l+HBJ?Kzku_j+oUALim0J}g05YM-aq?R^Sw{EC-13#uPW@h;ULP2Q8Kq|^
W&q$rA31+TfTF<qTYx8{ha7F-$tToI4
--
GitLab