Added DoubleLayerCapacity for incompressible mixture

examples/DoubleLayerCapacity/Lambda.py

+'''
+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_num, 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
+phi_R = 0.0
+p_R = 0
+
+# Solver settings
+number_cells = 1024
+relax_param = 0.03
+p_right = 0
+
+# phi^L domain
+Vol_start = -1.0
+Volt_end = 1.0
+n_Volts = 5_000
+
+Phi_Pot_Diff_dim = np.linspace(Vol_start, Volt_end, n_Volts)
+Phi_Pot_Diff_dimless = Phi_Pot_Diff_dim * e0/(k*T)
+
+
+
+
+# Lambda2
+Lambda2_vec = np.array([1e-6, 1e-7, 1e-8])
+
+Q_DL_dim_ = []
+Q_DL_dimless_ = []
+for Lambda2 in Lambda2_vec:
+    y_A_R, y_C_R = Molarity / nR_mol, Molarity / nR_mol
+    y_N_R = 1 - y_A_R - y_C_R
+    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, kappa))
+    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, nR_m, e0, LR, kappa))
+    
+C_DL_dim = [C_dl(q_dl, Phi_Pot_Diff_dim) for q_dl in Q_DL_dim_]
+C_DL_dimless = [C_dl(q_dl, Phi_Pot_Diff_dimless) for q_dl in Q_DL_dimless_]
+
+
+
+# Plotting
+# Charge
+fig = plt.figure()
+colors = ['tab:blue', 'tab:orange', 'tab:green']
+color_dimless = 'tab:purple'
+color_dim = 'tab:red'
+ax = fig.add_subplot(111, label="1")
+ax2 = fig.add_subplot(111, label="2", frame_on=False)
+
+# Plot dimensional data
+for i, q_dl in enumerate(Q_DL_dim_):
+    ax.plot(Phi_Pot_Diff_dim, q_dl, color=colors[i], label=f'$\lambda^2$: {Lambda2_vec[i]}')
+ax.grid()
+ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
+ax.set_ylabel('$Q$ [\u03bc$F/cm^3]$', color=color_dim)
+ax.tick_params(axis='x', colors=color_dim)
+ax.tick_params(axis='y', colors=color_dim)
+
+for i, q_dl in enumerate(Q_DL_dimless_):
+    ax2.plot(Phi_Pot_Diff_dimless, q_dl, color=colors[i])
+ax2.xaxis.tick_top()
+ax2.yaxis.tick_right()
+ax2.set_xlabel('$\delta \\varphi [-]$', color=color_dimless) 
+ax2.set_ylabel('$Q [-]$', color=color_dimless)       
+ax2.xaxis.set_label_position('top') 
+ax2.yaxis.set_label_position('right') 
+ax2.tick_params(axis='x', colors=color_dimless)
+ax2.tick_params(axis='y', colors=color_dimless)
+
+fig.legend()
+fig.tight_layout()
+fig.show()
+
+
+# Double Layer Capacity
+fig = plt.figure()
+color_dimless = 'tab:purple'
+color_dim = 'tab:red'
+ax = fig.add_subplot(111, label="1")
+ax2 = fig.add_subplot(111, label="2", frame_on=False)
+
+# Plot dimensional data
+for i, c_dl in enumerate(C_DL_dim):
+    ax.plot(Phi_pot_center(Phi_Pot_Diff_dim), c_dl, color=colors[i], label=f'$\lambda^2$: {Lambda2_vec[i]}')
+ax.grid()
+ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
+ax.set_ylabel('$C_{dl}$ [\u03bc$F/cm^2]$', color=color_dim)
+ax.tick_params(axis='x', colors=color_dim)
+ax.tick_params(axis='y', colors=color_dim)
+
+for i, c_dl in enumerate(C_DL_dimless):
+    ax2.plot(Phi_pot_center(Phi_Pot_Diff_dimless), c_dl, color=colors[i])
+# ax2.grid()
+ax2.xaxis.tick_top()
+ax2.yaxis.tick_right()
+ax2.set_xlabel('$\delta \\varphi [-]$', color=color_dimless) 
+ax2.set_ylabel('$Q [-]$', color=color_dimless)       
+ax2.xaxis.set_label_position('top') 
+ax2.yaxis.set_label_position('right') 
+ax2.tick_params(axis='x', colors=color_dimless)
+ax2.tick_params(axis='y', colors=color_dimless)
+
+fig.legend()
+fig.tight_layout()
+fig.show()
\ No newline at end of file

examples/DoubleLayerCapacity/Molarity.py

+'''
+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_num, 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
+Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
+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
+phi_R = 0.0
+p_R = 0
+
+# Solver settings
+number_cells = 1024
+relax_param = 0.03
+p_right = 0
+
+# phi^L domain
+Vol_start = -1.0
+Volt_end = 1.0
+n_Volts = 5_000
+
+Phi_Pot_Diff_dim = np.linspace(Vol_start, Volt_end, n_Volts)
+Phi_Pot_Diff_dimless = Phi_Pot_Diff_dim * e0/(k*T)
+
+
+
+
+# Molarity
+Molarity = np.array([0.01, 0.1, 1])
+
+Q_DL_dim_ = []
+Q_DL_dimless_ = []
+for mol in Molarity:
+    y_A_R, y_C_R = mol / nR_mol, mol / nR_mol
+    y_N_R = 1 - y_A_R - y_C_R
+    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, kappa))
+    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, nR_m, e0, LR, kappa))
+    
+C_DL_dim = [C_dl(q_dl, Phi_Pot_Diff_dim) for q_dl in Q_DL_dim_]
+C_DL_dimless = [C_dl(q_dl, Phi_Pot_Diff_dimless) for q_dl in Q_DL_dimless_]
+
+
+
+# fig, axs = plt.subplots(ncols=2, figsize=(10, 5))
+# ax1 = axs[0]
+# ax2 = axs[1]
+
+# # Plot the first set of data on the primary y-axis
+# colors = ['tab:blue', 'tab:orange', 'tab:green']
+# ax1.set_xlabel('Phi_Pot_Diff')
+# ax1.set_ylabel('$Q$ [\u03bc$F/cm^3]$')#, color=color)
+# for i, q_dl in enumerate(Q_DL_dim_):
+#     ax1.plot(Phi_Pot_Diff_dim, q_dl, color=colors[i], label=f'Molarity: {Molarity[i]}')
+# ax1.tick_params(axis='y')#, labelcolor=color)
+# ax1.grid()
+
+# # Create a secondary y-axis
+# ax1_twin = ax1.twinx()
+# ax1_twin.set_ylabel('$Q [-]$')#, color=color)
+# for i, q_dl in enumerate(Q_DL_dimless_):
+#     ax1_twin.plot(Phi_Pot_Diff_dimless, q_dl, color=colors[i])
+# ax1_twin.tick_params(axis='y')#, labelcolor=color2)
+# ax1_twin.grid()
+
+# ax2.set_xlabel('Phi_Pot_Diff')
+# ax2.set_ylabel('$C_{dl}$ [\u03bc$F/cm^2]$')#, color=color)
+# for i, c_dl in enumerate(C_DL_dim):
+#     ax2.plot(Phi_pot_center(Phi_Pot_Diff_dim), c_dl, color=colors[i])
+# ax2.tick_params(axis='y')#, labelcolor=color)
+# ax2.grid()
+
+# # Create a secondary y-axis
+# ax2_twin = ax2.twinx()
+# ax2_twin.set_ylabel('$C_{dl} [-]$')#, color=color)
+# for i, c_dl in enumerate(C_DL_dimless):
+#     ax2_twin.plot(Phi_pot_center(Phi_Pot_Diff_dimless), c_dl, color=colors[i])
+# ax2_twin.tick_params(axis='y')#, labelcolor=color2)
+# ax2_twin.grid()
+
+
+# fig.tight_layout()  # Adjust layout to prevent overlap
+# fig.legend()
+# fig.show()
+
+# # plt.figure()
+# # [plt.plot(Phi_pot_center(Phi_Pot_Diff_dim), c_dl) for c_dl in C_DL_dim]
+# # plt.show()
+
+# # plt.figure()
+# # [plt.plot(Phi_pot_center(Phi_Pot_Diff_dimless), c_dl) for c_dl in C_DL_dimless]
+# # plt.show()
+
+
+# Charge
+fig = plt.figure()
+colors = ['tab:blue', 'tab:orange', 'tab:green']
+color_dimless = 'tab:purple'
+color_dim = 'tab:red'
+ax = fig.add_subplot(111, label="1")
+ax2 = fig.add_subplot(111, label="2", frame_on=False)
+
+# Plot dimensional data
+for i, q_dl in enumerate(Q_DL_dim_):
+    ax.plot(Phi_Pot_Diff_dim, q_dl, color=colors[i], label=f'M: {Molarity[i]}')
+ax.grid()
+ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
+ax.set_ylabel('$Q$ [\u03bc$F/cm^3]$', color=color_dim)
+ax.tick_params(axis='x', colors=color_dim)
+ax.tick_params(axis='y', colors=color_dim)
+
+for i, q_dl in enumerate(Q_DL_dimless_):
+    ax2.plot(Phi_Pot_Diff_dimless, q_dl, color=colors[i])
+ax2.xaxis.tick_top()
+ax2.yaxis.tick_right()
+ax2.set_xlabel('$\delta \\varphi [-]$', color=color_dimless) 
+ax2.set_ylabel('$Q [-]$', color=color_dimless)       
+ax2.xaxis.set_label_position('top') 
+ax2.yaxis.set_label_position('right') 
+ax2.tick_params(axis='x', colors=color_dimless)
+ax2.tick_params(axis='y', colors=color_dimless)
+
+fig.legend()
+fig.tight_layout()
+fig.show()
+
+
+# Double Layer Capacity
+fig = plt.figure()
+color_dimless = 'tab:purple'
+color_dim = 'tab:red'
+ax = fig.add_subplot(111, label="1")
+ax2 = fig.add_subplot(111, label="2", frame_on=False)
+
+# Plot dimensional data
+for i, c_dl in enumerate(C_DL_dim):
+    ax.plot(Phi_pot_center(Phi_Pot_Diff_dim), c_dl, color=colors[i], label=f'M: {Molarity[i]}')
+ax.grid()
+ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
+ax.set_ylabel('$C_{dl}$ [\u03bc$F/cm^2]$', color=color_dim)
+ax.tick_params(axis='x', colors=color_dim)
+ax.tick_params(axis='y', colors=color_dim)
+
+for i, c_dl in enumerate(C_DL_dimless):
+    ax2.plot(Phi_pot_center(Phi_Pot_Diff_dimless), c_dl, color=colors[i])
+# ax2.grid()
+ax2.xaxis.tick_top()
+ax2.yaxis.tick_right()
+ax2.set_xlabel('$\delta \\varphi [-]$', color=color_dimless) 
+ax2.set_ylabel('$Q [-]$', color=color_dimless)       
+ax2.xaxis.set_label_position('top') 
+ax2.yaxis.set_label_position('right') 
+ax2.tick_params(axis='x', colors=color_dimless)
+ax2.tick_params(axis='y', colors=color_dimless)
+
+fig.legend()
+fig.tight_layout()
+fig.show()
\ No newline at end of file

examples/DoubleLayerCapacity/Solvation.py

+'''
+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_num, 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
+Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
+a2 = (pR)/(nR_m * k * T)
+K = 'incompressible'
+Molarity = 0.01
+y_R = Molarity / nR_mol
+z_A, z_C = -1.0, 1.0
+phi_R = 0.0
+p_R = 0
+
+# Solver settings
+number_cells = 1024
+relax_param = 0.03
+p_right = 0
+
+# phi^L domain
+Vol_start = -1.0
+Volt_end = 1.0
+n_Volts = 5_000
+
+Phi_Pot_Diff_dim = np.linspace(Vol_start, Volt_end, n_Volts)
+Phi_Pot_Diff_dimless = Phi_Pot_Diff_dim * e0/(k*T)
+
+
+
+
+# kappa
+kappa_vec = np.array([0, 5, 10, 15])
+
+Q_DL_dim_ = []
+Q_DL_dimless_ = []
+for kappa in kappa_vec:
+    y_A_R, y_C_R = Molarity / nR_mol, Molarity / nR_mol
+    y_N_R = 1 - y_A_R - y_C_R
+    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, kappa))
+    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, nR_m, e0, LR, kappa))
+    
+C_DL_dim = [C_dl(q_dl, Phi_Pot_Diff_dim) for q_dl in Q_DL_dim_]
+C_DL_dimless = [C_dl(q_dl, Phi_Pot_Diff_dimless) for q_dl in Q_DL_dimless_]
+
+
+
+# Plotting
+# Charge
+fig = plt.figure()
+colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:cyan']
+color_dimless = 'tab:purple'
+color_dim = 'tab:red'
+ax = fig.add_subplot(111, label="1")
+ax2 = fig.add_subplot(111, label="2", frame_on=False)
+
+# Plot dimensional data
+for i, q_dl in enumerate(Q_DL_dim_):
+    ax.plot(Phi_Pot_Diff_dim, q_dl, color=colors[i], label=f'$\kappa$: {kappa_vec[i]}')
+ax.grid()
+ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
+ax.set_ylabel('$Q$ [\u03bc$F/cm^3]$', color=color_dim)
+ax.tick_params(axis='x', colors=color_dim)
+ax.tick_params(axis='y', colors=color_dim)
+
+for i, q_dl in enumerate(Q_DL_dimless_):
+    ax2.plot(Phi_Pot_Diff_dimless, q_dl, color=colors[i])
+ax2.xaxis.tick_top()
+ax2.yaxis.tick_right()
+ax2.set_xlabel('$\delta \\varphi [-]$', color=color_dimless) 
+ax2.set_ylabel('$Q [-]$', color=color_dimless)       
+ax2.xaxis.set_label_position('top') 
+ax2.yaxis.set_label_position('right') 
+ax2.tick_params(axis='x', colors=color_dimless)
+ax2.tick_params(axis='y', colors=color_dimless)
+
+fig.legend()
+fig.tight_layout()
+fig.show()
+
+
+# Double Layer Capacity
+fig = plt.figure()
+color_dimless = 'tab:purple'
+color_dim = 'tab:red'
+ax = fig.add_subplot(111, label="1")
+ax2 = fig.add_subplot(111, label="2", frame_on=False)
+
+# Plot dimensional data
+for i, c_dl in enumerate(C_DL_dim):
+    ax.plot(Phi_pot_center(Phi_Pot_Diff_dim), c_dl, color=colors[i], label=f'$\kappa$: {kappa_vec[i]}')
+ax.grid()
+ax.set_xlabel('$\delta \\varphi [nm]$', color=color_dim)
+ax.set_ylabel('$C_{dl}$ [\u03bc$F/cm^2]$', color=color_dim)
+ax.tick_params(axis='x', colors=color_dim)
+ax.tick_params(axis='y', colors=color_dim)
+
+for i, c_dl in enumerate(C_DL_dimless):
+    ax2.plot(Phi_pot_center(Phi_Pot_Diff_dimless), c_dl, color=colors[i])
+# ax2.grid()
+ax2.xaxis.tick_top()
+ax2.yaxis.tick_right()
+ax2.set_xlabel('$\delta \\varphi [-]$', color=color_dimless) 
+ax2.set_ylabel('$Q [-]$', color=color_dimless)       
+ax2.xaxis.set_label_position('top') 
+ax2.yaxis.set_label_position('right') 
+ax2.tick_params(axis='x', colors=color_dimless)
+ax2.tick_params(axis='y', colors=color_dimless)
+
+fig.legend()
+fig.tight_layout()
+fig.show()
\ No newline at end of file

examples/DoubleLayerCapacity/ValidationAnalyticalMethod.py

@@ -11,12 +11,99 @@ import os
 src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
 sys.path.insert(0, src_path)
 
-from Eq04 import solve_System_4eq
+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
+
 
 # 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
+Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
+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
+phi_right = 0.0
+p_right = 0
+
+# Solver settings
+number_cells = 1024
+relax_param = 0.03
+p_right = 0
+
+
+# phi^L domain
+Vol_start = 0
+Volt_end = 0.75
+n_Volts = 20#0
+
+phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
+
+
+# Numerical calculations
+# Solution vectors
+y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
+for i, phi_bcs in enumerate(phi_left):
+    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, refinement_style='hard_log', rtol=1e-4, max_iter=2500, return_type='Vector', relax_param=relax_param)
+    y_S_ = 1 - y_A_ - y_C_
+    y_A_num.append(y_A_)
+    y_C_num.append(y_C_)
+    y_S_num.append(y_S_)
+    phi_num.append(phi_)
+    p_num.append(p_)
+    x_num.append(x_)
+    
+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 = np.array(Q_num)
+
+dx_ = phi_left[1] - phi_left[0] # [1/V], Assumption: phi^L is uniformly distributed
+C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
+C_DL_num = np.array(C_DL_num)
+C_dl_num = C_dl(Q_num, phi_left)
+
+
+
+# Analytical calculations
+Q_ana = Q_DL_dimless_ana(y_R, y_R, 1-2*y_R, z_A, z_C, phi_left, phi_right, p_right, Lambda2, a2, kappa)
+C_DL_ana = C_dl(Q_ana, phi_left)
+
+
+plt.figure()
+plt.plot(phi_left, Q_num - Q_ana, label='Difference')
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('$Q_{num} - Q_{ana} [-]$')
+plt.tight_layout()
+plt.show()   
 
-# Define Parameter and buondary conditions
+plt.figure()
+plt.plot(Phi_pot_center(phi_left), C_DL_num - C_DL_ana, label='Difference')
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('$C_{dl,num} - C_{dl,ana} [-]$')
+plt.tight_layout()
+plt.show()
\ No newline at end of file

src/Helper_DoubleLayerCapacity.py

+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+import numpy as np
+
+
+# Helper functions
+def Phi_pot_center(Phi_pot:np.ndarray) -> np.ndarray:
+    '''
+    Returns vector with the center of the electric potential values.
+
+    Parameters
+    ----------
+    Phi_pot : np.ndarray
+        Input vector with electric potential values.
+
+    Returns
+    -------
+    np.ndarray
+        Vector with the center of the electric potential values.
+    '''
+    return (Phi_pot[1:] + Phi_pot[:-1]) / 2
+
+def dx(Phi_pot:np.ndarray) -> float:
+    '''
+    Returns the difference between the first two electric potential values.
+    
+    Assumes that the electric potential values are equally spaced.
+
+    Parameters
+    ----------
+    Phi_pot : np.ndarray
+        Input vector with electric potential values.
+
+    Returns
+    -------
+    float
+        Difference between the first two electric potential values.
+    '''
+    return Phi_pot[1] - Phi_pot[0]
+
+def C_dl(Q_DL:np.ndarray, Phi_pot:np.ndarray) -> np.ndarray:
+    '''
+    Double Layer Capacity
+
+    Parameters
+    ----------
+    Q_DL : np.ndarray
+        Charge of the system
+    Phi_pot : np.ndarray
+        Electric potential values
+
+    Returns
+    -------
+    np.ndarray
+        Double Layer Capacity
+    '''
+    return (Q_DL[1:] - Q_DL[:-1]) / dx(Phi_pot)
+
+def n_num(p:np.ndarray, K:str|float) -> np.ndarray:
+    '''
+    Calculates the total number density
+
+    Parameters
+    ----------
+    p : np.ndarray
+        Pressure
+    K : str | float
+        Bulk modulus, use 'incompressible' for an incompressible mixture
+
+    Returns
+    -------
+    np.ndarray
+        Total number density
+    '''
+    if K == 'incompressible':
+        return np.ones_like(p)
+    n_Dimensionless = (p-1) / K + 1
+    return n_Dimensionless
+
+def 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:
+    '''
+    Calculates the charge of the system
+
+    Q = ∫_Ω n^F dΩ
+
+    Parameters
+    ----------
+    y_A : np.ndarray
+        Anion fraction
+    y_C : np.ndarray
+        Cation fraction
+    n : np.ndarray
+        Total number density
+    x : np.ndarray
+        Spatial discretization
+    z_A : float, optional
+        Charge number of anions, by default -1.0
+    z_C : float, optional
+        Charge number of cations, by default 1.0
+
+    Returns
+    -------
+    float
+        Charge of the system
+    '''
+    nF_dimensionless = (z_C * y_C + z_A * y_A) * n
+    nF_int = -np.trapz(nF_dimensionless, x)
+    return nF_int
+
+def 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, Lambda2:float, a2:float, kappa:float) -> float:
+    '''
+    Calculates charge of the system using the analytical method in dimensionless units
+
+    Q = sgn(φᴸ −φᴿ)λ√(2/(κ+1))(√(Λᴸ)−√(Λᴿ))
+
+    with 
+    Λ = ln(∑_α D_α exp(−z_α φ + Eₚ a²)
+
+    Parameters
+    ----------
+    y_A_R : float
+        Value of anion fraction at the right boundary
+    y_C_R : float
+        Value of cation fraction at the right boundary
+    y_N_R : float
+        Value of neutral fraction at the right boundary
+    z_A : float
+        Charge number of anions
+    z_C : float
+        Charge number of cations
+    phi_L : float
+        Electric potential at the left boundary
+    phi_R : float
+        Electric potential at the right boundary
+    p_R : float
+        Pressure at the right boundary
+    Lambda2 : float
+        Dimensionless parameter
+    a2 : float
+        Dimensionless parameter
+    kappa : float
+        Solvation number
+
+    Returns
+    -------
+    float
+        Charge of the system in dimensionless units
+    '''
+    z_N = 0
+    D_A = y_A_R / (np.exp(-a2*p_R-z_A*phi_R))
+    D_C = y_C_R / (np.exp(-a2*p_R-z_C*phi_R))
+    D_N = y_N_R / (np.exp(-a2*p_R-z_N*phi_R))
+    E_p = p_R
+    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))
+    CLambda_R = np.log(D_A * np.exp(-z_A * phi_R + E_p * a2) + D_C * np.exp(-z_C * phi_R + E_p * a2) + D_N * np.exp(-z_N * phi_R + E_p * a2))
+    Lambda = np.sqrt(Lambda2)
+    Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2/(kappa+1)) * (np.sqrt(CLambda_L) - np.sqrt(CLambda_R))
+    return Q_DL
+
+def 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, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, kappa:float) -> float:
+    '''
+    Calculates charge of the system using the analytical method in dimensionless units
+
+    Q = sgn(φᴸ −φᴿ)λ√(2/(κ+1))(√(Λᴸ)−√(Λᴿ))
+
+    with 
+    Λ = ln(∑_α D_α exp(−z_α φ + Eₚ a²)
+
+    Parameters
+    ----------
+    y_A_R : float
+        Value of anion fraction at the right boundary
+    y_C_R : float
+        Value of cation fraction at the right boundary
+    y_N_R : float
+        Value of neutral fraction at the right boundary
+    z_A : float
+        Charge number of anions
+    z_C : float
+        Charge number of cations
+    phi_L : float
+        Value of electric potential at the left boundary
+    phi_R : float
+        Value of electric potential at the right
+    p_R : float
+        Value of pressure at the right boundary
+    Lambda2 : float
+        Dimensionless parameter
+    a2 : float
+        Dimensionless parameter
+    nR_m : float
+        Reference number density in 1/m^3
+    e0 : float
+        Dielectric constant
+    LR : float
+        Reference length in m
+    kappa : float
+        Solvation number
+
+    Returns
+    -------
+    float
+        Charge of the system in µAs/cm³
+    '''
+    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, Lambda2, a2, kappa)
+    Q_DL *= nR_m * e0 * LR
+    Q_DL *= 1e+6 
+    Q_DL *= 1/(1e+4)
+    return Q_DL
\ No newline at end of file