Numerical Analysis of an Electroless Plating Problem in Gas–Liquid Two-Phase Flow

: Electroless plating in micro-channels is a rising technology in industry. In many electroless plating systems, hydrogen gas is generated during the process. A numerical simulation method is proposed and analyzed. At a micrometer scale, the motion of the gaseous phase must be addressed so that the plating works smoothly. Since the bubbles are generated randomly and everywhere, a volume-averaged, two-phase, two-velocity, one pressure-ﬂow model is applied. This ﬂuid system is coupled with a set of convection–diffusion equations for the chemicals subject to ﬂux boundary conditions for electron balance. The moving boundary due to plating is considered. The Galerkin-characteristic ﬁnite element method is used for temporal and spatial discretizations; the well-posedness of the numerical scheme is proved. Numerical studies in two dimensions are performed to validate the model against earlier one-dimensional models and a dedicated experiment that has been set up to visualize the distribution of bubbles.


Introduction
Electroless plating is an industrial chemical process aimed at forming a thin film or layer on a base substrate by reducing complex metal cations in a liquid solution [1][2][3]. This technique has been widely applied in various industries. For instance, surface decoration, hard-wearing coating, manufacture of hard-disc drive, printed circuit boards, etc. [3,4]. It is well known that many electroless plating processes produce parasite hydrogen bubbles, such as electroless nickel and electroless copper systems [5][6][7]. There are several works on the simulation of electroless processes that study the convection or migration of chemical species under a single-phase flow (e.g., [8,9]). As far as we know, there is no computational work on bubble generation in an electroless process. In large-scale electroless plating processes, the gas generation will not be a serious issue, in general, however, for microscale plating, hydrogen bubbles may prevent electrolytes from going into the region needed to be plated. The existence of relatively large bubbles has been an important issue in the study of micro-fluids [10][11][12]. As experiments are difficult, there is a need for a reliable numerical simulation tool.
From a theoretical physics point of view, electroless processes are complex; the entire physical system participates to the plating. Therefore, for the simulation, we chose a system that includes a gas-liquid two-phase flow, chemical species transport, surface reaction, and moving boundary due to deposition.

The Modeling
Because bubbles and sprays are frequent in engineering design with fluids, numerous papers on the modeling and simulation of gas-liquid two-phase flows have been published such as [13][14][15][16][17]. Working models to compute gas-fluid flows can be sorted into two Results show that the bubbles are not only appearing on the copper plate, but also near the top glass sheet. In the video, one can see that there were several bubbles going to the top from the center or the bottom side of the channel. The bubble size distribution was not measured and it ought to be done in the future and combined with numerical simulations as in [37]. The measurement of the bubble densities in the experiment is possible by means of acoustic, optical, and laser diffraction approaches [38]. However, these techniques are hard (or expensive) to be adapted on a microfluidic scale.
The computer simulations qualitatively arrive at the same conclusion. The experiment indicates that the clustering of bubbles happens on both the top side and the bottom side of the channel. Second, the numerical simulation predicts that most bubbles are generated at an early stage and near the inlet. The experiment shows that the bubble generation is more exuberant near the inlet. This observation coincides with that of Figures 2 and 3a. The region near the inlet at t = 20 is of the highest concentration of dissolving hydrogen gas. In addition, large bubbles were observed at the back end of the copper plates, which is also the case numerically as seen in Figure 3b.
The comparison is qualitative but sufficient to assert that the simulation software is a more powerful tool to prototype future industrial designs. Potentially, it seems more reliable than experiments and it gives detailed information on the free boundary and on the speeds and concentrations of the chemical, which is highly important for the design of commercial systems. Thus, we are confident in the future of the numerical method. The mathematical properties of robustness, stability, and convergence, verified here numerically, are a certification of the computer software for future use.

Modeling Equations for Liquid-Gas Flow
Let Ω(t) be the time-dependent physical domain which is a thin channel between a top and a bottom plate. The boundary of Ω consists of the inlet Γ in , the outlet Γ out , the solid wall Γ wall , and the reacting surface S(t) (see Figure 4). for the test problem in Section 5.2 is initially a rectangle of size 10 mm × 1 mm. We assume a fixed inflow velocity and given chemical concentrations from the left on Γ in , a solid wall on the top side with a no-slip condition for the velocity, and a traction-free outflow on Γ out . On the bottom side, S(t) is a free boundary and its motion is given by (23). However as the reaction site is active mostly for x ∈ (1.5 mm, 5.5 mm), we may block the chemical reactions for x < 1 mm to avoid a corner singularity at the entrance and also for x > 6 mm because experiments show that almost no plating occurs there. In the regions x ∈ (1.0 mm, 1.5 mm) ∪ (5.5 mm, 6.0 mm) the numerical simulations may not be accurate due to the singularity caused by the discontinuity in the boundary conditions (see Figure 3 for details).

Volume Averaging
We review the derivation proposed by Ni and Beckermann [39].
Let Ω 0 (x, t) be an small open set to be observed in Ω(t) and Ω k ⊂ Ω 0 the set occupied by phase k and bounded by the interface ∂Ω k . Assume that ∪ k Ω k = Ω 0 and Ω k ∩ Ω j = 0, k = j. Let n k be a outer normal to ∂Ω k and w k the normal velocity of ∂Ω k .
Let Ψ be a function of a slow variable x and a fast variable y due to the phase change. The volume average of Ψ in phase k is Ψ k (x, t) = 1 V 0 Ω 0 (x,t) χ k (y)Ψ(x, y)dy, where χ k is the indicator function of the domain of phase k and V 0 = Ω 0 dx, assumed constant. The intrinsic volume average is defined as The volume fraction r k = V k V 0 has the properties ∑ k r k = 1 and Ψ k = r k Ψ (k) k . Some useful averaging formulas are listed below [40,41]: In principle, one should also introduce fast and slow time variables but it is assumed that spatially averaged functions are no longer varying fast in time.

Mass Conservation
We consider a gas and a liquid phase. Let ρ g be the density of gas, ρ l the density of liquid. We have the mass conservation for both phases (l for liquid and g for gas): whereṠ g is the mass gained owing to the precipitation of dissolved gas,Ṡ l is the mass loss when liquid evaporates into gas; u g (x, t), u l (x, t) are the volume averaged fluid flow of gas and liquid, respectively. Since the mass gained in gas balances the mass lost in liquid, we haveṠ For chemical species, we assume that the ions are transported only by the liquid electrolyte. Let c s be the volume averaged concentration of metallic ions destined to be deposited on the reacting surface, c g the volume averaged concentration of dissolved gas, and c k , k = k 1 , . . . , k M the volume averaged concentration of other chemical species participating to the chemical reactions. The equations for the concentrations are ∂ t (r l ρ l c j ) + ∇ · (r l ρ l c j u l ) − ∇ · (r l ρ l D j ∇c j ) − G j = 0, j = s, g, k.
where G j , j = s, k are interfacial terms due to the phase change. By (3), we can rewrite the above equation as where D j are the diffusion coefficients. In particular, since the gas is consumed by the phase change, by assuming that the gas precipitation is linearly dependent on the dissolving gas concentration [42,43], we have In the above, w l is the interface velocity of ∂Ω l and K is a constant independent of r g , r l ; c sat is the saturated concentration of the gas, M g is the reciprocal of the molar mass of the gas,Ṡ Moreover, G j , j = s, k, g can be estimated by (see Appendix A) For incompressible fluids, a volume conservation is derived from (3): By (4), the above reduces to ∇ · (r g u g + r l u l ) =Ṡ g 1 Recall that the physical domain is occupied either by gas or liquid, therefore r g (t) + r l (t) = 1 at all times.

Equations of Motion
Let µ g , µ l be the viscosities of gas and liquid. The volume averaged Navier-Stokes equations are used for momentum balance (see [39]): where p j , M D,j , F j , j = l, g are pressure, drag and external force terms. Following [24,39] where C D is a drag coefficient, and the external force are in fact the interfacial terms · n j dA, j = l, g, and D(v) = ∇v + (∇v) T is twice the symmetric gradient of v; F g and F l can be estimated by (see Appendix A) In view of (11), following [44,45], we assume a constitutive relation p = p l = p g in order to close the system of equations. The velocity fields of both phases are assumed to be 0 outside their own single phase region, respectively. Consequently, and by (3) and (12), the momentum equations simplify to r j ρ j (∂ t u j + (u j · ∇)u j ) + r j ∇p − µ j ∇ · (r j D(u j )) + γ j C D r g |u g − u l |(u g − u l ) = 0, j = g, l, with γ g = 1 and γ l = −1 (15)

Boundary Conditions
We consider a fluid flow from an input boundary Γ in to an output boundary Γ out with a solid wall at the bottom, Γ wall : The boundary conditions for r j , j = g, l are where is a fixed positive small constant. The boundary conditions for the concentrations of chemicals are, with c j,in given: Referring to Figure 5, if S(t) is the reaction surface, we denote S l (t) ⊂ S(t) the region occupied by the liquid and S g (t) := S(t) \ S l (t) the region occupied by gas. Choos-ing an arbitrary subset W ⊂ S(t), the surface reaction takes place only on W ∩ S l (t). Assuming that the concentration profile is uniform near the small region W, we have: Therefore by dividing both sides by W 1dA: for a positive number β indicating the chemical equivalence for gaseous molecular generation; F is the Faraday constant, and z is the atomic number of the material. I j is the current density satisfying the Butler-Volmer equation where R is the gas constant, E j are the chemical potentials of species j, θ is the temperature, α j , β j ,L j , κ j are constants, and E mix is given by writing electrical neutrality : On S(t), the fluid velocity induced by the deposition is where V s is a constant. Finally, S(t) moves according tȯ

Single-Phase Flow
If there is no gaseous phase in the system and no dissolved gas in liquid, i.e., r g = c g = 0, then u = u l andṠ g = 0 and mass conservation reduces to ∇ · u = 0. The convection-diffusion of chemicals become, and the fluid system reduce to the Navier-Stokes equations:

Notations
If f ∈ R, we denote by f + := max( f , 0) and by f − := − min( f , 0). We denote by · L p := · L p (Ω(t)) the L p norm on Ω(t), · W k,p := · W k,p (Ω(t)) the W k,p norm, and · H k = · W k,2 , 0 ≤ k ≤ +∞, 1 ≤ p ≤ +∞. Remembering that c k (x, t) is a vector, let us denote C = (c s , c T k , c g ) T . We assume the densities ρ j constant and denote α j = r j ρ j the mass fractions and ν j = µ j /ρ j the kinematic viscosities. The system is In addition asṠ g = −Ṡ l = Kα l (c g − c sat ) + /ρ l , we may use the redundant Equation (11):

Semi-Discrete Schemes
Let T be the final time and δt a time step. We denote by φ m , m = 0, 1, . . . , N := T/δt, the numerical solution of any physical quantity φ at time mδt. Convection terms are approximated in time by the method of characteristics. Let X m j (x) ≈ x − u m j (x)δt. Then Consider the following scheme For electroless plating the domain is Because of the asymmetrical treatment of α g and α l the scheme (30) does not imply that However, by (30), (31) and (33), we have By a Taylor expansion at x, we obtain Plugging (37) into (36), we have So the scheme is consistent with the equation for α g .

Positivity
Positivity of α m+1 l holds only if δt is small enough. When positivity is required absolutely, an O(δt) modification of (30) forces the positivity of α l : Indeed assume that α m l is strictly positive, or more precisely that α m l ≥ > 0 for all x; then we have

1.
Let us show first that (30) generates a bounded sequence {α m l } m=1..N . For clarity we assume homogeneous data at the boundaries. With simplified notations A multiplication by α m+1 and an integration on Ω m+1 leads to By the Cauchy Schwarz inequality and the positivity of φ m , The inverse of the determinant of the Jacobian of the transformation , where C is a generic constant function of ∇ 2 u m L ∞ , the norm of the Hessian of u.

2.
Stability of the scheme for C is shown by the same argument. 3.
Stability of the scheme for u g and u l is a consequence of a similar argument combined with the Ladhyzenskaya-Babuska-Brezzi saddle point theory (LBB) [46]. We denote by (·, ·) the L 2 inner product. For tensor-valued functions such that With self explanatory notations, the equations for the velocities (32) and (33) are written in variational form as: Find u g , u l and p satisfying the Dirichlet conditions and such that, where, for j = g, l, α j := α m+1 j , β j := 1 . (42) and where !g = l, !l = g. Here H 1 0 is the subspace of H 1 of functions which are zero on the Dirichlet boundaries. Note that the above is a semi-linearization of (32) and (33). However in Algorithm 1 below, the nonlinear problem is solved by an iterative fixed point which uses (41) and (42). The LBB theorem says that the solution of (41) exists and is unique because, for every Let us show stability in the special case f = 0 because one can always subtract w from α g ρ g u g + α l ρ l u l so as to work with u g,in = u l,in = 0 and f = 0. Thus, setting v g = u g , v l = u l and q = p leads to By the same argument used above, it implies that u j , j = g, l is bounded. Indeed, assuming α m ≥ 0, for some generic constant C depending on c m g L ∞ and ∇ 2 u m l L ∞ . Finally, we obtain This estimate would be optimal if the constant C did not depend on the Hessian of the velocities. This is a drawback of the characteristic method and of the unsophisticated treatment of the nonlinearity. At the expense of long mathematical arguments it could be fixed as in [34], the scheme would be H 1 stable.

4.
Note that we have swept under the rug the fact that at level m the domain of definition of the functions is Ω m and at level m + 1 it is Ω m+1 . The problem can be solved but at the cost of difficult notations and iterations between y m+1 and u m+1 ; for details see [47].

Finite Element Implementation
For simplicity, the physical domain Ω(t) is assumed to be a two-dimensional polygonal domain.

Mesh
Let {K h (t)} h>0 be an affine, shape regular (in the sense of Ciarlet [48]) family of mesh conforming to Ω(t). The conforming Lagrange finite element space of degree p on Ω(t) is where P p is the space of polynomials of degree p of R 2 .
We assume that the connectivity of the mesh K h (t) never changes with time.

Spatial Discretization
We use the Hood-Taylor element: the velocities are in V h (t) := (X 2 h,t ) 2 and the pressure is in P h (t) := X 1 h,t . For the volume fractions and the concentrations we use also P h (t).
Recall that the nodes of X 2 h,t are the vertices and the middle of the edges. Denote by {ψ 1 , . . . , ψ N a } the nodal Lagrange basis of X 2 h,t . As before On the boundaries where Dirichlet conditions are set, the functions are known. We denote P m 0h and V m 0h , the corresponding spaces where basis functions attached to a Dirichlet node are removed.

Mass Fractions
h satisfying the Dirichlet boundary conditions and such that

Remark 2.
A modification similar to (39) will insure the positivity of α m+1 l,h .

Concentration Profiles
subject to ∑ j=s,k where w g,h is the last component of w h and (20).
Let L g , L l , f be defined by (42). Data: Set u j = u m j,h , j = g, l. for n = 1 . . . N do Find u g , u l and p sastifying the Dirichlet conditions and such that, end Set u n+1 j = u j , j = g, l.

Consistence and Stability
Variational formulations discretized by finite element methods inherit the stability and consistency of the continuous equations. The LBB theorem applies also to the Hood-Taylor element for velocity pressure problems. Therefore, as in the continuous case, the H 1 norms of α m+1 j , u m+1 j , C m+1 j are less than (1 + C( ∇ 2 u m l L ∞ )δt) times the H 1 norms of α m j , u m j , C m j . If we could show that C( ) is bounded, then it would imply that the scheme converges when δt → 0.
Problem (53) can be formally expressed as a system of linear equations: where Φ is a (4N a + N q ) × (4N a + N q ) matrix, U m+1 and F m are (4N a + N q ) vectors. Note that Φ has the form with, for i, j = 1, . . . , 2N a , n = j = 1, . . . , N q , k = g, l,

Proposition 1. The linear system (56) is uniquely solvable.
Proof. According to the LBB theorem [46] the saddle point problem (52) is well posed when for p ∈ P m+1 Therefore Φ has full rank and is non singular. Q.E.D.

Iterative Process
At each time step, (47), is solved first, then (48) and (49) is solved iteratively by using a semi-linearization of the nonlinear boundary terms. Then (50) and (51) is solved iteratively by a semi-linearization of the nonlinear terms; each block involves the solution of a well posed symmetric linear system. Finally S m is updated by (23). Algorithm 2 summarizes the procedure.

One-Dimensional Electroless Nickel Plating Problem
Here we reproduce, with a two-dimensional computation, the one-dimensional study by Kim and Sohn [8]. In their work, the electroless nickel plating process on a rotating disk with constant angular velocity is considered. In this situation, the velocity field near the surface of the rotating disk can be approximated by a uniformly distributed flow towards the plating surface. In addition, the thickness of the diffusion layer is assumed uniform on the surface. Consequently, for the modeling, the domain becomes one-dimensional. Given that the gas generation is not considered and only the steady state is computed in [8], a single-phase recovery r l = 1, c g = 0 is applied. Finally, four partial reactions in the electroless nickel plating process are considered: Update the mesh by y m+1 = y m + δtu 2 m+1 g,h . A Gauss-Seidel smoother is applied if oscillations occur: y m+1 (x) is averaged with its neighbors, y m+1 (x − dx) and y m+1 (x + dx). end All chemical species are labeled as follows: c 1 is the concentration of the anodic hypophosphite (H 2 PO − 2 ), c 2 the concentration of the cathodic hypophosphite, c 3 the concentration of the nickel ion (Ni 2+ ), and c 4 the concentration of the hydrogen ion (H + ). Now the two-dimensional analog can be formulated: Let Ω = (0, δ 3 ) × (0, ), where δ 3 is the thickness of the diffusion layer for nickel and << δ 3 is a small positive number. The thickness of the diffusion layer for species j is given in [32]: The governing equation for the concentration profile is given by subject to the boundary conditions with the electron balance constraint The velocity field can be expressed as in [32]: where a = 0.51023 is an experimental constant, r = 0.995 is the ratio between the hypophosphite anodic part and the cathodic part on the reacting surface. The equilibrium potential E 0j for species j can be approximated by the Nernst equation, with pH = log(c 04 ): (68) By simulating system (64), (65), and (66), with the physical constants given in Tables 1 and 2, until a steady state is reached, the numerical tests show that the present model agrees well with the previous 1D studies of Kim a Estimated from [49]. b Assumed in this study. c Taken from [50]. d Calculated based on [51]. All values except   Figure 6. In red, the mixed potential E mix computed by the one dimensional system (65) versus pH= log c 04 . In black, the same but computed with the full two dimensional system.  Regarding (59)-(62), atomic nickel and phosphorus are deposited on the surface during the electroless process. The deposition thickness can be estimated in terms of the current densities: where V P , V Ni are molar volumes of phosphorus and nickel, respectively, and t is the deposition time.

Two Species in a Gas-Liquid Two-Phase Flow
Let the initial domain Ω be a rectange of size 0.01 × 0.001 (in meters). We consider complexed (by tartrate, denoted by L) copper ions, formaldehyde, and hydrogen dissolved in water, which are denoted by the subscriptions s, k, g, respectively, for the chemical species transport equations.
The chemical reaction can be expressed as the following two partial reactions: Given the above equations, we also use the subscriptions s and k to represent the quantities corresponding to (70) and (71), respectively.
The values of the physical constants are listed in Table 3.
For convenience, the following scalings are applied: The initial conditions are set to: constant phase ratio and Poiseuille flow: plus the first equation in (76) subject to (21). The inflow values are The boundary conditions on S(t m+1 ) are where α = 0.0005 and Boundary conditions on Γ out and Γ wall are as in Section 2.4. See also Figure 4.

Remark 3.
We note that α = 0.0005 is much larger than the experimental values; the numerical simulations produce u g 2 (and u l 2 ) of magnitude in the order O(10 −4 ). On the other hand, the deposition rate in a typical experiment is of order 1 µm per hours [52], which is not larger than O(10 −6 ). Yet the numerical test is conducted to validate the numerical method when the evolution of the domain is larger than real-life values.

Convergence
First, we conduct the convergence test for different time step with a fixed mesh. To obtain a "reference solution", the system (47)-(51) is solved with a 50 × 10 uniform mesh and a small time step δt = 0.01 and T = 10. The convergence with respect to δt is studied without changing the mesh; results are given in Table 4 and the rate of convergence for each variable is presented in Figure 8. Numerical tests for solving two-phase flow problem and volume fraction problem present a linear decay of L 2 error with respect to the time step. However, the convergence for solving the concentration profiles does not reach the expectation due to extremely low diffusion coefficients and large current densities. Table 4. L 2 error with respect to the reference solution provided with time step δt = 0.01, 50 × 10 uniform mesh, and α = 0.0005 for numerical simulation in Section 5.2 at T = 10. Second, we conduct the convergence tests for different time steps and mesh pairs; the time step is always proportional to the mesh size. The reference solution is obtained with 200 × 20 uniform mesh and δt = 0.05 at T = 10. Figure 9 and Table 5 present a linear decay of L 2 error with respect to the time step for each variable. Figure 9. Convergence with respect to δt and mesh size for the case described in Section 5.2 (see Table 5 for the time step and mesh size pair): log-log plot of the error for each unknown (note that the curves for u l and u l overlap). Third, a convergence test similar to the second test above is performed with different time steps and mesh pairs, but the errors are computed at T = 120. The reference solution is obtained with 200 × 20 uniform mesh and δt = 0.3. Figure 10 and Table 6 show a nearly linear decay of L 2 error with respect to the time step for each variable. In this test, r g and c k vary significantly as the grid is refined. The intensity maps of r g are given in Figure 11. The intensity maps of all variables present no significant differences with the reference solutions when δt ≤ 0.6. Figure 10. Convergence with respect to δt and mesh size for the case described in Section 5.2 (see Table 5 for the time step and mesh size pair): log-log plot of the error for each unknown (note that the curves for u g and c k overlap and the curve of u l is closed to them).

Robustness for Large Time Steps
With a large time step δt = 1 and 100 × 10 uniform mesh, the product of the maximal liquid fluid speed with the time step is around 1.5 times of the mesh size, which is optimal for the Galerkin-characteristic method. Solutions are displayed in Figures 12-16. (a) Intensity map of r g at t = 120 with δt = 1.2 and 50 × 5 uniform mesh.
(b) Intensity map of r g at t = 120 with δt = 0.6 and 100 × 10 uniform mesh.
(c) Intensity map of r g at t = 120 with δt = 0.4 and 150 × 15 uniform mesh.
(d) Intensity map of r g at t = 120 with δt = 0.3 and 200 × 20 uniform mesh.

CPU Time
With δt = 1 and 100 × 10 uniform mesh, it took 5832 s to reach the final time T = 180 with an Intel Core i7-8750H @ 2.20 GHz. During the computation, it took 0.086% of the total CPU to solve the volume fraction problem, 7.66% to solve the chemical species transport problem, and 91.93% to solve the two-velocities/pressure flow problem.
The computer program is written using the FreeFEM++ toolkit [35].

Results
In Figure 12a,b the velocity vector fields u g and u l are seen to be almost parabolic in y (Poiseuille flow), but the phase change and moving boundary induce a non-zero asymmetric vertical component u 2g (see Figure 15); both play important roles for the bubble distribution. Bubble density can be inferred by analyzing c g and r g (see Figure 3). The color maps of Figure 13 display a high gas volume fraction area near the top and bottom plates. Figure 14 shows how the steady state is established and how the electrolyte disappears in the plating region due to the plating. Figure 2 displays a high volume fraction of the gaseous phase near the reacting surface. The deposition-induced movement of S is presented in Figure 17. Figure 3b shows that the region of highest bubble density is moving away from the inlet as the electroless plating proceeds. Figure 17. The thickness of the deposition is given by the motion of S(t), plotted here at 5 instants of time, with respect to the x-axis (in mm). Notice that the motion t → S(t) is very small; the oscillations are blown-out of proportion by the scaling used in the graphic.

Comparison with Experimental Results
To validate the numerical method on a real-life problem, an experiment for reproducing the numerical study in Section 5.2 is conducted. Here, we shall show that the experimental result can be qualitatively fitted by the numerical simulation.
The experimental setting is described as the following: A micro-channel is enclosed by two sheet glasses of size 8 mm × 8 mm and another two of size 8 mm × 1 mm, which form a rectangular channel. The electrolyte goes in the channel from the left and exit on the right. One piece of the square sheet glasses is partially glued on a copper plate of size 8 mm × 4 mm, where the longer side of the copper plate coincides with an edge of the inlet (see Figure 1 for geometry setting). The inflow is set to be of average velocity 0.115 mm/s. At inlet, the copper ion concentration is c s0 = 39.34 mol/m 3 and the formaldehyde concentration is c k0 = 77.5883 mol/m 3 . Here, the inlet concentrations c g0 and c k0 are the reference concentrations for copper ion and formaldehyde, respectively. We further define the reference concentration of the hydrogen gas to be c g0 = 1 mol/m 3 . Other physical parameters are given in Table 3. Some parameters, for example, reference current densities i s , i k , and i g , may not be exactly same as what are given in Table 3. Nevertheless, they are acceptably closed to reality, or at least in the same order.

Experimental
To fabricate the test vehicle, a 4 inch glass wafer was first sputtered with 30 nm chromium and 200 nm copper which served as an adhesion layer and seed layer, respectively. The wafer was then diced into each 8 mm × 8 mm glass dies. To ensure a significant comparison between the regions being plated or not, each test die was half immersed in SPS (Na 2 S 2 O 8 ) solution and hydrochloric acid to remove the copper and chromium layer. The glass die turned out half transparent and half coated with copper where the electroless copper plating took place. Thereafter, a fully transparent glass, which was identical to the size of the test die, was face-to-face aligned and bonded via using a flip-chip die-bonder in order to obtain a clear observation view. Two tungsten wires which were 8 mm in length and 2 mm in diameter were glued by UV gel and placed on the periphery of the test die for the purpose of restricting the flow direction and defining the height between the dies (see Figure 18). The test vehicle was then subjected to a micro-fluidic system composed of a PDMS mode containing a micro-fluidic channel and a bottom glass. Clips were used to seal the micro-fluidic system and prevent the leakage of electrolyte. A peristatic pump was used to control the flow and connect the micro-fluidic system with a silicone tube. Prior to the electroless plating, the test vehicle was immersed in 10% sulfuric acid to remove copper oxide. Finally, the electroless copper plating was conducted in a water tank controlled at 50 • C with in situ recording via stereomicroscope (charged coupled device digital camera CCD). The electrolyte PHE-1 Uyemura possessing the given reference concentrations c s0 of (complexed) copper ion and c k0 of formaldehyde was used for the experiment. The complete equipment setup is described in Figure 19.

Results
Experimental results (see Figure 20) show that the bubbles are not only appearing on the copper plate, but also appearing on the top. In the video, one can see that there were several bubbles going to the top from the center or the bottom side of the channel. The region above the glass becomes darker with time. The simulation results (see Figure 13) qualitatively arrive at the same conclusion. The experiment indicates that the clustering of bubbles happens on both the top side and the bottom side of the channel. Second, the numerical simulation predicts that most bubbles are generated at an early stage and near the inlet. The experiment shows that the bubble generation is more exuberant near the inlet in comparison with other regions at t = 20. This observation coincides with that of Figures 2 and 3a. The region near the inlet at t = 20 is of the highest concentration of dissolving hydrogen gas. In addition, large bubbles were observed at the back end of the copper plates (i.e., region near (x,y) = (6,0) corresponding to Figure 4), which is also the case in Figure 3b.

Discussion
For an electroless plating process accompanying gas generation, the bubble distribution with respect to time, in the micro-channel, is the most important index for evaluating the quality of deposition. To measure it quantitatively, a high-quality optical system installed in the micro-channel is indispensable. For example, several types of fiber optical probes have been used to measure the particle (or bubble) size and distribution in a channel flow (or micro-channel flow) [53][54][55][56]. However, such an optical system is difficult to install in our case because there is no appropriate place to set up the light source and the detector in the micro-channel. The signal interference caused by the copper plate or glued gel on two sides is almost inevitable.

Conclusions
The numerical simulation of electroless plating is difficult for two reasons: multi-phase modeling and nonlinearities. We have proposed a phase averaged liquid-gas two-fluidvelocity/one-pressure system combined with phase densities and chemical concentration equations. The nonlinearities being similar to those of the Navier-Stokes equations, we have used a semi-Eulerian time discretization leading to a generalized Stokes operator for the two-velocity/one-pressure system; the inf-sup saddle point theorem has lead to a proof of stability and well-posedness of the discretized system by the Hood-Taylor finite element method. The two-phase flow model is compatible with single phase models when the volume fraction of gas and the concentration of the gas in the liquid phase are set to zero. The model is also compatible with the one-dimensional model proposed in [8]. The numerical results confirm the robustness of the method. To validate the model, a real-life experiment has been performed. The numerical results agree qualitatively with the experiment for the repartition of bubbles near the plating boundary. We believe that in the future the computer code will be used to design industrial and experimental systems. However, as to the measurement of the deposition rate, It takes at least one hour to obtain an observable thickness of plating. In this case, bubbles have accumulated everywhere in the micro-channel and there is ground for an extension of the present code with a level set or phase field model which tracks the main liquid to gas interface. To establish a mathematical model suitable for a larger time simulation is left as future work. the growth rate of bubbles governed by the local mass loss prescribed by Equation (3) can be computed by the relation where N q is the amount of bubbles in a local volume V 0 . Therefore, we have the following formulae on A g and A l , respectively The quantity R B is useful when the fluid velocity is large enough so that each bubble will not stay at the observed physical domain, because every bubble has not been far from the state that is just after nucleation.
Given a small cube V 0 of size |V 0 | = d × d × d and a typical radius R B , the ratio of its surface area and volume is 4πN q R 2 B d 3 , where N q can be estimated by Therefore, if d is small enough so that the physical quantities in F α defined in Section 2.3 can be assumed uniform, then we have the approximation Similarly, F g ≈Ṡ g u g (A6)