diff --git a/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py b/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
index f2cfffa6ce66efb6811e721d5cce97ea18a13d14..c17780a490e4e75201db5ab285f5e42d77e07d05 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
@@ -36,7 +36,7 @@ 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
+LR = 20e-9
chi = 80 # [-]
# Parameter and bcs for the electrolyte
@@ -54,12 +54,13 @@ p_right = 0
number_cells = 1024
relax_param = 0.03
p_right = 0
+rtol = 1e-4 # ! Change back to 1e-8
# phi^L domain
-Vol_start = 0
+Vol_start = 0.1 # ! Change back to 0
Volt_end = 0.75
-n_Volts = 5#0
+n_Volts = 25#0
phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
@@ -68,7 +69,7 @@ phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
# 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, solvation=kappa, refinement_style='hard_log', rtol=1e-10, max_iter=2500, return_type='Vector', relax_param=relax_param)
+ 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, solvation=kappa, refinement_style='hard_log', rtol=rtol, 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_)
@@ -82,9 +83,9 @@ for j in range(len(phi_left)):
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
-# C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
-# C_DL_num = np.array(C_DL_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)
@@ -115,4 +116,8 @@ 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
+plt.show()
+
+print('phi_left:', phi_left)
+print('Q_num:', Q_num)
+print('Q_ana:', Q_ana)
\ No newline at end of file
diff --git a/examples/ReproducableCode/ValidationSimplifiedSystem.py b/examples/ReproducableCode/ValidationSimplifiedSystem.py
new file mode 100644
index 0000000000000000000000000000000000000000..d8a2d835df8ee8ef65c2d3c00197c22d9b4c2a19
--- /dev/null
+++ b/examples/ReproducableCode/ValidationSimplifiedSystem.py
@@ -0,0 +1,84 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+
+This script is used to validate the simplification of the system of four equations to a system of two equations in the one-dimensional equilibrium case.
+'''
+
+# 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 Eq04 import solve_System_4eq
+from Eq02 import solve_System_2eq
+
+# Remove the src directory from sys.path after import
+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
+phi_right = 0.0
+p_right = 0.0
+y_A_R = 0.01
+y_C_R = 0.01
+z_A = -1.0
+z_C = 1.0
+K = 'incompressible'
+Lambda2 = 8.553e-6
+a2 = 7.5412e-4
+solvation = 15
+number_cells = 1024#*4
+refinement_style = 'log'
+rtol = 1e-8
+
+# 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()
\ No newline at end of file
diff --git a/src/Eq02.py b/src/Eq02.py
index 45901a10e5dd669a394fb25cb7aa9c36e061a452..7ecd7d0e72574f5207edd068d100190aeb56e7c9 100644
--- a/src/Eq02.py
+++ b/src/Eq02.py
@@ -8,7 +8,7 @@ from mpi4py import MPI
from dolfinx import mesh, fem, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
-from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Mesh
+from ufl import TestFunctions, split, dot, grad, dx, inner, Mesh, exp
from basix.ufl import element, mixed_element
import matplotlib.pyplot as plt
@@ -51,7 +51,7 @@ def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
return msh
-def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float, z_C:float, y_A_R:float, y_C_R:float, K:float|str, Lambda2:float, a2:float, number_cells:int, solvation:float = 0, PoissonBoltzmann:bool=False, relax_param:float=None, x0:float=0, x1:float=1, refinement_style:str='uniform', return_type:str='Scalar', rtol:float=1e-8, max_iter:float=500):
+def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float, z_C:float, y_A_R:float, y_C_R:float, K:float|str, Lambda2:float, a2:float, number_cells:int, solvation:float = 0, relax_param:float=None, x0:float=0, x1:float=1, refinement_style:str='uniform', return_type:str='Scalar', rtol:float=1e-8, max_iter:float=500):
'''
Solve the simplified dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
@@ -96,8 +96,6 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
Number of cells in the mesh
solvation : float, optional
solvation number, by default 0
- PoissonBoltzmann : bool, optional
- Solve classical Nernst-Planck model with the use of the Poisson-Boltzmann formulation if True, else solve the presented model by Dreyer, Guhlke, Müller, by default False
relax_param : float, optional
Relaxation parameter for the Newton solver
xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
@@ -122,6 +120,8 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
If return_type is 'Vector', the solution is returned as numpy arrays
'''
+ if return_type == 'Scalar':
+ raise NotImplementedError('Scalar return type is not implemented yet')
# Define boundaries of the domain
x0 = 0
x1 = 1
@@ -143,20 +143,18 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
CG1_elem = element('Lagrange', msh.basix_cell(), 1)
# Define Mixed Function Space
- W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
+ W_elem = mixed_element([CG1_elem, CG1_elem])#, CG1_elem, CG1_elem])
W = fem.functionspace(msh, W_elem)
# Define Trial- and Testfunctions
u = fem.Function(W)
- y_A, y_C, phi, p = split(u)
- (v_A, v_C, v_1, v_2) = TestFunctions(W)
+ phi, p = split(u)
+ (v_1, v_2) = TestFunctions(W)
# Collapse function space for bcs
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()
- W2, _ = W.sub(2).collapse()
- W3, _ = W.sub(3).collapse()
-
+
# Define boundary conditions values
def phi_left_(x):
return np.full_like(x[0], phi_left)
@@ -164,42 +162,36 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
return np.full_like(x[0], phi_right)
def p_right_(x):
return np.full_like(x[0], p_right)
- def y_A_right_(x):
- return np.full_like(x[0], y_A_R)
- def y_C_right_(x):
- return np.full_like(x[0], y_C_R)
-
+
# Interpolate bcs functions
- phi_left_bcs = fem.Function(W2)
+ phi_left_bcs = fem.Function(W0)
phi_left_bcs.interpolate(phi_left_)
- phi_right_bcs = fem.Function(W2)
+ phi_right_bcs = fem.Function(W0)
phi_right_bcs.interpolate(phi_right_)
- p_right_bcs = fem.Function(W3)
+ p_right_bcs = fem.Function(W1)
p_right_bcs.interpolate(p_right_)
- y_A_right_bcs = fem.Function(W0)
- y_A_right_bcs.interpolate(y_A_right_)
- y_C_right_bcs = fem.Function(W1)
- y_C_right_bcs.interpolate(y_C_right_)
-
+
# Identify dofs for boundary conditions
# Define boundary conditions
- facet_left_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Left)
- facet_right_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Right)
- bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(2))
- bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(2))
-
- facet_right_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Right)
- bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(3))
-
+ facet_left_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Left)
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
- bc_right_y_A = fem.dirichletbc(y_A_right_bcs, facet_right_dofs, W.sub(0))
+ bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(0))
+ bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(0))
facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
- bc_right_y_C = fem.dirichletbc(y_C_right_bcs, facet_right_dofs, W.sub(1))
+ bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(1))
+
# Combine boundary conditions into list
- bcs = [bc_left_phi, bc_right_phi, bc_right_p, bc_right_y_A, bc_right_y_C]
-
+ bcs = [bc_left_phi, bc_right_phi, bc_right_p]
+
+ def y_A(phi, p):
+ D_A = y_A_R / exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
+ return D_A * exp(-(solvation + 1) * a2 * p - z_A * phi)
+
+ def y_C(phi, p):
+ D_C = y_C_R / exp(-(solvation + 1) * a2 * p_right - z_C * phi_right)
+ return D_C * exp(-(solvation + 1) * a2 * p - z_C * phi)
# Define variational problem
@@ -207,30 +199,6 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
# total free charge density
def nF(y_A, y_C):
return (z_C * y_C + z_A * y_A)
-
- # Diffusion fluxes for species A and C
- def J_A(y_A, y_C, phi, p):
- return ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi
-
- def J_C(y_A, y_C, phi, p):
- return ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi
-
- # Variational Form
- A = (
- inner(grad(phi), grad(v_1)) * dx
- - 1 / Lambda2 * nF(y_A, y_C) * v_1 * dx
- ) + (
- inner(grad(p), grad(v_2)) * dx
- + 1 / a2 * nF(y_A, y_C) * dot(grad(phi), grad(v_2)) * dx
- ) + (
- inner(grad(J_A(y_A, y_C, phi, p)), grad(v_A)) * dx
- + inner(grad(J_C(y_A, y_C, phi, p)), grad(v_C)) * dx
- )
- if PoissonBoltzmann:
- A += (
- inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_A)) * dx
- + inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_C)) * dx
- )
else:
# total number density
def n(p):
@@ -239,37 +207,16 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
# total free charge density
def nF(y_A, y_C, p):
return (z_C * y_C + z_A * y_A) * n(p)
-
- # Diffusion fluxes for species A and C
- def J_A(y_A, y_C, phi, p):
- return ln(y_A) + a2 * (solvation + 1) * K * ln(1 + 1/K * (p-1)) + z_A * phi
-
- def J_C(y_A, y_C, phi, p):
- return ln(y_C) + a2 * (solvation + 1)* K * ln(1 + 1/K * (p-1)) + z_C * phi
-
- A = (
- inner(grad(phi), grad(v_1)) * dx
- - 1 / Lambda2 * nF(y_A, y_C, p) * v_1 * dx
- ) + (
- inner(grad(p), grad(v_2)) * dx
- + 1 / a2 * nF(y_A, y_C, p) * dot(grad(phi), grad(v_2)) * dx
- ) + (
- inner(grad(J_A(y_A, y_C, phi, p)), grad(v_A)) * dx
- + inner(grad(J_C(y_A, y_C, phi, p)), grad(v_C)) * dx
- )
+ # Variational Form
+ A = (
+ inner(grad(phi), grad(v_1)) * dx
+ - 1 / Lambda2 * nF(y_A(phi, p), y_C(phi, p)) * v_1 * dx
+ ) + (
+ inner(grad(p), grad(v_2)) * dx
+ + 1 / a2 * nF(y_A(phi, p), y_C(phi, p)) * dot(grad(phi), grad(v_2)) * dx
+ )
F = A
- # Initialize initial guess for u
- y_C_init = fem.Function(W1)
- y_A_init = fem.Function(W0)
- y_C_init.interpolate(lambda x: np.full_like(x[0], y_C_R))
- y_A_init.interpolate(lambda x: np.full_like(x[0], y_A_R))
-
- with u.vector.localForm() as u_loc:
- u_loc.set(0)
- u.sub(0).interpolate(y_A_init)
- u.sub(1).interpolate(y_C_init)
-
# Define Nonlinear Problem
problem = NonlinearProblem(F, u, bcs=bcs)
@@ -294,20 +241,22 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
print(f"Number of interations: {n:d}")
# Split the mixed function space into the individual components
- y_A, y_C, phi, p = u.split()
+ phi, p = u.split()
# Return the solution
if return_type=='Vector':
x_vals = np.array(msh.geometry.x[:,0])
- y_A_vals = np.array(u.sub(0).collapse().x.array)
- y_C_vals = np.array(u.sub(1).collapse().x.array)
- phi_vals = np.array(u.sub(2).collapse().x.array)
- p_vals = np.array(u.sub(3).collapse().x.array)
+ phi_vals = np.array(u.sub(0).collapse().x.array)
+ p_vals = np.array(u.sub(1).collapse().x.array)
+
+ # Calculate the atomic fractions
+ D_A = y_A_R / np.exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
+ y_A_vals = D_A * np.exp(-(solvation + 1) * a2 * p_vals - z_A * phi_vals)
+
+ D_C = y_C_R / np.exp(-(solvation + 1) * a2 * p_right - z_C * phi_right)
+ y_C_vals = D_C * np.exp(-(solvation + 1) * a2 * p_vals - z_C * phi_vals)
return y_A_vals, y_C_vals, phi_vals, p_vals, x_vals
- elif return_type=='Scalar':
- return y_A, y_C, phi, p, msh
-
if __name__ == '__main__':
# Define the parameters
diff --git a/src/Helper_DoubleLayerCapacity.py b/src/Helper_DoubleLayerCapacity.py
index d8bd6970498e7bb1f73f9a367cfd10f903346359..4da01ebab88716abe73041b31c0f973a3470e9c4 100644
--- a/src/Helper_DoubleLayerCapacity.py
+++ b/src/Helper_DoubleLayerCapacity.py
@@ -157,10 +157,15 @@ def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float
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))
E_p = p_R
- CLambda_L = np.log(D_A * np.exp(-z_A * phi_L - (solvation+1) * E_p) + D_C * np.exp(-z_C * phi_L - (solvation+1) * E_p) + D_N * np.exp(-z_N * phi_L - (solvation+1) * E_p))
- 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))
+
+ 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))
+
+ # dx_Phi_L = np.sqrt((CLambda_L - (solvation+1) * a2 * E_p) * 2 / (solvation + 1))
+ # dx_Phi_R = np.sqrt((CLambda_R - (solvation+1) * a2 * E_p) * 2 / (solvation + 1))
+ # Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
Lambda = np.sqrt(Lambda2)
- Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2/(solvation+1)) * (np.sqrt(CLambda_L) - np.sqrt(CLambda_R))
+ Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2/(solvation+1)) * (np.sqrt(CLambda_L + E_p * a2) - np.sqrt(CLambda_R + 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))
diff --git a/src/__pycache__/Eq02.cpython-312.pyc b/src/__pycache__/Eq02.cpython-312.pyc
index 0a43f005f31a0dc1160dbe3957614a55dd0419bd..99acdd38417e0f045265675a551c6f42a3df0a2d 100644
Binary files a/src/__pycache__/Eq02.cpython-312.pyc and b/src/__pycache__/Eq02.cpython-312.pyc differ
diff --git a/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc b/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc
index 0520aeca0b59a17e6054ae07e36af289020f6477..ec8a8993fff2cb48986394c354c67b2ed42ec555 100644
Binary files a/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc and b/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc differ