# Finite Element Method - 3D Mesh Generator - Metfem3D

## Results

### The electrodiffusion equation

The system of coupled equations describing process of electrodiffusion (\ref{electrodiffusion}) written in terms of the FEM method (see Eq.~(\ref{pnp_final})) where k+ = k- = k, D+ = D- = D and the time step equals 0.01 has been numerically solved up to the final time Tend by using the Newton algorithm. The parameter k denotes (ke)/(ε0ε) where k = (De)/(kBT). Computations have been performed for below-given sets of boundary conditions and parameters. The first considered case shows boundary layers in x, y - directions. The second one presents the range of k values in which the one stable solution to this nonlinear problem has been obtained for the assumed values of D = 0.05 and Δt = 0.01. Moreover, solution's convergence depends on h0 - the mesh parameter and the order of FEM approximation. The third case is computed for the physical values of parameters which locate solutions within the stable regime.

• The boundary values of n+, n-, φ are set as 1 at x = 0, x = π, y = 0, y = π, and z = 0 and equal 2 at z = π. Parameters take the values: D = 0.05, k = 0.05, and Tend = 0.39 [units]. The chosen mesh size h0 is 0.44 and the ,,linear'' approximation has been used.
• The boundary values of n+, n-, φ are set as 0 at y = 0, y = π, z = 0, z = π, and x = 0 and equal 1 at x = π. Parameters take the values: D = 0.05, k varies from 0.04 up to 0.43, and Tend = 3.99 [units]. The computations are done on uniform meshes of h0 = 0.44 (k ≤ 0.43) and 0.35 (k ≤ 0.04) in the ,,linear'' FEM approximation and of h0 = 0.55 (k ≤ 0.04) in the ,,quadratic'' one. This analysis forms some preliminary study of the stable and a chaotic behavior of the considered system of nonlinear equations. A closer insight into the problem goes beyond this article.
• The boundary values of n+, n- are set as 0 at y = 0, y = π m, z = 0, z = π m, and x = 0 and equal 6 ✗ 107 [n.u.] at x = π m where 1 n.u. = 1019 ions in m3 (whereas 6 ✗ 1026 ions/m3 gives a concentration circa 1M). The applied electric potential φ is set in Volts. It only at x = π m differs from 0 and equals 6V. Parameters take the values: D = 2 ✗ 10-9 m2/s, k = 1.76 ✗ 10-17 ≈ D ✗ 10-8 m3/s (and are the physical values for ions K+ and Cl- in the room temperature T = 298 K). The time step Δt equals 0.01 s and Tend = 3.99 s. The chosen mesh size h0 is 0.35 and the ,,linear'' approximation has been used.

In all cases an initial guess of n+, n- and φ distributions has been chosen as 0 everywhere in the domain apart from its boundaries.
Fig. 7 presents obtained profiles of cations n+ and the potential φ at the center of the domain i. e. at z = π/2 and at the final time Tend computed for the first and the third case. There is no visible difference between cations n+ and anions n- distributions in both cases so the latter is not shown.
In the first case, the maximum of n+, i+1 - n+, i is decaying from 0.023 for i = 0 (t = 0) to 0.0093 for i = 38 (t = 0.38) and in the latter from 0.17 (✗ 107 [n.u.]) at t = 0 to 0.0033 (✗ 107 [n.u.]) for i = 398 (t = 3.98 s). The system of particles tends to its stationary state. Moreover, the maximum of difference between n+ and n- distributions computed at each node at the time Tend equals 1.3e-09 in the first case and 2.3e-10 (✗ 107 [n.u.]) in the latter one. Employing physical notion it means that the system of charged particles is electroneutral.

The picture shows FEM approximation (linear) of the system of electrodiffusion equations being the boundary value problem. The solutions of a) and c) the distribution of cations n+ and b) and d) the profile of the potential φ are depicted in the center of the cubic domain i. e. at z = π/2 and at times Tend = 0.39 with the values of parameters Δt = 0.01 [units], k+ = k- = 0.05 [units], D+ = D- = 0.05 [units] (upper cases) and the physical ones (lower cases), respectively; computations have been performed on the uniform mesh with the volume of tetrahedron V0 = 0.01 and the element size h0 = 0.44 (upper cases) and V0 = 0.005 for h0 = 0.35 (lower ones).