Smoothed Particle Hydrodynamics Modeling of Electrodeposition and Dendritic Growth Under Migration- and Diffusion-Controlled Mass Transport

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.

is used in this work where α = 63/478πh² 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].

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 m₀ and the original solid particle’s mass becomes (m_i−m₀).

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²) of 16. The domain is divided between fluid and solid particles.

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.

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

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

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

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

There are no conflicts of interest.

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