DLKap Update for compressible case, solvatin not working, as y_n not coming...

DLKap Update for compressible case, solvatin not working, as y_n not coming from expressions for y_alpha

examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py

@@ -75,12 +75,13 @@ nR_m = nR_mol * NA * 1/(1e-3)# [1/m^3]
 pR = 1.01325 * 1e+5 # [Pa]
 LR = 20e-8
 chi = 80 # [-]
-
+K_vec = ['incompressible', 50_000, 15_000, 1_500, 500] # Working
+# ! Not working: 100_000, 20_000
 # 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 = 1_000#'incompressible'
 kappa = 0
+K = 50_000
 Molarity = 0.01
 y_R = Molarity / nR_mol
 z_A, z_C = -1.0, 1.0
@@ -103,26 +104,26 @@ phi_left_dimless = 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_dimless):
-#     y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, 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_)
+y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
+for i, phi_bcs in enumerate(phi_left_dimless):
+    y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, 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_dimless)):
-#     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)
+Q_num = []
+for j in range(len(phi_left_dimless)):
+    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_dimless[1] - phi_left_dimless[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_dimless)
+dx_ = phi_left_dimless[1] - phi_left_dimless[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_dimless)
 
 # Analytical calculations
 # K = K_vec[0]
@@ -137,7 +138,7 @@ plt.figure()
 # plt.plot(phi_left, Q_num[1] - Q_ana, label='Difference')
 # plt.plot(phi_left_dimless, Q_num - Q_ana, label='Difference')
 # plt.plot(phi_left_dimless, Q_num[1], label='Numerical')
-# plt.plot(phi_left_dimless, Q_num, label='Numerical')
+plt.plot(phi_left_dimless, Q_num, label='Numerical')
 plt.plot(phi_left_dimless, Q_ana, label='Analytical')
 plt.grid()
 plt.legend()
@@ -150,7 +151,7 @@ plt.figure()
 # plt.plot(Phi_pot_center(phi_left), C_DL_num[1] - C_DL_ana, label='Difference')
 # plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num - C_DL_ana, label='Difference')
 # plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num[1], label='Numerical')
-# plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num, label='Numerical')
+plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num, label='Numerical')
 plt.plot(Phi_pot_center(phi_left_dimless), C_DL_ana, label='Analytical')
 plt.grid()
 plt.legend()

examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py

@@ -44,7 +44,7 @@ chi = 80 # [-]
 Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
 a2 = (pR)/(nR_m * k * T)
 K = 'incompressible'
-kappa = 0
+kappa = 5
 Molarity = 0.01
 y_R = Molarity / nR_mol
 z_A, z_C = -1.0, 1.0
@@ -61,7 +61,7 @@ rtol = 1e-4 # ! Change back to 1e-8
 # phi^L domain
 Vol_start = 0.1 # ! Change back to 0
 Volt_end = 0.75
-n_Volts = 5#0
+n_Volts = 10#0
 
 phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
 
@@ -99,7 +99,7 @@ C_DL_ana = C_dl(Q_ana, phi_left)
 # Plotting
 plt.figure()
 # plt.plot(phi_left, Q_num - Q_ana, label='Difference')
-# plt.plot(phi_left, Q_num, label='Numerical')
+plt.plot(phi_left, Q_num, label='Numerical')
 plt.plot(phi_left, Q_ana, label='Analytical')
 plt.grid()
 plt.legend()
@@ -110,7 +110,7 @@ plt.show()
 
 plt.figure()
 # plt.plot(Phi_pot_center(phi_left), C_DL_num - C_DL_ana, label='Difference')
-# plt.plot(Phi_pot_center(phi_left), C_dl_num, label='Numerical')
+plt.plot(Phi_pot_center(phi_left), C_dl_num, label='Numerical')
 plt.plot(Phi_pot_center(phi_left), C_DL_ana, label='Analytical')
 plt.grid()
 plt.legend()
@@ -120,5 +120,5 @@ plt.tight_layout()
 plt.show()
 
 print('phi_left:', phi_left)
-print('Q_num:', Q_num)
+# print('Q_num:', Q_num)
 print('Q_ana:', Q_ana)
\ No newline at end of file

examples/ReproducableCode/ValidationSimplifiedSystem.py

@@ -22,7 +22,6 @@ del sys.path[0]
 # Further imports
 import matplotlib.pyplot as plt
 import numpy as np
-import pandas as pd
 
 # Define the parameters and boundary conditions
 phi_left = 4.0
@@ -36,59 +35,11 @@ Lambda2 = 8.553e-6
 a2 = 7.5412e-4
 solvation = 15
 number_cells = 1024#*4
-refinement_style = 'log'
+refinement_style = 'hard_log'
 rtol = 1e-8
 
-# # Incompressible
-# K = 'incompressible'
-
-# # solve the complete system
-# y_A_4eq, y_C_4eq, phi_4eq, p_4eq, x_4eq = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
-
-# # solve the simplified system
-# y_A_2eq, y_C_2eq, phi_2eq, p_2eq, x_2eq = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
-
-# # Evaluate the difference
-# plt.figure()
-# # plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
-# # plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
-# plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
-# plt.grid()
-# plt.xlim(0, 0.05)
-# plt.legend()
-# plt.show()
-
-# plt.figure()
-# # plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
-# # plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
-# plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
-# plt.grid()
-# plt.xlim(0, 0.05)
-# plt.legend()
-# plt.show()
-
-# plt.figure()
-# # plt.plot(x_4eq, phi_4eq, label='phi_4eq')
-# # plt.plot(x_2eq, phi_2eq, label='phi_2eq')
-# plt.plot(x_4eq, phi_4eq - phi_2eq, label='phi_4eq - phi_2eq')
-# plt.grid()
-# plt.xlim(0, 0.05)
-# plt.legend()
-# plt.show()
-
-# plt.figure()
-# # plt.plot(x_4eq, p_4eq, label='p_4eq')
-# # plt.plot(x_2eq, p_2eq, label='p_2eq')
-# plt.plot(x_4eq, p_4eq - p_2eq, label='p_4eq - p_2eq')
-# plt.grid()
-# plt.xlim(0, 0.05)
-# plt.legend()
-# plt.show()
-
-
-
-# Compressible
-K = 10_000
+# Incompressible
+K = 'incompressible'
 
 # solve the complete system
 y_A_4eq, y_C_4eq, phi_4eq, p_4eq, x_4eq = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
@@ -98,20 +49,20 @@ y_A_2eq, y_C_2eq, phi_2eq, p_2eq, x_2eq = solve_System_2eq(phi_left, phi_right,
 
 # Evaluate the difference
 plt.figure()
-plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
-plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
-# plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
+# plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
+# plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
+plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
 plt.grid()
-plt.xlim(0, 0.1)
+plt.xlim(0, 0.05)
 plt.legend()
 plt.show()
 
 plt.figure()
-plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
-plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
-# plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
+# plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
+# plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
+plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
 plt.grid()
-plt.xlim(0, 0.1)
+plt.xlim(0, 0.05)
 plt.legend()
 plt.show()
 
@@ -120,7 +71,7 @@ plt.figure()
 # plt.plot(x_2eq, phi_2eq, label='phi_2eq')
 plt.plot(x_4eq, phi_4eq - phi_2eq, label='phi_4eq - phi_2eq')
 plt.grid()
-plt.xlim(0, 0.1)
+plt.xlim(0, 0.05)
 plt.legend()
 plt.show()
 
@@ -129,6 +80,54 @@ plt.figure()
 # plt.plot(x_2eq, p_2eq, label='p_2eq')
 plt.plot(x_4eq, p_4eq - p_2eq, label='p_4eq - p_2eq')
 plt.grid()
-plt.xlim(0, 0.1)
+plt.xlim(0, 0.05)
 plt.legend()
-plt.show()
+plt.show()
\ No newline at end of file
+
+
+
+# # Compressible
+# K = 10_000
+
+# # solve the complete system
+# y_A_4eq, y_C_4eq, phi_4eq, p_4eq, x_4eq = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
+
+# # solve the simplified system
+# y_A_2eq, y_C_2eq, phi_2eq, p_2eq, x_2eq = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
+
+# # Evaluate the difference
+# plt.figure()
+# plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
+# plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
+# # plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
+# plt.grid()
+# plt.xlim(0, 0.1)
+# plt.legend()
+# plt.show()
+
+# plt.figure()
+# plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
+# plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
+# # plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
+# plt.grid()
+# plt.xlim(0, 0.1)
+# plt.legend()
+# plt.show()
+
+# plt.figure()
+# # plt.plot(x_4eq, phi_4eq, label='phi_4eq')
+# # plt.plot(x_2eq, phi_2eq, label='phi_2eq')
+# plt.plot(x_4eq, phi_4eq - phi_2eq, label='phi_4eq - phi_2eq')
+# plt.grid()
+# plt.xlim(0, 0.1)
+# plt.legend()
+# plt.show()
+
+# plt.figure()
+# # plt.plot(x_4eq, p_4eq, label='p_4eq')
+# # plt.plot(x_2eq, p_2eq, label='p_2eq')
+# plt.plot(x_4eq, p_4eq - p_2eq, label='p_4eq - p_2eq')
+# plt.grid()
+# plt.xlim(0, 0.1)
+# plt.legend()
+# plt.show()
\ No newline at end of file

src/Helper_DoubleLayerCapacity.py

@@ -4,6 +4,7 @@ Jan.Habscheid@rwth-aachen.de
 '''
 
 import numpy as np
+from scipy.optimize import fixed_point, fsolve
 
 
 # Helper functions
@@ -157,66 +158,51 @@ def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float
         # Assume E_p = p_R = 0        
         D_A = y_A_R / (np.exp(-(solvation+1)*a2*p_R-z_A*phi_R))
         D_C = y_C_R / (np.exp(-(solvation+1)*a2*p_R-z_C*phi_R))
-        D_N = y_N_R / (np.exp(-(solvation+1)*a2*p_R-z_N*phi_R))
-
-        D_tilde_A = np.log(D_A)
-        D_tilde_C = np.log(D_C)
-        D_tilde_N = np.log(D_N)
-
-        CLambda_L = - z_A * phi_L + D_tilde_A - z_C * phi_L + D_tilde_C - z_N * phi_L + D_tilde_N
-        CLambda_R = - z_A * phi_R + D_tilde_A - z_C * phi_R + D_tilde_C - z_N * phi_R + D_tilde_N
-
-        Lambda = np.sqrt(Lambda2)
-        a = np.sqrt(a2)
-
-        N = 3
-
-        Q_DL = Lambda * a * (np.sqrt(-CLambda_L/(-(N-1)*(solvation+1)-1)) - np.sqrt(-CLambda_R/(-(N-1)*(solvation+1)-1)))
-
-        # CLambda_L = np.log(D_A * np.exp(-z_A * phi_L) + D_C * np.exp(-z_C * phi_L) + D_N * np.exp(-z_N * phi_L))
-        # CLambda_R = np.log(D_A * np.exp(-z_A * phi_R) + D_C * np.exp(-z_C * phi_R) + D_N * np.exp(-z_N * phi_R))
-
+        D_N = y_N_R / (np.exp(-a2*p_R))
+
+        Q_DL = []
+        dx_phi_L = 0
+        for phi_L_ in phi_L:
+            def func(p_L):
+                # A_Term = D_A * np.exp(-z_A*phi_L)
+                A_Term = D_A * np.exp(-(solvation+1)*a2*p_L-z_A*phi_L_)
+                C_Term = D_C * np.exp(-(solvation+1)*a2*p_L-z_C*phi_L_)
+                N_Term = D_N * np.exp(-a2*p_L)
+
+                return (A_Term + C_Term + N_Term) - 1
+            p_L = fsolve(func, dx_phi_L)
+            dx_phi_L = np.sqrt(p_L - E_p) * np.sqrt(2*a2/Lambda2)
+            Q_DL_ = Lambda2 * dx_phi_L
+            Q_DL.append(Q_DL_)
+        Q_DL = np.array(Q_DL)
+
+        # def func_L(p_L):
+        #     CLambda = D_A * np.exp(-z_A * phi_L - solvation*a2*p_L) + D_C * np.exp(-z_C * phi_L - solvation*a2*p_L) + D_N
+        #     return np.log(CLambda)/a2
+        
+        # p_L = fixed_point(func_L, 0, maxiter=5_000)
         # Lambda = np.sqrt(Lambda2)
-        # Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2) * (np.sqrt(CLambda_L/(solvation+1) - E_p * a2) - np.sqrt(CLambda_R/(solvation+1) - E_p * a2))
-    else:
-        C_A = y_A_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_A*phi_R))
-        C_C = y_C_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_C*phi_R))
-        C_N = y_N_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_N*phi_R))
-
-        # # Lambda_tilda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
-        # # Lambda_tilda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
-
-        # # # Lambda_hat_L = np.exp(np.log(Lambda_tilda_L) / (a2 * K))
-        # # # Lambda_hat_R = np.exp(np.log(Lambda_tilda_R) / (a2 * K))
-
-        # # dx_Phi_L = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_L))
-        # # dx_Phi_R = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_R))
-
-        # # Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
-
-        # Lambda_tilde_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
-        # Lambda_tilde_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
-
-        # Lambda_hat_L = (1/Lambda_tilde_L)**(-1/((solvation+1)*a2*K))
-        # Lambda_hat_R = (1/Lambda_tilde_R)**(-1/((solvation+1)*a2*K))
-
         # a = np.sqrt(a2)
-        # Lambda = np.sqrt(Lambda2)
-        # K_tilde = 1 - K + E_p
+        # Q_DL = np.sqrt(2) * Lambda * a * np.sqrt(p_L - E_p)
 
-        # Q_DL = Lambda * a * np.sqrt(2) * (np.sqrt(Lambda_hat_L - K_tilde) - np.sqrt(Lambda_hat_R - K_tilde))
+        # ! ToDo: Implement full analytical for solvation == 0
+        
+    else:
+        C_A = y_A_R / ((K + p_R  - 1)**(-1*a2*K)*np.exp(-z_A*phi_R))
+        C_C = y_C_R / ((K + p_R  - 1)**(-1*a2*K)*np.exp(-z_C*phi_R))
+        C_N = y_N_R / ((K + p_R  - 1)**(-1*a2*K)*np.exp(-z_N*phi_R))
 
-        Lambda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
-        Lambda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
+        CLambda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N
+        # CLambda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N
 
-        B = -a2 * K
-        A = 1/2 * Lambda2/a2
-        K_tilde = K + E_p - 1
+        Left = np.sqrt((1/CLambda_L)**(-1/(a2*K)) + 1 - K - E_p)
+        # Right = np.sqrt((1/CLambda_R)**(-1/(a2*K)) + 1 - K - E_p)
 
-        dx_Phi_L = (np.power(1/Lambda_L, B) - K_tilde) / A 
-        dx_Phi_R = (np.power(1/Lambda_R, B) - K_tilde) / A
+        Lambda = np.sqrt(Lambda2)
+        a = np.sqrt(a2)
 
-        Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
+        Q_DL = Lambda * a * np.sqrt(2) * (Left)# - Right)
+                     
     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, K:str|float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, solvation:float) -> float: