Source Code

Warning

If the Newton solver diverges for any of the solutions, you may try to reduce the relaxation parameter.

Eq04.create_refined_mesh(refinement_style: str, number_cells: int) Mesh

Creates a one-dimensional mesh with a refined region at the left boundary

Parameters:
  • refinement_style (str) – How the mesh should be refined. Options are ‘log’, ‘hard_log’, ‘hard_hard_log’

  • number_cells (int) – Number of cells in the mesh

Returns:

One-dimensional mesh, ready for use in FEniCSx

Return type:

Mesh

EqN.solve_System_Neq(phi_left: float, phi_right: float, p_right: float, z_alpha: list, y_R: list, 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 = 'Vector', rtol: float = 1e-08, max_iter: float = 500)

Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024

System of equations:

λ²Δ φ =−L²n^F

a²∇p=−n^F∇ φ

div(J_α)=0  α∈ {1,…,N−1}

with φ the electric potential, p the pressure, n^F the total free charge density, J_α the diffusion fluxes of species α, λ² a dimensionless parameter, L²=1, a² a dimensionless parameter, N the number of species, and α the species index.

! If the Newton solver diverges, you may try to reduce the relaxation parameter.

Parameters:
  • phi_left (float) – Value of φ at the left boundary

  • phi_right (float) – Value of φ at the right boundary

  • p_right (float) – Value of p at the right boundary

  • z_alpha (list) – Charge numbers for species α = 1,…,N-1

  • y_R (list) – Atomic fractions at right boundary for species α = 1,…,N-1

  • K (float | str) – Dimensioness bulk modulus of the electrolyte. If ‘incompressible’, the system is solved for an incompressible electrolyte

  • Lambda2 (float) – Dimensionless parameter

  • a2 (float) – Dimensionless parameter

  • number_cells (int) – Number of cells in the mesh

  • solvation (float, optional) – solvation number, not implemented yet, 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, Not implemented yet, by default False

  • relax_param (float, optional) – Relaxation parameter for the Newton solver xₙ₊₁ = γ xₙ f(xₙ)/f’(xₙ) with γ the relaxation parameter , by default None -> Determined automatically

  • x0 (float, optional) – Left boundary of the domain, by default 0

  • x1 (float, optional) – Right boundary of the domain, by default 1

  • refinement_style (str, optional) – Specify for refinement towards zero Options are ‘uniform’, ‘log’, ‘hard_log’, ‘hard_hard_log’ by default ‘uniform’

  • return_type (str, optional) – ‘Vector’ or ‘Scalar’ (not implemented yet, should be implemented in a later version), ‘Scalar’ returns dolfinx.fem type and ‘Vector’ numpy arrays of the solution, by default ‘Vector’

  • rtol (float, optional) – Relative tolerance for Newton solver, by default 1e-8

  • max_iter (float, optional) – Maximum number of Newton iterations, by default 500

Returns:

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 Only return_type ‘Vector’ is implemented yet

Return type:

y_A, y_C, phi, p, msh

Eq04.solve_System_4eq(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-08, max_iter: float = 500)

Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024

System of equations:

λ²Δ φ =−L²n^F

a²∇p=−n^F∇ φ

div(J_A)=0

div(J_C)=0

with φ the electric potential, p the pressure, n^F the total free charge density, J_α the diffusion fluxes of species α, λ² a dimensionless parameter, L²=1, a² a dimensionless parameter, N the number of species, and α the species index.

! If the Newton solver diverges, you may try to reduce the relaxation parameter.

Parameters:
  • phi_left (float) – Value of φ at the left boundary

  • phi_right (float) – Value of φ at the right boundary

  • p_right (float) – Value of p at the right boundary

  • z_A (float) – Charge number of species A

  • z_C (float) – Charge number of species C

  • y_A_R (float) – Atomic fractions of species A at right boundary

  • y_C_R (float) – Atomic fractions of species C at right boundary

  • K (float | str) – Dimensioness bulk modulus of the electrolyte. If ‘incompressible’, the system is solved for an incompressible electrolyte

  • Lambda2 (float) – Dimensionless parameter

  • a2 (float) – Dimensionless parameter

  • number_cells (int) – 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 , by default None -> Determined automatically

  • x0 (float, optional) – Left boundary of the domain, by default 0

  • x1 (float, optional) – Right boundary of the domain, by default 1

  • refinement_style (str, optional) – Specify for refinement towards zero Options are ‘uniform’, ‘log’, ‘hard_log’, ‘hard_hard_log’ by default ‘uniform’

  • return_type (str, optional) – ‘Vector’ or ‘Scalar’, ‘Scalar’ returns dolfinx.fem type and ‘Vector’ numpy arrays of the solution, by default ‘Scalar’

  • rtol (float, optional) – Relative tolerance for Newton solver, by default 1e-8

  • max_iter (float, optional) – Maximum number of Newton iterations, by default 500

Returns:

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

Return type:

y_A, y_C, phi, p, msh

Eq02.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-08, 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

System of equations:

λ²Δ φ =−L²n^F

a²∇p=−n^F∇ φ

with the space charge

n^F = z_A y_A(φ, p) + z_C y_C(φ ,p)

if the mixture is compressible:

y_alpha = C_alpha * (K+p−1)^(−κ+1)a²K exp(−z_α φ)

if the mixture is incompressible:

y_alpha = D_alpha * exp(−(κ+1)a²p−z_α φ)

with φ the electric potential, p the pressure, n^F the total free charge density, J_α the diffusion fluxes of species α, λ² a dimensionless parameter, L²=1, a² a dimensionless parameter, N the number of species, and α the species index.

! If the Newton solver diverges, you may try to reduce the relaxation parameter.

Parameters:
  • phi_left (float) – Value of φ at the left boundary

  • phi_right (float) – Value of φ at the right boundary

  • p_right (float) – Value of p at the right boundary

  • z_A (float) – Charge number of species A

  • z_C (float) – Charge number of species C

  • y_A_R (float) – Atomic fractions of species A at right boundary

  • y_C_R (float) – Atomic fractions of species C at right boundary

  • K (float | str) – Dimensioness bulk modulus of the electrolyte. If ‘incompressible’, the system is solved for an incompressible electrolyte

  • Lambda2 (float) – Dimensionless parameter

  • a2 (float) – Dimensionless parameter

  • number_cells (int) – Number of cells in the mesh

  • solvation (float, optional) – solvation number, by default 0

  • relax_param (float, optional) – Relaxation parameter for the Newton solver xₙ₊₁ = γ xₙ f(xₙ)/f’(xₙ) with γ the relaxation parameter , by default None -> Determined automatically

  • x0 (float, optional) – Left boundary of the domain, by default 0

  • x1 (float, optional) – Right boundary of the domain, by default 1

  • refinement_style (str, optional) – Specify for refinement towards zero Options are ‘uniform’, ‘log’, ‘hard_log’, ‘hard_hard_log’ by default ‘uniform’

  • return_type (str, optional) – ‘Vector’ or ‘Scalar’, ‘Scalar’ returns dolfinx.fem type and ‘Vector’ numpy arrays of the solution, by default ‘Scalar’

  • rtol (float, optional) – Relative tolerance for Newton solver, by default 1e-8

  • max_iter (float, optional) – Maximum number of Newton iterations, by default 500

Returns:

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

Return type:

y_A, y_C, phi, p, msh

ElectrolyticDiode.ElectrolyticDiode(Bias_type: str, phi_bias: float, g_phi: float, z_A: float, z_C: float, y_A_bath: float, y_C_bath: float, K: float | str, Lambda2: float, a2: float, number_cells: list, solvation: float = 0, PoissonBoltzmann: bool = False, relax_param: float = None, Lx: float = 2, Ly: float = 10, rtol: float = 1e-08, max_iter: float = 500, return_type: str = 'Vector')

Solves the system of equations for the example of an electric diode

Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024

! If the Newton solver diverges, you may try to reduce the relaxation parameter.

Parameters:
  • Bias_type (str) – ForwardBias, NoBias, BackwardBias

  • phi_bias (float) – Bias in φ

  • g_phi (float) – Neumann boundary condition for φ

  • z_A (float) – Charge number of species A

  • z_C (float) – Charge number of species C

  • y_A_bath (float) – Bath concentration of anions

  • y_C_bath (float) – Bath concentration of cations

  • K (float | str) – Dimensioness bulk modulus of the electrolyte. If ‘incompressible’, the system is solved for an incompressible electrolyte ! Only implemented for incompressible mixtures, yet

  • Lambda2 (float) – Dimensionless parameter

  • a2 (float) – Dimensionless parameter

  • number_cells (int) – 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 , by default None -> Determined automatically

  • Lx (float, optional) – Length of domain in x-direction, by default 2

  • Ly (float, optional) – Length of domain in y-direction, by default 10

  • rtol (float, optional) – Relative tolerance for Newton solver, by default 1e-8

  • max_iter (float, optional) – Maximum number of Newton iterations, by default 500

  • return_type (str, optional) – Type of return value, by default ‘Vector’ ‘Vector’ -> Returns the solution as a numpy array for triangle elements ‘Extended’ -> Returns the solution as both, a numpy array for triangle elements and the fenicsx function u

Returns:

Returns atomic fractions for species A and C, electric potential, pressure, and the mesh as numpy arrays for triangle elements. x_vector is a list of both dimensions [x, y] If return_type is ‘Extended’, also returns the fenicsx function u

Return type:

y_A_vector, y_C_vector, phi_vector, p_vector, x_vector

Jan Habscheid Jan.Habscheid@rwth-aachen.de

Helper_DoubleLayerCapacity.C_dl(Q_DL: ndarray, Phi_pot: ndarray) ndarray

Double Layer Capacity

Parameters:
  • Q_DL (np.ndarray) – Charge of the system

  • Phi_pot (np.ndarray) – Electric potential values

Returns:

Double Layer Capacity

Return type:

np.ndarray

Helper_DoubleLayerCapacity.Phi_pot_center(Phi_pot: ndarray) ndarray

Returns vector with the center of the electric potential values.

Parameters:

Phi_pot (np.ndarray) – Input vector with electric potential values.

Returns:

Vector with the center of the electric potential values.

Return type:

np.ndarray

Helper_DoubleLayerCapacity.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

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

  • K (str | float) – Bulk modulus, use ‘incompressible’ for an incompress

  • 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

  • solvation (float) – Solvation number

Returns:

Charge of the system in µAs/cm³

Return type:

float

Helper_DoubleLayerCapacity.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

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 [V]

  • phi_R (float) – Electric potential at the right boundary

  • p_R (float) – Pressure at the right boundary

  • K (str | float) – Bulk modulus, use ‘incompressible’ for an incompressible mixture

  • Lambda2 (float) – Dimensionless parameter

  • a2 (float) – Dimensionless parameter

  • solvation (float) – Solvation number

Returns:

Charge of the system in dimensionless units

Return type:

float

Helper_DoubleLayerCapacity.Q_num_(y_A: ndarray, y_C: ndarray, n: ndarray, x: 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:

Charge of the system

Return type:

float

Helper_DoubleLayerCapacity.Q_num_dim(y_A: ndarray, y_C: ndarray, n: ndarray, x: ndarray, z_A: float, z_C: float, nR_m: float, e0: float, LR: float) 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

  • nR_m (float) – Reference number density in 1/m^3

  • e0 (float) – Dielectric constant

  • LR (float) – Reference length in m

Returns:

Charge of the system

Return type:

float

Helper_DoubleLayerCapacity.dx(Phi_pot: 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:

Difference between the first two electric potential values.

Return type:

float

Helper_DoubleLayerCapacity.n(p: ndarray, K: str | float) ndarray

Calculates the total number density

Parameters:
  • p (np.ndarray) – Pressure

  • K (str | float) – Bulk modulus, use ‘incompressible’ for an incompressible mixture

Returns:

Total number density

Return type:

np.ndarray