Modeling and simulation of electrodeposition: effect of electrolyte current density and conductivity on electroplating thickness

Electroplating or electrodeposition is a process carried out in an electrochemical cell where a current is used to form a coating on a metal surface. Developing and optimizing conditions for electroplating is time consuming and modeling and simulation could be used to optimize the electrodeposition process. Electrolyte current density and conductivity are important parameters for an electrodeposition system as they dictate the overall efficiency of flow of ions in the electrolyte system and thus optimization of these parameters is necessary. In this manuscript we report the development of a mathematical model to predict the electrodeposition of copper on cobalt chrome alloy in an electrochemical cell with copper and cobalt chrome alloy as the electrodes and copper sulfate as the electrolyte. The developed model was validated using experiments. The coating thickness of the samples was characterized using scanning electron microscope (SEM) and a thickness gage. At 30 min the model predicted the copper thickness to be 11.7 μm while experimentally the coating thickness was found to be 9.445+/-1.79 (mean +/SD) using SEM and 12.375+/-1.36 (mean +/SD) using thickness gauge. When predicting effect of current density the model accurately predicts general trends however the model seems to vary from experimental values in regions where there is significant effect of the electrochemical double layer that the model does not account for. The model accurately predicts the trend of effect of electrolyte conductivity on coating formation. The model can thus be used as a starting point to predict effect of process parameters on electrodeposition thickness *Correspondence to: Anil Mahapatro, Department of Biomedical Engineering, Wichita State University, 1845 Fairmount Street, Wichita, KS 67260, USA, Tel: 316-978-5912; E-mail: anil.mahapatro@wichita.edu


Introduction
Simulation using mathematical modeling is a tool that is used to predict, evaluate or optimize the performance of a proposed / current system under study over time [1,2]. Simulation is performed under specific input conditions and output of the model is compared with that of the actual system [3]. Physical processes could be modeled to minimize experimental trials necessary for optimization of the physical process [4]. Simulation modeling helps designers and engineers to understand the ways and conditions, in which a part could fail, the loads it can withstand and helps to avoid recurring usage of physical prototypes to analyze designs for new and existing parts [5]. There are various types of simulation modeling such as stochastic [6], dynamic [7] and Multiphysics modelling [8]. A stochastic model is performed from a source of randomness as it is based on certain assumptions of the system under study [9]. Statistical modeling is a stochastic model [10]. Examples of this type of modeling are linear regression [11], multiple regression [12], etc. Monte-Carlo method is a type of stochastic modelling [13]. Dynamic modelling [14] represents time aspect of a system. Simulation modeling used for optimizing material handling in a plant would be an example under this category [15]. Multiphysics modeling incorporates principle of physics, chemistry biology and engineering in a mathematical statement that describes the phenomena or system under consideration [16].
Electroplating or electrodeposition is a process carried out in an electrochemical cell where a current is used to form a coating on a metal [17]. In the electrochemical cell the metal to be coated is the cathode and the anode can be one of the two: sacrificial anode (dissolvable anode) or permanent anode (inert anode) [18]. Electrolyte in the electrochemical cell acts as a medium for the movement of electrons and forms the electric circuit between the electrodes [19]. Oxidation occurs at the anode while reduction reaction occurs at the cathode resulting in electrodeposition [20]. Electroplating is used for various applications including corrosion protection [20]. For example, electroplating of Palladium is used to manufacture catalytic converters because it has the ability to absorb excess hydrogen. Fasteners are electro-plated to have a better corrosive resistance [21]. Majority of the electrical parts and components are used after electrodeposition process. Silver electroplating has been used on copper or brass to enhance its conductivity and also used in silicon solar cells to increase its operating efficiency by 0.4% [22]. Nickel plating, tin plating and various alloys are used for corrosion protection on nuts, bolts, housings, brackets, other metal parts and components [23]. Though expensive, gold electroplating provides not only corrosion, but also tarnish protection [24].
Developing and optimizing conditions for electroplating is time consuming and modeling and simulation could be used to optimize the electrodeposition process. Various examples exist in the literature Where T is the thermodynamic temperature, S is the entropy and n i the number of moles of species i. Considering the exchange of matter between V' and V'', which takes place without energy exchange with their surroundings, it is satisfied that Each individual term of the sum is positive when the transport of species i is not coupled to transport of other species. In this case, i dn′ is determined only by (  and species i moves towards the region in which its electrochemical potential is lower, that is,  and vice versa. For example, this takes place in the case of ionic species in diluted solutions. When the transport of different species is coupled, one or more terms in the sum could be negative, but the sum is always positive.
The transport of the species i is described in terms of either its velocity i V  or its flux density i i i . If the area of the surface between the two volume elements is dA and its orientation is given by the unit vector n (from V ′ to V ′′ ), the number of moles of species i crossing the surface in a time dt is When the difference between µ′′  and i µ′  is not very large, it can be assumed that the rate of change of the amount of species i, ( i dn ′ )⁄dt, is proportional to the difference µ′′ i µ′  or to the gradient of this potential normal to the surface, The velocity of species i in linear approximation is expressed as Where u i is its mobility Flux density takes the form Where R is the universal gas constant and the Einstein relation between mobility and diffusion coefficient, D i =u i RT, has been used.
At constant temperature and pressure, the gradient i µ ∇  is caused by the changes in composition and electrical potential ø so that Where F is the Faraday constant and z i the charge number of species i. Substituting Equation (7) in Equation (6), the Nernst-Planck flux equation is obtained, where f denotes the ratio F/RT. The terms in the righthand side of this equation represent the transport mechanisms of diffusion and migration, respectively. Diffusion is a consequence of the random thermal motion of the particles which makes the concentration of all species uniform. Migration causes the influence of the electric field, , on the random motion of the charged particles, and Equation (8) shows that the particles a component of their velocity along the direction of the electric field as a result of this influence.
Flux of species i due to bulk flow in a moving fluid in X-direction, Therefore, the flux expression for each species i can be written as Where the ionic mobility u i , is assumed to be related to the diffusion coefficient D i by the Nernst-Einstein equation [32] that demonstrates the benefits of modeling of electrodeposition to predict plating outcomes [25][26][27]. For example Obaid et al., modeled the electroplating of hexavalent chromium to optimize the electrode spacing and anode height to obtain uniform thickness of the coating [28]. Loss of coating uniformity was observed when the electrode spacing decreased while greater electrode distance increased the uniformity of the formed coating [28]. The simulation results proved that an ideal electrode separation distance was necessary to obtain a uniform coating. It was also proved that, larger size of anode as compared to cathode resulted in a non-uniform coating [28]. Hughes et al., demonstrated that the electrode kinetics played an important role in the electrodeposition process [29]. Electrode kinetics defined the deposition process using the rate determining step and current distributions. Their simulation results proved that variables like surface electrode potential and the ion concentration in the electrolyte also influenced the kinetics associated with the deposition process [29].
Optimization of electrodeposition of new systems would require optimization of process conditions relevant to that system. This would require significant changes to a model or developing a new model specific to that system. Electrolyte current density and conductivity are important parameters for an electrodeposition system as they dictate the overall efficiency of flow of ions in the electrolyte system and thus optimization of these parameters is necessary. In this manuscript we report the development of a mathematical model to predict the electro deposition of copper on cobalt chrome alloy in an electrochemical cell with copper and cobalt chrome alloy as the electrodes and copper sulfate as the electrolyte.

Governing equations
A mathematical model was developed to predict the electrodeposition of copper on cobalt chrome alloy and to evaluate the effect of electrolyte current density and conductivity on electroplating thickness. Electroplating occurs due to mass transport in the solution as a result of migration in electric field, diffusion in concentration gradient and convection in a flow field [30]. The transport equation which was applied in the diffusion layer was based on the flux equation of the ionic species in the solution. The general mass balance equation is given below Where N i -molar flux, Z i -charge number, U i -mobility, F -Faradays constant, C i -concentration, D i -Diffusivity, ∅ -potential in the electrolyte, V -velocity The flux equations explaining the behavior of electrochemical systems are related to the Second Law of Thermodynamics [31]. Consider two neighboring volume elements V' and V'' of a solution that have same temperature and pressure but different electrochemical potentials of their constituents. The difference between the electrochemical potential of species i, µ′  and µ′′  , in these volume elements implies that this species tends to move from one volume element to the other as there is no distribution equilibrium. This motion of species i from one volume element to its neighbor is generically called transport of species i. As the transport of the different species in solution takes place under thermal and mechanical equilibrium, the change in internal energy U of these two volume elements is More energy is required by the reactants when there is a larger barrier for the activation energy.
is a probabilistic feature of the energy barrier component which has to crossed.
Pre exponential factor A, also known as frequency factor, gives a number of times the attempt was made to overcome the energy barrier consider the electrode reaction: Where k f and k b are the forward and backward reaction rate constants respectively. This is a general redox reaction, where O represents oxidized state and R represents reduced state. For this reaction, the equilibrium state is governed by the Nernst equation, which relates the equilibrium potential of the electrode (E eq ) to the concentration of the reactants and products (O and R) Overvoltage: A given current will require a penalty that should be paid in terms of electrode potential-penalty called overvoltage [34] due to irreversibility.
E_eq is the expected electrode potential and E is the electrode potential Where a and b are constants and I is the current density. This is the Tafel equation.

Rates of forward and backward reaction
For the above reaction, the rate of the forward reaction is given by: Rate of the backward reaction is given by: Reaction rate and current are correlated. Reduction occurs at the cathode and oxidation at the anode.
Net reaction rate: The net reaction rate or net current is given by

Potential dependence of k f and k b
Both k f and k b are potential dependent functions The forward reaction, which is a reduction, is an electron accepting process. The rate of reaction increases when the electrode potential reaches higher negative because the electrode loses electrons more easily. The opposite happens in the backward reaction, i.e., oxidation reaction.
The material balance equation for each ionic species at every point within the diffusion layer is as follows Where R i is the production rate of species i due to homogeneous chemical reactions.
Under the assumption, variations in composition are negligible in the electrolyte and migration of ions gives the net contribution to current in the electrolyte equation.
The concentration gradients in the above equations are neglected in the Secondary Current deposition interface and the current density is obtained from Ohm's law Where σ l denotes the electrolyte conductivity.
Since the electrolyte composition is assumed to be constant, the material balances are unsolved for the Secondary current distribution, The electroneutrality condition is given by the following expression:

Electrode kinetics
Electrodeposition is a combination of oxidation and reduction reactions occurring at the electrodes making the electrode kinetics very significant for the modeling of electrodeposition.
Consider the reaction given below: K f and K b are the rate constants of the forward and backward reactions respectively.

Reaction rates:
The reaction rate, also known as rate of reaction or speed of reaction, for a reactant or product in a particular reaction is defined as how fast or slow a reaction occurs.
The rate of forward reaction can be calculated by: The rate of backward reaction can be calculated by: The net rate of reaction can be calculated by Substituting equations (15) and (16) in (17), we get Equilibrium: Equilibrium is point at which the net reaction rate is zero. From the above equations, we can obtain the equilibrium concentration ratio as follows Where K is the equilibrium constant (K) Rate constant varies with temperature, generally it increases with T. The rate constant (K) and temperature are related by: Where, E a =activation energy, R=gas constant and A=pre exponential factor Activation energy: Activation energy [33] is the barrier that has to be overcome by the reactants before they are converted to product.
Where n is the number of electrons transferred. This is the Butler -Volmer [35] formulation of electrode kinetics. There are two components, of current i.e., anodic and cathodic, and it is exponentially dependent on potential. The net reaction rate is given by: Substituting the values for k f and k b This formulation is called the Butler -Volmer formulation of electrode kinetics At Equilibrium: At equilibrium, i=0

Thickness of deposition
The deposition process is assumed to take place through the following simplified mechanism: The first step is rate determining step (RDS), and the second step is assumed to be at equilibrium, which gives the following relation for the local current density as a function of potential and copper concentration: where n denotes the over potential defined as where ∅ (s,0) denotes the electronic potential of the respective electrode.
Equation at the cathode is given by: where n denotes the normal vector to the boundary.

Equation at the anode is
The amount of deposition is understood by the Faraday's laws of electrolysis [36] which states that the amount of a material deposited on an electrode is proportional to the amount of electricity used. For reduction of one mole of a given metal ion (charge of n + ), n moles of electrons are used for reduction. The total cathodic charge for the coating, Q(C), is the product of the number of gram moles of the metal coated, m, and the number of electrons required for the reduction reaction, n, Avogadro's number, N a (number of atoms in a molecule), and the electrical charge per electron, Qe(C). Thus, charge required to reduce m moles of metal is given by At equilibrium: The electrode potential and oxidation and reduction concentrations make net reaction rate zero. Thus: It can also be written as: Therefore, calling in the Nernst equation: Upon differentiating the above equation with respect to E, we get: The terms on the left hand side of the reaction sum up to 1 and called symmetry factors (for one electron transfer in the reaction).

Reductive and oxidative symmetry factors
The reductive symmetry factor is associated with the forward reaction which is represented by α The oxidative symmetry factor therefore becomes (1-α) The term α is the measure of the symmetry of energy barrier. If the change in the potential is same on both sides of the barrier, then α=0.5=1-α. Any asymmetry in the change causes fractional values of α.
Standard rate constants: If, k f = k f°, then E=Eº Similarly: k f° and k b° are termed as standard rate constants If the concentrations of oxidation and reductions are same, and the potential is maintained at Eº to cause any current flow: From above equation, k f° = k b°Th e larger the value of Kº, the faster is the equilibrium. Reactions with small standard rate constants are slow. The standard rate constant is large for simple redox couples The Butler-Volmer theory ( The product of Avogadro's number, N a (the number of atoms in a mole), and the electrical charge per electron, Qe(C) gives the Faraday constant, F. The number of moles of the metal reduced by charge Q is: The total charge used in the deposition can be calculated by the product of the current, I (A), and the time of deposition, t (sec), where the deposition current is constant. When current varies with time during the deposition process, Q can be calculated by, The weight of the deposit, W(g), can be calculated by product of the number of moles of metal reduced and atomic weight, Mw, of the deposited metal is given by: The thickness of the deposition, d (cm), can be solved by: Where ⍴ is the density of the metal (g = cm 3 ) and A is the area of deposition (cm 2 ).

Boundary conditions
All boundaries (B1 and B4) were insulating:

Initial conditions
The initial conditions of the electrolyte were:  Figure 1 shows the schematic representation of electrodposition and the corresponding model geometry for simutaltion of the electrodeposition process. Finite element methods were used to discretize the multiphysics governing equations. A finite element solver software COMSOL was used to simulate and predict plating thickness under varied process conditions. The model was set as a two diemnsional (2D) time dependent model (Figure 1). The deposition at cathode and the dissolving of the anode was set to take place at 100% current yield. During the electroplating process, changes in the electrolytic density of the electrochemical cell occurs which would result in the change in the density at anode and cathode. These changes could induced free convection in the cell, however based on our assumption that the variation in composition is small, the free convection component was neglected.

Computer implementation and simulation
The process was assumed to be time dependent as the boundary of the cathode was moving as the deposition of the metal was taking place. The model is governed by mass conservation for the copper ions Cu 2+ and sulphate SO 4 -2 and the electroneutrality condition. The upper boundary represented anode, and cathode was placed at the bottom. The vertical walls were assumed to be insulated (Figure 1).

Experimental methods
The developed mathematical model was validated experimentally. The electroplating of Copper on Cobalt-Chrome was conducted in a electrolytic cell. The electrolytic solution was prepared with 100 gms. of CuSO 4 in 100 ml of distilled water and 7 ml of 1N of H 2 SO 4 . The solution was kept on the stirrer to mix homogenously. A cobalt chrome strip of 14 cm × 14 cm was taken and an immersion area 1.4 cm × 1.4 cm cm was used and copper strip of 14 cm × 14 cm was taken and an immersion area 1cm × 1cm was used. The copper strip was connected to the anode and the cobalt chrome was connected to the cathode were immersed in the electrolyte and then current was passed. Electroplating was conducted under varying electrolyte conductivities such as; 4.23 S/ m 2 ,1.9 S/m 2 , 0.93 S/m 2 and 0.54 S/m 2 and varying current densities such as; 2.52 × 102(A/m 2 ), 3.57 × 102(A/m 2 ) and 6.12 × 102(A/m 2 ). Post deposition the samples were air dried and characterized for coating thickness. The samples were characterized with Scanning electron microscope (SEM) and a thickness gauge. For thickness analysis using SEM across sectional view of the coated sample was taken and the thickness of the coating determined. Four replicates were carried out and the coating thickness values were reported as mean +/-standard deviation (SD). An analysis of variance, ANNOVA (two way/one way) was used to determine statistical significant differences at the 95% confidence level.  for the first 500 sec the coating hasn't started, as the time goes by we see a steady coating of Cu on Co-Cr. At the end of 3000 sec we see the predicted coating thickness on Co-Cr to be 0.019 mm. To validate the model an electrodeposition experiment was carried out with same parameters as the model prediction for duration of 30min (1800 secs). Electrolyte test parameters for this experiments included an electrolyte conductivity of 4.23 S/m 2 , current density of 3.57 × 102(A/m 2 ) and a 1.4 cm × 1.4 cm area of contact. Four replicates were carried out to evaluate the variance of test results. The samples were characterized with Scanning electron microscope (SEM) and a thickness gauge to measure the thickness of the coatings formed. Figure 3 shows a representative SEM image used for calculation of the coating thickness while the thickness gauge was directly used on the sample to measure the thickness. For each sample measurement readings were taken at four different points and the average reported. Figure 4 shows the comparison of modeled (simulated) thickness and experimentally obtained thickness values.

Results and Discussion
From figure 4 and table 1, it can be seen that at 30 min the model predicted the copper thickness to be 11.7 µm (0.0117 mm) while experimentally the average coating thickness was found to be 9.445+/-1.79 (mean +/-SD) using SEM and 12.375+/-1.36 (mean +/-SD) using thickness gauge. Figure 4 shows that there was lot of variation of thickness values when using SEM as compared to the thickness gauge. It can be seen (Figure 4) that the thickness gauge measurement values (12.375+/-1.36) were found to be closer to the predicted values (11.7) as compared to that of SEM thickness values (9.445+/-1.79). This variation could be attributed to the challenges associated with using the SEM for measuring coating thickness. In SEM a precise cross sectional view was required to measure the thickness accurately. Deviation from an accurate cross sectional view (90O) could result in variation of thickness values. Also if the layer was not uniform it could cause variation in the thickness values measured. Thickness gauge measurements were relatively consistent in its measurements. However the relative accuracy of the predicted and experimental values validates our model. This demonstrates the preliminary applicability of the model in predicting electro deposition of copper on cobalt chrome substrate.

Effect of current density on coating thickness
To evaluate the effect of current density on coating thickness the model was used to simulate and predict the electroplating thickness at three different current densities. This experiment was conducted with different current densities at an electrolyte conductivity of 4.23 S/m 2 on 1.  Figure 5 showed experimental values to be lower than simulated conditions at lower current densities 2.50 × 102(A/m 2 ) and 3.75 × 102(A/m 2 ) while experimental values were higher than simulated values at higher current density of 6.12 × 102(A/m 2 ). This can be explained by the formation of an electrochemical double layer surrounding the electrodes which plays a critical role in the diffusion of ions from and towards the electrodes. The phenomenon of electrochemical double layer and its effect on diffusion of ions has been well documented in the literature [37][38][39]. At lower current density there will be a higher resistance from the double layer resulting in fewer ions going through the double layer resulting in lower coating values. While at higher current density the ions overcome the resistance of the double layer which results in higher coating thickness. The developed simulated model currently does not account for the double layer     formation and its effect on coating thickness. This is something that will be needed to be incorporated to simulate more accurate predictions over a wide range on current densities. Comparing SEM and gauge thickness values we see similar higher difference with larger coating thickness as compared to lower thickness due to larger variances SEM measurements as discussed previously.

Effect of electrolyte conductivity on coating thickness
This experimental test was conducted with different electrical conductivities of the electrolyte keeping the other parameters constant (duration 46 mins and current density 3.57 × 102(A/m 2 ). The electrolyte conductivity evaluated were 4.23 S/m 2 , 1.9 S/m 2 , 0.93 S/m 2 and 0.54 S/m 2 on a 1.4 cm × 1.4 cm area of contact. The results from the experiment are tabulated in Table 3. Table 3 shows a decreasing trend in the simulated thickness coatings with decrease in the electrical conductivity of the electrolyte. Experimental determination of the coating under similar conditions resulted in a coating under one condition (4.23 S/m 2 ), while the coating either was not uniformly formed or not formed at all under other conditions (1.9, 0.93 and 0.54 S/m 2 ). Thus although the model accurately predicted the general trends, experimental observations indicated difficulties in forming a uniform coating at those electrolyte conductivities (1.9, 0.93 and 0.54 S/m 2 ). Thus both the model and experiments indicates the requirement of appropriate electrolyte conductivity for the electrodeposition to occur. Lack of appropriate conductivity would hinder flow of ions within the electrolyte solution leading to difficulties in formation of a uniform coating.

Summary and Conclusions
In summary we have developed a mathematical model to simulate the electrodeposition of copper on cobalt-chrome alloy. At 30 min the  Table 3. Electrodeposition coating at varying electrolyte conductivities NA: coating not uniformity formed model predicted the copper thickness to be 11.7 µm while experimentally the coating thickness was found to be 9.445+/-1.79 (mean +/-SD) using SEM and 12.375+/-1.36 (mean +/-SD) using thickness gauge. The relative accuracy of the predicted and experimental values validates our model. This demonstrates the preliminary applicability of the model in predicting electro deposition of copper on cobalt chrome substrate.
When predicting effect of current density the model accurately predicts general trends however the model seems to vary from experimental values in regions where there is significant effect of the electrochemical double layer that the model does not account for. The model accurately predicts the trend of effect of electrolyte conductivity on coating formation. The model can thus be used as a starting point to predict effect of process parameters on electrodeposition thickness however modifications are needed to incorporate experimental realities not accounted for in the model.