## Abstract

In many electrochemical processes, the transport of charged species is governed by the Nernst–Planck equation, which includes terms for both diffusion and electrochemical migration. In this work, a multi-physics, multi-species model based on the smoothed particle hydrodynamics (SPH) method is presented to model the Nernst–Planck equation in systems with electrodeposition. Electrodeposition occurs when ions are deposited onto an electrode. These deposits create complex boundary geometries, which can be challenging for numerical methods to resolve. SPH is a particularly effective numerical method for systems with moving and deforming boundaries due to its particle nature. This paper discusses the SPH implementation of the Nernst–Planck equations with electrodeposition and verifies the model with an analytical solution and a numerical integrator. A convergence study of migration and precipitation is presented to illustrate the model’s accuracy, along with comparisons of the deposition growth front to experimental results.

## Introduction

Electrochemical systems are ubiquitous and found in everything from commercial electronics to alternative fuel vehicles. In electrochemical processes, mass transport and surface reactions are driven by complex, coupled phenomena at the interfaces, such as electrodeposition [1]. Electrodeposition is the deposition of metallic ions (disassociated in an electrolyte) onto an electrically conductive surface (the electrode) in the presence of an electric field [2,3]. This process is used to plate materials such as automobile parts and beverage containers [4,5] or to create thin metallic films for devices such as fuel cells and water splitting [6,7], and also describes the process of charging lithium metal batteries (LMBs) [8–11].

One of the primary challenges associated with the electrodeposition process is non-uniform deposition and subsequent dendritic growth. Control over dendrite morphology is an important aspect of industrial processes [12]. Additionally, dendritic growth can cause safety and performance issues in LMBs [8,10,13]. Experimentally observing the electrodeposition process and dendrite morphology in situ at an interface is challenging due to its embedded nature and small scale. Alternatively, computational fluid dynamics (CFD) methods can isolate this region and study the critical physics and driving forces at the interface during electrodeposition [14,15].

Modeling geometrically and physically complex and moving reactive boundaries, i.e., non-uniform deposition or dendritic growth, can be challenging for CFD methods, especially in mesh-based methods. The interfacial region is difficult to continuously track and the mesh must be updated at each time-step leading to computationally inefficient models [16]. Moreover, the predictions of complex structures, such as dendrites, in mesh-based methods struggle to capture experimental data. These challenges make mesh-based CFD methods less than ideal for modeling electrodeposition.

To simplify modeling, some approaches assume that the concentration of charged species is constant at the interface [11,17]. However, fluctuations of ionic transport, which occur when the ionic concentration varies spatially and temporally [18], are the primary cause of non-uniform deposition [19]. Tracking the ionic concentration variation at the interface can enable greater understanding of the causes of non-uniform deposition. Recent computational modeling has focused on alternative methods such as Lagrangian methods, diffusion limited aggregation, and phase field methods (PFM). PFM has received the most attention in recent years; however, PFM has struggled to accurately capture the dendrite morphology and produces artificial symmetries and morphologies without modifications [20–22]. Chen and Pao [23] have made improvements to the PFM modeling of dendrite growth by using a rate modification factor. Also, it is challenging to model complex structures in the diffusion region surrounding the anode with PFM, such as separator microstructures or anode protective layers. PFM also requires an additional differential equation (the phase-field equation) to resolve the interface which increases computational time.

In this study, we present a mesh-free model of the electrode-electrolyte interface for the study of electrodeposition based on smoothed particle hydrodynamics (SPH). The SPH method provides inherent solutions to many of the challenges of electrodeposition modeling [24]. SPH is a mesh-free, Lagrangian particle method that uses an interpolation scheme to solve the governing partial differential equations (PDE). The PDEs can be solved explicitly and exactly which simplifies their implementation. SPH conserves mass and does not need to explicitly track interfaces so that complex geometric structures, such as dendrites, can be handled without undue computational resources [25].

The model presented in this work builds upon a previous SPH model of diffusion and reactive interfaces with precipitation [26–29]. As discussed in the following sections, two main additions to the model have been included that enable more robust modeling of the electrodeposition process: (1) calculation of the spatially and temporally varying electric field, which is dependent on the local ionic concentrations, and (2) electrochemical migration mass transport. These additional physics allow the model to represent the electrodeposition process more accurately at the electrode-electrolyte interface. The implementation of the model is verified through comparison to an analytical solution and a numerical integrator. A convergence study of ionic migration and deposition growth front is presented. Finally, a qualitative validation of dendrite morphology is presented in comparison to experimental results for three operating conditions: high reaction rate, low reaction rate, and pulse plating.

## Governing Equations

**r**

_{f}in the electrolyte, Ω

_{f}, is governed by the Nernst–Planck equation, which includes mass transport mechanisms by both diffusion and migration [18]

*C*is the concentration of species

_{i}*i*,

*D*is the diffusion coefficient of species

_{i}*i*,

*μ*is the migration mobility coefficient of species

_{i}*i*, and

*ϕ*is the potential.

*a*and

*c*for anions and cations respectively, are dissolved in the electrolyte solution and remain electrically neutral [1]. If an electric potential is applied to the solution, the cations migrate toward the negatively charged electrode and the anions migrate away from the negatively charged electrode. As a result, the electric potential is altered by the charge disparity in the solution. The electric potential distribution is governed by the electrostatic Poisson equation [1]

*ρ*is the electric charge density given by

_{e}*e*is the elementary charge,

*z*and

_{c}*z*are the cation and anion electric charge, and $\epsilon $ is the permittivity constant of the electrolyte. Note that in this form of Eq. (2), the permittivity constant of the electrolyte is assumed to be spatially independent.

_{a}^{+}) are reduced at the electrode/dendrite surface

*k*, and the ionic concentration at the electrode-electrolyte interface, Γ. The reduction reaction rate is mainly controlled by the operating conditions of the process, although secondary reactions at the surface can also play a role in regulating the reaction rate [30,31]. A general first-order reaction equation is implemented [1]

^{+}).

## Smoothed Particle Hydrodynamics Discretization

*a*}, and a resultant scalar field

*A*(

*r*) approximated by

*, and*

**r**_{i}, m_{i}*ρ*are the position, mass, and density of particle

_{i}*i*, respectively.

*W*is the SPH smoothing function, which is compact and non-zero up to the distance

*h*from particle

*i*. The distance between particle

*i*and

*j*is

**r**

_{ij}. Similarly, the gradient of the scalar field, $\u2207A$, can be approximated by

*M*

_{6}smoothing function [32]

is used in this work where *α* = 63/478*πh*^{2} for two spatial dimensions. *W*(|**r**_{ij}|, *h*) will be denoted as *W*_{ij} in the following equations.

The SPH model discretizes the computational domain into two sub-domains of particles: Ω_{f} is discretized with “fluid” particles that make up the electrolyte and Ω_{s} is discretized with “solid” particles that form the electrode and precipitates. Since there is no convection in the current model, both particle sets do not move but fluid particles can precipitate into solid particles and the solid particles can dissolve into fluid particles as described in Tartakovsky et al. [25].

*C*

_{c,}_{i}describes the cation concentration at position

*r**. The SPH-continuum surface reaction (CSR) formulation developed by Ryan et al. [33] reformulates the heterogeneous boundary condition of Eq. (4) as a homogeneous boundary condition and a volumetric source term,*

_{i}*S*. The volumetric source term is calculated from

_{v}*S*(Eq. (3)), a characteristic function (

_{s}*β*), and the surface normal vectors (

*n*) as

_{i}*β*is defined as

*a*refers to anion properties and

*C*

_{a,}_{i}describes the anion concentration at position

*r**.*

_{i}*i*and

*j*refer to particles

*i*and

*j*. The potential at the reactive surface (electrode and dendritic growth) is the reference potential, which is equal to the ground potential

*ϕ*

_{0}relative to the reference potential

The potential throughout the fluid domain is calculated using the discretized Poisson equation, Eq. (16).

The process of precipitation and dendrite growth is simulated by tracking the mass, *m _{i}*, in the solid. When the mass of a particle surpasses twice the initial mass,

*m*, the nearest fluid particle precipitates and becomes a solid. The new solid particle has a mass of

_{0}*m*

_{0}and the original solid particle’s mass becomes (

*m*−

_{i}*m*

_{0}).

The SPH model is implemented into the LAMMPS code base, as part of an SPH module [34]. The simulation domain (Fig. 1) is two-dimensional and includes 262,144 particles, discretizing a square domain of 128 by 128 units of h (except in the convergence cases) with an average particle density (number of particles within area *h*^{2}) of 16. The domain is divided between fluid and solid particles.

## Verification and Convergence Studies

While the implementation of the full Nernst–Planck equation has not been previously verified, the implementation of the diffusion and reaction terms has been verified [26,28]. As such, the verification cases presented here focus on verifying the electrochemical potential and migration equations for two oppositely charged species. Two separate verification cases are presented to assess the accuracy of the implemented mathematical models (i.e., Eqs. (11)–(19)); the first case considers the cation concentration change due to migration under a constant electric field. While the second case considers the concentration change for both cations and anions due to migration subject to a varying electric field. The electric field is concentration-dependent; it changes due to differences in local concentrations according to Eq. (16). For the second case, an analytical solution does not exist and so a comparison is made to numerical integration of the governing equations using the Runge–Kutta method. For both cases, the SPH simulations were conducted in a two-dimensional square domain (0 *< x < L*, 0 < *y* < *L*) and the one-dimensional concentration through the electrolyte was calculated by dividing *L* into *n* bins in the *x* dimension and then taking the mean of each bin.

*> x > L*= 1 and

*μ*E = 1 and the initial concentration is given as

*is 0.5 and (*

_{m}/L*σ/L)*

^{2}is 0.01 and the boundary conditions are

*C*

_{c}(

*x*= 0,

*t*) =

*C*

_{a}(

*x*= 0,

*t*) = 1 and

*C*

_{c}(

*x*=

*L*,

*t*) =

*C*

_{a}(

*x*=

*L*,

*t*) = 0. The initial concentration was selected so that the concentration varied smoothly; large disparities in concentration between the two species are non-physical.

*E*

_{L1}, is calculated as

*C*

_{A,i}and

*C*

_{B,i}are the concentrations predicted at location

*i*for methods

*A*and

*B*over all particles,

*n*. Equation (25) is used to calculate the L1 relative error between the SPH method, analytical solution, and Runge–Kutta numerical integrator.

The effects of particle ordering on the SPH results were also considered. In both cases 1 and 2, the SPH simulations were completed using particles placed with equal, ordered spacing. To consider the effects of disordered particles, case 2 was also run with randomly spaced particles. The particles are disordered by shifting them randomly up to a maximum of 20% of their initial spacing. As shown in Fig. 3, the ordered particle SPH simulation and the disordered particle SPH simulation compare well with a maximum difference between the ordered and disordered simulations of less than 0.11%.

Additionally, the effects of spatial resolution are compared for the SPH simulation of case 2 in terms of both the error and the computational cost. Figure 4 depicts the L1 error of the SPH simulations compared to the numerical integrator at different spatial resolutions as well as the time for the simulation to run ∼1.3 million time-steps. Based on the relative low error and fast run time, *L* = 128 *h* was chosen for the particle spacing.

A dendritic growth convergence study was conducted to demonstrate the model’s ability to accurately simulate dendritic growth. In Fig. 5, simulations are presented at differing resolutions (Fig. 5(a): *L* = 64 *h* and 65,536 particles, Fig. 5(b): *L* = 128 *h* and 262,144 particles).

The average growth fronts of the low- and high-resolution simulations were compared in Fig. 6. The average growth fronts were calculated using bin averaging in the *x* dimension. The simulations at the two resolutions predict similar average growth fronts with less than 5% relative difference (Fig. 7).

## Application to Dendritic Growth in Electrodeposition

As mentioned previously, non-uniform deposition can lead to dendritic growth, which is problematic in a variety of electrodeposition processes where precision is required. To further evaluate the SPH model’s ability to capture dendritic growth and morphology, the SPH model was used to simulate experimental data of dendrite growth under high and low current scenarios and pulsed plating conditions.

Experimental data from Schneider et al. [35] on electrodeposition of copper under galvanostatic conditions under various growth rates and pulsed plating conditions are used for comparison to the SPH model. The growth front lengths are plotted based on the experimental imaging (Fig. 8(b)). The growth front length is a normalized length, to measure surface roughness. It is normalized by the length for a smooth surface such that a smoother deposition will have a value close to one and dendritic growth will have a value higher than one.

The first two experimental tests are conducted at lower and higher growth rates. In the experimental work, the higher growth rate condition used an applied current that was five times higher than that of the lower growth rate condition. These were conducted to provide baseline evidence for different electrodeposition regimes where the lower growth rate from lower applied current produces reaction rate limited conditions and relatively uniform electrodeposition (Fig. 8(a), top), and the higher growth rate from higher applied current produces mass transport limited conditions and high dendritic growth (Fig. 8(b), bottom).

At lower growth rates (Fig. 8, top), the deposition is relatively uniform along the electrode surface as demonstrated by the normalized growth front length close to one. The growth front moves in unison suggesting the system is reaction rate limited and the local current density is also relatively uniform. At higher growth rates (Fig. 8, bottom), the applied current is increased 5× higher than the low growth rate conditions. The ionic deposition forms dendrites suggesting mass transport limited conditions. Under the high growth rate, the deposition morphology is a long, thin branching structure. There is rapid growth at the tips of the dendrites and hardly any deposition in the regions without dendrites. The normalized growth front length (Fig. 8(b), bottom) increases rapidly once there is mass transport limited condition for the reaction at the anode.

The SPH model is able to capture similar dendritic behavior. The reaction rate is the model parameter that is comparable to the applied current in the experimental tests. Initially, there is a 1 *h* thick layer of solid particles at *y* = 0 and the rest of the domain is fluid particles. The concentration at *y* = *L* is held constant throughout the simulations to represent the concentration outside of the diffusion layer, which is a common assumption in modeling electrochemical systems, and has been shown to accurately represent the physical system [1,36]. In Fig. 9, a simulation of dendritic growth (gray) is presented at high (*a*) and low (*b*) reaction rates. The high reaction rate is five times higher than the lower reaction rate, matching the conditions in the experimental tests. Long, thin dendritic branching occurs in the high reaction rate case and a more uniform deposition occurs at the lower reaction rate.

The normalized growth front length at different times is shown in Fig. 10. The curves for the low and high reaction rates simulations compare well with the trends of the growth front curves from the experimental tests. At the low reaction rate, the normalized growth front length remains between one and two indicating uniform deposition but at the higher reaction rate, the growth front length increases exponentially as expected when there are long dendrite growths and regions without growth.

Yang et al. also conducted experiments on the effects of pulse plating on dendrite growth (Fig. 11). Pulse plating has been suggested as a method for rapidly charging batteries while concurrently maintaining smoother deposition [37]. Yang et al. conducted experimental tests using the high growth rate of Fig. 8 for 1 s followed by a 5 s rest period, allowing the ionic concentration near the interface to equilibrate. The SPH model was also used to test pulse plating with a matching high reaction rate to rest period ratio. As demonstrated by both the experimental data and SPH simulations, pulse plating allows the deposition to occur more uniformly along the electrode surface resulting in a smaller normalized growth front length, as shown in Fig. 10. Note that the experiments by Yang et al. [37] focused on surface roughness and plating in electrochemical deposition on a copper plate and did not consider battery applications. In recent work, we have expanded these initial pulse plating studies to consider pulse plating effects for fast charging of batteries [38].

## Conclusion

Electrodeposition is used in a variety of industrial applications and is critical to the fabrication of electronics and the performance of batteries. Experimentally observing the electrodeposition process in situ at the electrolyte-electrode interface is challenging. However, computational models are able to isolate this region to understand the critical phenomena occurring at the interface during electrodeposition. With this understanding, the electrodeposition process can be better controlled leading to improved deposition in fabrication processes, and improved performance in battery applications.

A numerical model based on the SPH method was presented for simulating mass transport and electrodeposition near an interface. The Nernst–Planck equations were implemented in the model to track concentration changes of two oppositely charged species due to diffusion and ionic migration. The Poisson equation is solved to calculate the local potential throughout the domain. The expanded physics of the SPH model presented here allows more accurate modeling of the physical processes at the interface than previous SPH models that did not include electric field effects and is better able to simulate experimental conditions. The implementation of the Nernst–Planck equation was verified with an analytical solution and a numerical integration method. The SPH model was shown to accurately reproduce the verification cases.

Dendritic morphology predicted by the SPH model qualitatively matches experimental data for both low and high reaction rate scenarios. Additionally, the growth front length of the dendrites predicted by SPH agrees with experimental data and shows that for the high reaction rate the growth front increases exponentially, while the growth front length is stable for the low reaction rate. Pulse plating simulations were also presented as a viable technique for reducing dendritic growth. The SPH model accurately captures the behavior of ionic deposition for pulsed plating and predicts a growth front comparable to the experimental data. Further work on pulsed plating for fast charging applications, and current-voltage relations with the SPH model have been considered and are presented in Morey et al. [39] and Melsheimer et al. [38].

Explicitly tracking conditions, i.e., concentration and potential, near the electrode-electrolyte interface can lend insight into advanced methods for controlling the deposition morphology, such as structured electrolytes [40,41], novel separator designs [42–44], and surface treatments [39,45]. The results of the computational studies demonstrate the model’s ability to capture the complex interfacial physics occurring during the electrodeposition process. Additionally, the model could be used to understand the mechanical effects of dendrite growth in non-liquid electrolytes with the addition of solid mechanics equations, which have been used in SPH modeling of other multiphase systems for biological and geological applications [46–48].

## Acknowledgment

Financial support for this research is provided by the National Science Foundation through Award Nos. 1911698, 1727316, and 2034154.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.