Confrontation of the Ohmic approach with the ionic transport approach for modeling the electrical behavior of an electrolyte

Two mathematical approaches have been compared for modeling the electrical behavior (e.g., electric/current density field) of an electrolyte. The first approach uses the traditional Ohmic method (Laplace equation) based on the assumption of a uniform concentration of ions in the electrolyte. The second approach utilizes the ionic transport method that takes into account the effect of electrochemical transport/reaction of ions assuming a quasi-electro neutral bulk for the electrolyte. For the latter, the Poisson–Nernst–Planck (PNP) equations are employed to include impacts of the non-uniformity of the concentration of ions on the electrical pattern of the electrolyte. Demonstratively, a 2D case was studied. It comprises two unequally sized electrodes, which are separated by the electrolyte. Findings show that the electrical behavior of the electrolyte remains unchanged using an Ohmic approach, regardless of variations in operational parameters such as applied voltage and electrode polarity. In contrast, the electrical behavior in response to the operational parameters is found to be significantly influenced using the ionic transport approach. Therefore, Ohm’s law is invalid in the bulk of an electrolyte, when a non-uniform concentration field of ions exists. Ultimately, the obtained results helped to propose an explanation on a commonly observed phenomenon (melt rate-polarity relationship) during the DC operation of the electroslag remelting (ESR) process.


Introduction
Diverse application areas for electrometallurgical processes exist in liquid metal battery (LMB) [1], extraction and refining of metals [2], copper magneto-electrolysis [3], and aluminum smelting [4]. SiO 2 , CaO, or Al 2 O 3 in their molten form are used as source materials for metal extraction [5]. The calcium fluoride CaF 2 -based electrolytes, typically made of CaF 2 , Al 2 O 3 , CaO, SiO 2 , FeO, and MgO, are widely used for refining alloys with a so-called electroslag remelting (ESR) process. Removing sulfur and non-metallic inclusions from the alloy, and providing heat through Joule heating into the ESR process are key functionalities of the molten electrolyte. Several electrochemical (Faradaic) reactions take place at the electrolyte-metal interface aiming at obtaining an alloy which is as cleaned and chemically refined as possible. For instance, Faradaic reactions of alloying elements such as Ti, S, and Al were reported [2,6]. The melt rate of the ESR electrode with positive polarity (anodic) was observed in situ to be higher than that for the ESR electrode with negative polarity (cathodic) during DC operation of the process [7,8]. Apparently, the aforementioned phenomena are attributed to the ionic properties of the molten electrolyte. Thereby, a dominantly ionic mechanism for the conduction of the electric current was noticed [9].
It is well-known that the distribution of electric current density significantly impacts the electromagnetic field and consequently the hydrodynamic behavior of the molten electrolyte [1][2][3][4]. Thus, an accurate prediction of the electric field/current distribution is a crucial step toward modeling those processes.
Customarily, three classes of the current distribution in the electrolyte are considered: primary, secondary, and tertiary. The primary current distribution ignores the electrode kinetics and concentration-dependent effects. The electrical resistance in the electrolyte is assumed to obey Ohm's law [1][2][3][4]10]. Similarly, the secondary current distribution neglects composition variations in the electrolyte that allows us to use Ohm's law. However, the potential drop at the electrode-electrolyte interface as a consequence of electrochemical reactions, known as activation overpotential, is taken into account [11,12]. The tertiary current distribution considers effects of variations in the electrolyte composition as well as electrode kinetics on the behavior of the electrolyte. The Nernst-Planck equations are solved subject to the approximation of electroneutrality [13,14].
Traditionally, based on primary current distribution, an effective electrical conductivity is used as the diffusion coefficient for the electric potential to model the electric field assuming a uniform concentration of ions in the electrolyte solution. Additionally, it is assumed that Faradaic reactions are extremely fast so that electrode kinetics/activation overpotentials are negligible [1-4, 10, 15]. These are rather vague assumptions, and require further examination. That is why we decided to perform a numerical study.
A typical molten electrolyte for the ESR process composed of CaF2 (%wt. 94)-FeO (%wt. 6) is considered. Here, only the ferrous ion (Fe 2+ ) participates in Faradaic reactions at the metal-electrolyte interface. In order to conduct the numerical study, two modeling approaches were investigated. In the first approach, the electric potential field was calculated using Ohm's law, and an effective electrical conductivity for the uniform electrolyte solution taken into consideration. In the second approach, the electric potential field was calculated using the Poisson-Nernst-Planck (PNP) equations and Butler-Volmer formula to take into account the transport/ reaction of ions [16]. Unlike the tertiary current distribution, we solve the full set of PNP equations. Therefore, the bulk electro-neutrality is obtained as a result rather than being imposed. This paper discusses influences of the operational parameters such as applied voltage and electrode polarity on the distribution of electric current density in the electrolyte. The computed results for the aforementioned approaches are compared. The goal is to identify the role of the transport/reaction of ions in the electrical behavior (e.g., electric/current density field) of the molten electrolyte. Based on the modeling results, we put forward a possible explanation for the relation between the polarity (positive or negative) and melt rate of the electrode during DC operation of the ESR process.

Mathematical model
Recently, we developed a one-dimensional model including two planar, parallel electrodes which were separated by a completely dissociated CaF 2 -FeO electrolyte subjected to a DC electric field [17]. The model enabled us to calculate vital parameters of the electrochemical system including activation and concentrated overpotential. We successfully managed to validate our model against experiments [17]. In the present study, we extended the model considering a 2D axisymmetric domain involving two iron electrodes with different sizes, separated by a CaF 2 -FeO electrolyte. The system is also operated under DC voltage.

Model description/assumptions
All symbols used in this manuscript are described in BNomenclature.^In order to achieve the goal of this study (confronting two modeling approaches: Ohmic versus ionic) and to avoid extra complexity, the following assumptions for the model are made: (i) It is assumed that the electrolyte remains stagnant. Thus, the transport of ions by the advection of flow is neglected. Generally, the Ohmic approach cannot be applied for a stagnant electrolyte. In fact, the Ohmic approach is a particular case of the ionic transport approach. The Ohmic approach (Laplace equation) can be deduced from the ionic transport approach (PNP equations) under the conditions of equal diffusion coefficients for involving ions and in the absence of concentration gradients. As such, a rigorous stirring of the electrolyte solution by the flow is required to achieve a uniform concentration field (no concentration gradient) in the bulk of electrolyte. To apply the Ohmic approach for a motionless electrolyte, one can assume that concentration gradients remain negligible at sufficiently low-current density. Here, a uniform concentration filed ( ∇ ! c i ≈0: ) is considered to apply the Ohmic approach as reported in Table 1. (ii) The electrolyte is assumed to be fully dissociated into its component ions. Therefore, all molecules split into ions,  [18].
At anode ð Þ: At cathode ð Þ: In terms of mathematical modeling and to ensure the mass conservation, the ferrous ion is injected into the electrolyte at the anode (positive mass flux), and removed from the electrolyte at the cathode (negative mass flux). We ignore the neutral iron (Fe) in our calculations. In other words, the electrodeelectrolyte interface is the boundary of the computational domain. As such, no calculation is carried out in the bulk of electrode.
(iv) Complexation of metal ions and formation of oxyfluorides, which may occur in the molten CaF 2 -FeO electrolyte at elevated temperatures, are not taken into account [19]. In addition, we neglect any other chemical reactions such as generation/recombination which may occur in the bulk of electrolyte. (v) The well-known Stern model is used to describe the formation of the electric double layer (EDL) at the electrode-electrolyte interface [16]. The structural change at the EDL (anomalous electrical capacitance) is ignored [20][21][22]. (vi) Here, the standard model (PNP equations) for a dilute electrolyte solution which contains point-like ions is applied. Correspondingly, steric effects are ignored [23].
According to the Ohmic approach, we solve the Laplace equation (Table 1) in which the effective electrical conductivity for the uniform electrolyte solution is the diffusion coefficient. Based on the ionic approach, the electric potential is dependent on the concentrations of all involving ions (Ca 2+ , Fe 2+ , O 2− , and F − ) via the Poisson equation (Table 1). A total number of four conservation equations are solved in a coupled manner to calculate concentration fields for each ion. The concentrations of ions are used to determine the charge density that in turn appears as the source term in the Poisson equation (Table 1).
At the surface of electrodes (anode or cathode), the total flux including migration and diffusion for nonreacting ions (Ca 2+ , O 2− , and F − ) is zero. However, the total flux is related to the current density for the reacting ferrous ion through the Butler-Volmer formula. The computational domain also contains insulating walls where the flux of electric potential is zero. In addition, a value of zero is set to describe the total flux (diffusion + migration) for all involved ions (Ca 2+ , Fe 2+ , O 2− , and F − ) at insulated walls.
All governing equations, pertaining to boundary conditions and system parameters for the two mathematical modeling approaches, Ohmic and ionic, are listed in Table 1. They are described in detail elsewhere [17].
Once the electrolyte is subjected to the applied voltage, an electric double layer (EDL) forms at metal-electrolyte interface where the net charge density exists. The EDL comprises multilayers: the compact stern layer, diffuse charge layer, and diffusion layer [16]. Ions are adsorbed to the electrode's surface where they participate in Faradaic reactions in the extremely thin (~0.1-1 nm) Stern layer. The diffuse charge layer (~1-100 nm) contains the net charge density, which has a characteristic size known as the Debye screening length (λ D ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ).
Furthermore, the diffusion layer remains electrically neutral where the concentration of ions may differ from the bulk concentration. It is both logical and effective to formulate the PNP equations together with their boundary conditions in the dimensionless form [17,[24][25][26]. Subsequently, the Poisson equation (Table 1) in the dimensionless form is expressed as where . According to this formula, the value of ζ for a typical mesoscale process (L ref~1 cm) is extremely small (ζ~10 −7 ). As described by Bazant [24][25][26], it is impossible to satisfy all the boundary conditions assuming a negligible screening length compared to the domain size (ζ~0). Physically, this assumption renders the bulk electro neutrality condition (ρ ¼ ∑ i z i c i ∼0 ), that is always violated near electrodes within the diffuse charge layer. Therefore, it is necessary to use a rather large value for ζ (e.g., in this study ζ0 .0005) to emphasize the role of the diffuse charge layer and to achieve a numerical solution [24][25][26]. In other words, a quasi-electro-neutral bulk is assumed that allows for capturing asymptotic behavior of the electrochemical system.

Numerical solution
In this case, the numerical simulation was carried out for a medium-sized computational domain, including the electrolyte (1 cm), which was confined between two unequally sized electrodes (1 and 2 cm). Configurations of the system are shown in the left half of Fig. 1(a1, b1, c1), and corresponding boundary conditions are listed in Table 1. Principally, the model is not limited to the size of the system, and it can be effectively applied for largescale electrochemical systems (several centimeters). However, the total number of required mesh elements increases with the increased size of computational domain. Therefore, a large-scale system may contain several million mesh elements requiring an enormous amount of computational time. Here, a domain with 0.3 million volume elements was considered. A mesh of variable element size was generated considering a very fine mesh in the vicinity of electrodes (minimum cell size~2 μm) to resolve the boundary layer. The size of the elements was smoothly and incrementally increased away from the electrodes toward the bulk of electrolyte, considering a successive ratio of 1%. All equations were implemented into the commercial CFD software, FLUENT-ANSYS v.14.5. The numerical solution was carried out using the wellestablished finite volume method (FVM), which guarantees mass conservation of ions [27]. Furthermore, the model can be efficiently used to investigate complex multi-ion electrochemical systems as no restriction is required for the number of involved ions (reacting or nonreacting). The temporal numerical discretization is based on first-order implicit method. The spatial discretization is obtained by the third-order MUSCL scheme that enables us to accurately handle variable grid sizes [28]. Table 2 summarizes all required parameters for the calculations.
Transient calculations were made, but only the steady-state results were evaluated. In other words, the equations were solved by means of the pseudo-transient computation technique in which the transient terms (temporal derivatives) were retained [27]. In a mathematical point of view, the pseudotransient approach is a homotopy that embeds the steadystate problem in a space-time setting [33]. Capturing the transient behavior (temporal accuracy) of the system is not the objective. As such, the temporal accuracy is sacrificed in favor of rapid convergence to the steady state [27,33]. This technique is robust/effective to achieve the steady-state solution for complex problems involving highly non-linear equations when the initial iterate is far from the final solution.

Results and discussions
Several calculations were performed to explore the influence of modeling approach (Ohmic versus ionic) on the field structures (e.g., current density, concentration, etc.) in the electrolyte. For the ionic approach, the response of the system to the applied voltage (amount of imposed current) and polarity of electrodes are investigated. A summary of the results is illustrated in Fig. 1. Configurations of the system (left half) and the calculated electric field (right half) are shown in Fig.  1(a1, b1, and c1). The concentration fields for anions including O 2− (left half) and F − (right half) are illustrated in Fig.  1(a2, b2, and c2). The concentration fields for cations including Ca 2+ (left half) and Fe 2+ (right half) are illustrated in Fig.  1(a3, b3, and c3). In a nutshell, column 'A' made of (a1, a2, and a3) in Fig. 1 shows field structures at low voltage (1 V) with positive polarity for the small electrode, whereas column 'B' made of (b1, b2, and b3) shows results at low voltage (1 V) with negative polarity for the small electrode. Finally, column 'C' made of (c1, c2, and c3) indicates the response of system at high applied voltage (5 V) with the positive polarity for the small electrode. In all cases, the electric field is found to Column (a) made of (a1), (a2), and (a3) shows the aforementioned parameters for the system involving the small electrode with positive polarity at low voltage (V app = 1 V). Column (b) made of (b1), (b2), and (b3) shows the results for the system involving the small electrode with negative polarity at low voltage (V app = 1 V). Column (c) made of (c1), (c2), and (c3) shows the results for the system at high voltage (V app = 5 V) be remarkably strong near electrodes as a consequence of concentration overpotential [11]. The behavior of nonreacting ions including Ca 2+ , O 2− , F − follows the same trend.
Once the system is subjected to the applied voltage, the cation (Ca 2+ ) accumulates near the cathode, whereas the anions (O 2− and F − ) depart from the cathode and move toward the anode. In contrast, the operational parameters like polarity of the electrode and applied voltage can significantly influence the spatial variation of the concentration field for the reacting ferrous ion (Fe 2+ ). As shown in Fig. 1(a1), the ferrous ion is injected into the electrolyte at the anode. However, the concentration of Fe 2+ at the anodic electrolyte-electrode interface is low (Fig. 1(a3)) as a consequence of a strong, local electric field which impels the ferrous ion to migrate to the bulk. Eventually, a region of high concentration of ions appears near the anodic surface in the bulk, as well as at the cathode surface, for the system involving the small electrode with a positive polarity operating under a low voltage. Similarly, the concentration of Fe 2+ is small at the anode surface for the system involving the small electrode with negative polarity operating under low voltage as shown in Fig. 1(b3). However, the concentration of the Fe 2+ remarkably decreases due to high removal rate (current density) near the cathode. Figure 1(c3) shows concentration of Fe 2+ for the system involving the small electrode with positive polarity operating under high voltage. As a consequence of the high injection rate (current density) at the anode and vigorous electromigration in the bulk, an intense deposition of Fe 2+ at the surface of cathode is observed. The influence of electrode polarities on the electrical behavior of the electrolyte at constant applied voltage (e.g., here 1 V) is further elucidated. As shown in Fig. 2, variations of electric potential, charge density, and diffusion/migration flux of the reacting ferrous ion are plotted along the axis. The potential field is determined by the distribution of ions (reacting and non-reacting) through the Poisson equation (Eq. 1). As shown in Fig. 2a, the intense deposition of the cations/anions notably increases the potential drop near the surface of the cathode/anode electrode. Figure 2b indicates that the accumulation of charge density occurs in the vicinity of electrodes, whereas the bulk expectedly remains electroneutral.
As shown in Fig. 2c, the magnitude of diffusion flux for the ferrous ion, N diff ¼ D Fe 2þ ∇ c Fe 2þ k k , is plotted along the axis. Furthermore, Fig. 2d illustrates the calculated magnitude of the migration flux for the ferrous ion, N mig ¼ , across the electrolyte. A comparison is made between Fig. 2c, d to explore the relative strength of both flux parts (diffusion and migration). The gradient of concentration/electric potential and accordingly diffusion/ migration flux is larger near electrodes than that in the bulk. As a consequence of small gradient of the electric potential, the migration flux becomes very low in the bulk of electrolyte where the diffusion flux is dominant. Independent of the electrode size (large or small), the migration flux is larger than the diffusion flux near the anode. However, a relatively similar strength in flux parts (diffusion and migration) is obtained near the cathode. Figure 3 illustrates the calculated distribution of current density in the bulk of electrolyte for four cases. The highest amount of current density (J max ) is observed under the edge of small electrode, whereas a very low current density flows near the upper edge of the domain. Considering the Ohmic approach, Fig. 3a simultaneously shows the current density and electric fields. They are linearly dependent with the correlation called electrical conductivity (Table 1). An elliptical profile of the current density/electric field develops under the shadow of the electrode. It is notable that the ion concentrations which vary in the electrolyte must be ignored to be able to apply Ohm's law [16]. Consequently, the pattern of current density remains unchanged, no matter how operational  parameters such as voltage and electrode polarity change. In contrast, the current density field is significantly influenced by operational parameters, considering a non-uniform concentration field of ions as shown in Fig. 3b-d. At a steady state, the current density is determined by the flux of the reacting ferrous ion (Table 1), that in turn is dependent on the voltage and electrode polarity. As shown in Fig. 3b, c, either a semicircular profile or a dumbbell-shaped profile of the current density may develop by altering the polarity of the small electrode from positive (anodic) to negative (cathodic), respectively. The parabolic profile is obtained as the voltage increases, as shown in Fig. 3d. The obtained results reveal that Ohm's law is not valid in the bulk of the electrolyte with non-uniform ionic concentrations. Generally, the distribution of Joule heating (Q) significantly impacts the thermal field in the electrolyte. According to Table 1, the amount of released Joule heat to the electrolyte correlates with the current density. A comparison is made between Fig. 3b, c to explore the effect of electrode polarity on the Joule heat near the small electrode. Obviously, the amount of current density (~Joule heating) is much greater near the small electrode, with the positive polarity compared to that of the negative polarity. The aforementioned findings can well explain the in situ observation of higher melt rates for an anodic ESR electrode (positive polarity) than the cathodic one (negative polarity) in the DC-operated ESR process [7].
Our work demonstrates that the non-uniformity of ionic distribution can significantly impact an electrochemical system. Traditionally, the presence of excess of non-reacting ions (called supporting electrolyte) permitted a simplification in the analysis. Concentration gradients for the non-reacting ions had been thus far neglected, allowing for the assumption of a constant electrical conductivity (using Ohm's law) for the wellmixed electrolyte [16,34]. In the presence of flow, concentrations of ions near electrodes differ significantly from their bulk values. The flow may sufficiently promote stirring in the bulk of electrolyte, so that the concentration field becomes uniform. As such, it is necessary to include the advection effect of ions, caused by the flow, to improve our model. Nevertheless, the presented results can already explain the observed phenomenon (melt rate-polarity relationship) in the ESR process.

Summary
In certain industrial processes such as electroslag remelting (ESR), liquid metal battery (LMB), copper magneto-electrolysis, and aluminum smelting, the electric current is mainly conducted by the transport of ions in the high-temperature molten electrolyte. The hydrodynamics of those processes is dependent on the electric field that in turn is determined by the transport of ions in the molten electrolyte. Thus, a precise prediction of the electric/current density field is an essential step to model the molten electrolyte-metal processes. Here, a numerical study has been performed to compare two modeling approaches. Firstly, we applied the classical Ohmic approach to solve the Laplace equation, in which an effective electrical conductivity for the uniform electrolyte solution was considered. Secondly, we put forward an ionic approach, in which the Poisson-Nernst-Planck (PNP) equations were solved, for a quasielectro neutral bulk of the electrolyte. Previously, we proposed a one-dimensional model employing PNP equations to study the fully dissociated CaF 2 -FeO electrolyte at elevated temperature (1723 K). The model was successfully validated against experiments. Here, we refined the model to calculate a 2D axisymmetric domain including two unequally sized iron electrodes separated by the molten CaF 2 -FeO electrolyte. It was found that regardless of variations in any operational parameter of the system like the applied voltage or electrode polarity, the pattern of current density remains always unchanged if the Ohmic approach is used. In contrast, significant differences in the electrical behavior were observed in response to the operational parameters for the non-uniform concentration of the electrolyte solution according to the ionic transport approach. This implies that Ohm's law is invalid in the bulk of an electrolyte where the electric current is affected by the movement of ions. Finally, modeling results help to explain the dependency of melt rate and polarity for an ESR electrode on the operational parameters. This dependency is an observed phenomenon during DC operation of the ESR process.