Hydrodynamic Simulation of Gas – Particle Injection Into Molten Lead

Departamento de Ingeniería en Metalurgia y Materiales, Instituto Politécnico Nacional, Escuela Superior de Ingeniería Química e Industrias Extractivas – ESIQIE UPALM, 07738, México D.F., México Centro de Investigación en Metalurgia y Materiales, Universidad Autónoma del Estado de Hidalgo, Carretera Pachuca-Tulancingo Km 4.5, 42084, Pachuca-Hgo, México Departamento de Materiales, Universidad Autónoma de Coahuila, Boulevard V. Carranza y José Cárdenas Valdés 25280, Saltillo, Coah, México


Introduction
Lead is a non-ferrous metal, which is used in soil pipes, as antiknock in gasoline, cable sheathings, paints, ammunition and different alloys; however, by far the biggest use of lead, worldwide, is for the lead-acid battery due of its characteristics (conductivity, corrosion resistance and the special reversible reaction between lead oxide and sulfuric acid) 1 .Primary lead is produced from lead ores.Nowadays, lead recovered from exhausted batteries is carried out by the pyrometallurgical route, which may cause environmental problems like the emission of considerable amounts of dust containing lead particulate and SO 2 emission into the atmosphere 2 .The initial smelting and refining of secondary lead is generally done in a blast furnace, rotary furnace or reverberatory furnace, which produces different qualities of lead bullion.The process of refining consists in the removal of impurities in these furnaces or in a separate operation in individual kettles.The refining stage involves three main processes: a) the copper drossing; b) softening, and c) desilvering.The copper drossing is a simple cooling process in a large kettle in which impurity compounds (mostly antimonies, arsenides, and stannides) separate as either solid or immiscible liquid phases due to the decreasing solubility of the phases at lower temperatures.The impurities are removed by skimming during cooling, and the molten lead is finally cooled to about 340 °C.To remove the last amount of copper between 20-200 ppm, sulfur may be stirred into the melt, and cooling is continued to very near the freezing point of the melt (~330 °C).The phases to be removed float on the lead surface and are called mattes or speiss , depending on the impurities they contain.They are skimmed off, tapped off, or otherwise removed [3][4][5] .At this point antimony, arsenic, tin, bismuth and any noble metals (primarily silver) are the major impurities remaining.The first three of these elements are removed in an oxidation process called softening.Silver is removed last and is generally removed only in primary lead refining.Desilvering can be done by successive cooling and crystallization steps and/or by the Parkes process 6,7 .In the latter process, zinc is stirred into the melt that is then cooled.A separate zincrich phase containing the silver, any gold and some lead separates to the top and is removed.The zinc remaining in the melt is then removed by an oxidation or a vacuum dezincing (distillation) process.At this point, the bullion is relatively pure lead.The procedure of adding powder reagents, sulfur and zinc to remove copper and silver onto the molten lead surface and further stirring is an inefficient step in the lead refining 8,9 .A technique used to increase the efficiency of lead refining is the submerged injection of solid reagents into the melts.Powder injection has become, in recent years, an important tool in metallurgical operations focused primarily for use in the steelmaking processes.The injection efficiency depends mainly on the interaction (solid particle and conveying gas) with the liquid metal.Important factors to consider are the real contacting interfacial area and the residence time between the liquid metal and solid particle 10 .The powder injection process considers two possible reaction zones 11 .These are, the transitory reaction zone (reaction between the particle injected and liquid metal) and the permanent reaction zone (reaction between the metal-slag interfaces).The residence time of the particle is determined by their velocity relative to the melt and by the bulk flow velocity of the melt.If the particles have low settling velocities, their motion will be controlled by the bulk circulation generated by the bubble plume.If they remain in the region of the lance tip upon injection, they will be rapidly carried by the plume to the upper surface of the melt with low efficiency 12 .Some studies have been carried out on the refining metals with mathematical simulation.Ohguchi et al. 13 studied the desulfurization process in torpedo cars.They establish a kinetic model in terms of the mixing time to study the kinetics of powder injection and the metal-slag reaction considering a fixed volume of slag.Schwarz 14 established mathematical models to simulate the air-water interaction and the bath deformation in gas injected melts.Vargas et al. 11 analyzed the pre-treatment of pig iron by powder injection of calcium oxide as a reactant and established a mathematical model for multicomponent reaction kinetics.Plascencia et al. 8 carried out experiments to remove copper from lead bath by conventional top addition and injection of sulfur powder.In addition, they used thermodynamic and kinetic models to simulate the lead refining.Scheepers et al. 15 modelled the desulfurization process of molten iron through the injection of calcium carbide by a submerged lance.The modelling considered the kinetic, thermodynamic and transport processes to predict the sulfur levels in pig iron for the Mittal Saldanha Steel plant.In the literature described above, modelling of metallurgical processes for either multi-phase or multicomponent flows have been carried out with the use of computational commercial software.When two immiscible fluids are in contact, the interface often plays a central role in the fluid dynamics and rheology of the mixture system.In general, the mathematical treatment of such moving boundary problems is complicated by the fact that the position of the interface is not known a priori.Instead it evolves according to the flow within both fluid components.The phase-field method stands out in its treatment of the interface as a physically diffuse thin layer in two immiscible fluids.This formalism systematically incorporates complex fluid rheology as well as interfacial dynamics into a unique CFD tool for multi-phase and multi-component complex fluids 16 .Yue et al. 17 adopted a finite-element method which easily allows complex domains and unstructured meshing schemes using the phase field parameter as criteria for the dynamic refinement and coarsening of the grid.Chen et al. 18 implemented an accurate and efficient semi-implicit spectral method for solving the phase field equations.The time variable was discretized by using semi-implicit schemes which allow much larger time step sizes than explicit schemes; the space variables were discretized by using a Fourier-spectral method whose convergence rate is exponential in contrast to second order by a usual finite difference method.Accurate and efficient numerical methods have been developed to solve the coupled Cahn Hilliard-Navier Stokes system with the resolution of very thin layers (diffuse interfaces) to capture the physics of the problems studied 19,20 .A multi-phase simulation was carried out with hydrodynamic behavior involved in the particle injection with nitrogen as conveying gas through a lance into molten lead.The model consists of a Navier Stokes system coupled with a Cahn Hilliard equation based on the phase-field method.A Petrov-Galerkin method for the numerical approximation was used in this work.The effects of operating parameters like gas flow rate, lance height, and different kettle and lance dimensions on the residence and mixing time were obtained.The understanding of the hydrodynamic behavior of the lance-kettle system will be useful for further studies on the mass transfer efficacy on the impurity transfer in the copper drossing and desilvering process for lead refining.

Hydrodynamic behavior
The phase field method offers an attractive alternative to more established methods for solving multiphase flow problems.Instead of directly tracking the interface between two fluids, the interfacial layer is controlled by a phase field variable, (Φ).A detailed development of the phase field model was given previously [16][17][18][19][20] .Only a brief summary will be presented here.In a mixture of two incompressible fluids, the interface has a small but finite thickness.Inside of this interface the two components are mixed and store a mixing energy.The mixing-energy density (f mix ) is represented as a function of (Φ).

(
) ( ) If a double-well potential is used, f 0 is expressed as follows: ( ) where f 0 is the bulk energy, c is a capillarity width representing the interfacial thickness and λ is the energy density parameter.(Φ) =-1 at liquid region, Φ=1 at gas region and Φ=0 at the interface.The bulk energy shows total separation from the phases in domains of pure components instead of the complete mix of the components.The profile of Φ across the interface is determined by the competition between the two effects.Since the mixing energy represents molecular interaction between the two phases, the interfacial tension is considered equal to the surface tension coefficient σ T (N/m) which is calculated as follows If the interfacial thickness (c) shrinks towards zero, so should the energy density parameter (λ); their ratio gives the interfacial tension in the sharp-interface limit.If the diffuse interface is not at equilibrium but it is relaxing according with the evolution on Equation 6given below, the interfacial tension is not constant, and in fact it reflects the reality that the interface has its own dynamics which cannot be embodied by a constant interfacial tension except under limiting conditions.The hydrodynamic system describing the mixture of two Newtonian fluids with a free interface can be done with the Navier-Stokes equations and the phase field method.In the solution domains, we consider the continuity and momentum equations in the usual form: where σ 1 is a deviatoric stress tensor and the field variables, velocity u and pressure p.The surface tension force is added to Navier-Stokes equations as a body force by multiplying the chemical potential of the system by the gradient of the phase field variable (as it is shown in Equation 16).The evolution of the phase field variable is controlled by the Cahn-Hilliard equation.(6)   where δF/δΦ = G is the chemical potential (Pa) and g is the mobility (m 3 s kg -1 ).The mobility determines the time scale of the Cahn-Hilliard diffusion and must be large enough to retain a constant interfacial thickness but small enough so that the convective terms are not overly damped.The mobility is determined by a mobility tuning parameter that is a function of the interface thickness 16 .
where c is a capillary width that scales with the thickness of the diffuse interface χ (m).The mixing energy considered to produce the interface structure is obtained with Equation 8 16 .
( ) The parameter gλ determines the relaxation time of the interface and it is considered according with Jacqmin 21 .For this case, the two phases have different densities and viscosities; therefore to obtain an incompressible flow solution, the volume of each fluid should remain constant in time.It should not be allowed to change because of diffuse flux.When an amount of fluid 1 flows out of an infinitesimal volume element due to interfacial diffusion, there will be also an amount of fluid 2 of the same volume that would enter the volume element at the same time, and viceversa.To obtain the continuity equation of a divergence-free velocity , a volume average velocity u is defined 16 .
Where n denotes the mass flow rate (per unit volume).Equation 9represent a mixture of the components by advection without diffusion, thus, the velocity at a special point is defined as that of the component occupying that point; it is spatially continuous and remains solenoidal.The Cahn-Hilliard equation forces the phase field variable (Φ) to take a value of 1 or -1 except in a very thin region on the fluid-fluid interface.The application model decomposes the Cahn-Hilliard equation into two second-order partial differential equations 22 .
where f ext is a user defined free energy (J/m 3 ).In this work, the external free energy was zero.In order to model the flow of two different immiscible fluids, a two-phase flow in the phase field method is required.It must be taken into account differences in the two fluids such as densities and viscosities and includes the effect of surface tension and gravity.For this simulation, molten lead and nitrogen gas were considered as fluid 1 and 2, respectively.Therefore, the volume fraction of fluid 2 is computed by the following equation.
[ ] ( ) The min and max operators are used to consider that the volume fraction has a lower limit of 0 and an upper limit of 1.The density and viscosity are computed by Equation 13and 14, respectively.
( ) where ρ 1 and ρ 2 are the densities and η 1 and η 2 are the viscosities of fluid 1 and fluid 2, respectively.For incompressible fluids ( 0) u ∇ ⋅ = , in addition of Equations 13 and 14, the Equation 15 was also considered and solved simultaneously with the phase field method.
( ) ( ) The two forces on the right-hand side of Equation 15 are due to gravity and surface tension.The surface tension and gravity forces are represented in Equation 16 and 17, respectively.
where g is the gravity vector.The mathematical simulation considered the injection of an inert tracer into the molten bath to determine its time dependence.The mixing time is determined from a convection-diffusion model represented by Equation 18.
where b is the concentration (kg m -3 ), D is the diffusion coefficient (m 2 s -1 ) and u is the velocity (m s -1 ).The velocity was obtained from the phase field model solution (Equation 8).The Reynolds number was in the range from 4000 to 26000 for the experimental conditions of this work.Therefore, the laminar diffusion of tracer was excluded due to its low concentration, so the turbulent transport predominate and the diffusion term (D) was obtained from the phase field simulation as a turbulent viscosity (v T ) represented by Equation 19.
where ρ is the density, C µ is the turbulence modelling constant (0.09), ε is the turbulent dissipation rate, k is the turbulent kinetic energy 23 .

Numerical solution
Equations given above were solved with the commercial software Comsol Multiphysics 4.1 22 .The software considers a finite element flow solver that Hu et al. 24 have used for simulating particle motion in Newtonian and viscoelastic fluids and an adaptive meshing scheme based in the work of Freitag et al. 25 .The Cahn Hilliard equation was solved with a combination of the Petrov-Galerkin/Compensated (P-G/C) method 22 to discretize spatially to improve stability and the backward difference schemes.P-G/C method uses a modified intrinsic time scale that completely eliminates the streamline diffusion contributions in regions dominated by diffusion.A semi-implicit scheme in which the principal elliptic operator is treated implicitly to reduce the associated stability constraint was used for the time step, while the nonlinear terms are still treated explicitly to avoid the expensive process of solving non linear equations at each time step.The accuracy in time can be significantly improved by using higher-order semi-implicit schemes.For this case, a second-order backward difference formula (BDF) was used for the time derivatives and a second order Adams-Bashforth (AB) scheme was used for the explicit treatment of non linear transport terms 18,19 .The simulations developed in this work considered grid and time-step refinements to ensure convergence.For the time step, the Courant-Friedricks-Lewy condition (CFL condition) was used 16,26 and it is represented in Equation 20.
where u is the velocity (m s -1 ), ∆t is the time step (s) and h is the length interval (m).For this simulation, a time step such that 2 ≤ CFL ≤ 10 was used to obtain a good relation between efficiency and accuracy.A moderated refined mesh was used to achieve the mass conservation at a discrete level with the finite element method.The phase field theory introduces three model parameters, λ, c, g which must be pragmatically chosen according to the particular situation to be studied.The simulation must consider the accuracy in the prediction of the interface dynamics and also must preserve the volume of fluids.Yue et al. 17 have validated values of these parameters with an adaptive mesh, which allows handling large interfacial morphological changes.

Model condition
The experimental conditions and the relevant physical data used for the model resolution are given in Table 1.The lance was located at the centerline of the kettle which was symmetrical as can be seen in Figure 1, so the domain solution was represented as a two dimension (2D) slice.The boundary conditions of the kettle contours are shown in Figure 1.The wall in contact with the fluid interface considers a no-slip boundary condition.Therefore, the Equations 21 to 28 were considered.
where n is the unit vector normal to the wall.The volume fraction V f must be specified at the inlet boundaries, normally either 0 or 1.
( ) (25)   For outflow boundaries, the following conditions were imposed.During the initialization step, the boundary condition sets the phase field function to 0.
The 2D geometry was meshed in 929 fine triangular elements and 8329 degrees of freedom for greater detail of the phenomena resulting from the simulation.The phase field model allows the use of a fixed Eulerian mesh to capture moving internal boundaries.A mesh with dense grids covering the interfacial region and coarser grids in the bulk are required.As the interface moves out of the fine mesh, the mesh in front needs to be refined while that left behind needs to be coarsened.Such adaptative meshing is achieved by using a general-purpose mesh generator GRUMMP 17,25 .The mesh generator provides triangular and tetrahedral elements for 2D and 3D, respectively.The mesh quality is improved through point insertion at the circumcenter of badly shaped cells having and angle smaller than a threshold θ min .If a proposed new point encroaches on a boundary edge, that vertex is not inserted.Instead, the encroached boundary edge is bisected.This process is repeated until all cells are well shaped.GRUMMP controls the spatial variation of grid size using a length scale (L s ), which specifies the intended grid size at each location in the domain.For this work, the grid size distribution is dictated by the need to resolve thin interfaces.The phase field variable Φ is constant (±1) in the bulk and varies steeply across the interface.A prescribed small grid size h 1 is imposed on the interface by making L s depend on ∇φ on every node.
( ) where h ∞ is the mesh size in the bulk, and the constant C controls the mesh size in the interfacial region: being the capillarity width.Values of C between 0.5 and 1 were used to ensure a thickness of the interface to contain the order of 10 grid points.Furthermore, the mesh size can be set to different values for the two bulk fluids.

Model validation
In order to validate the model results, the residence time was also calculated based on the efficiency of the transitory reaction given by Ohguchi 13 with the Equation 30and experimental data reported in the desilvering process during lead refining 28 .
where E is the efficiency of the transitory reaction, A pow is the particle surface area, k e is the mass transfer coefficient, t e is the residence time, L eq is the equilibrium coefficient of silver concentration in the particle and V pow is the particle volume.Rearranging terms and solving for the residence time, Equation 30 is represented as follows.
log (1  )   eq pow e pow e

L V E
A k The value of the efficiency of the transitory reaction is in the range from 0 to 1, when E = 1 indicates that the silver particle reacts totally with the molten lead in the residence time.The efficiency of the transitory reaction was considered as 0.25 28 .The equilibrium coefficient of silver concentration in the particle (L eq ) was obtained with Equation 32.
where (%Ag) is the silver particle concentration when it reaches the liquid surface, [%Ag] eq is the silver content in thermodynamic equilibrium at the process temperature in the bath metal.For this case L eq was considered in 230.4147 at 480 °C28 .The mass transfer coefficient 29 was determined with Equation 33.where r is the radius of the particle (m), D is the diffusion coefficient (m 2 s -1 ), pε is the specific mixing power (watt kg -1 ) and ν is the kinematic viscosity (m 2 s -1 ).The specific mixing power 11 is represented by Equation 34.
where M is the mass of molten metal (kg), Q is the volumetric flow of conveying gas (m 3 s -1 ), T is the bath temperature (K) and H is the bath depth (m).Substitution of Equations 33 and 34 into 31 gives the residence time.
Table 2 shows the parameters considered in the residence time determination calculated with the Ohguchi model.

Gas flow rate
The mathematical simulation was carried out to a flow rate injection of 3.3155, 6.6310 and 9.7887 m s -1 , and it was validated with the results reported by Rodríguez 28 .The hydrodynamic results are shown in Figures 2, 3 and 4 for a flow rate injection of 3.3155, 6.6310 and 9.7887 m s -1 , respectively for 0.01, 0.1, 0.2 and 0.3 s.The range of time studied in the simulation considers only the effect of the transitory reaction.Figures 2, 3 and 4 presents the flow pattern in the central vertical plane of a lance -kettle symmetrical system in two dimensions.The results of Figures 2 to 4 show the gas bubble formation according to the evolution of the volume fraction of nitrogen gas as time proceeds.The volume fraction of nitrogen gas was increased when time was increased and shows the interfacial layer formed between liquid lead and nitrogen gas which is represented by the phase field variable (Φ) in the Cahn-Hilliard equation.The results show the velocity vectors and black flow lines that indicate the flow trajectory and   direction of fluid velocity, as well as a vortex formation zone produced by a metal recirculation flow.This vortex formation was clearly observed located at the central part of the kettle, caused by the rotary motion of molten metal in contact with nitrogen gas.The vortex formation zone was increased when the rate injection was increased.A dead zone formation was observed also at the bottom of the kettle and it was increased to the lower rate injection.

Residence and mixing time
In order to determine the residence time, the partial differential equations were solved time dependent.The residence time indicates the time when the particles are in touch with the liquid metal and therefore represents the time when the transitory reaction is carried out.Figure 5 shows the residence time determination in 0.312, 0.280 and 0.262 s for the injection rate of 3.3155, 6.6310 and 9.7887 m s -1 , respectively.As can be observed, the residence time was decreased when the injection rate was increased.Figure 5 also shows the trajectory of the particle, as can be observed; increasing the gas flow rate increased the injection velocity of the particle and the amount of potential energy dissipated, which contributes to a deeper penetration of the jet. Figure 6 shows the residence time obtained with the Ohguchi model from experimental data and with simulation applying the phase field model.An absolute difference of 0.012, 0.001  and 0.005 s was obtained for both models for the injection rate of 3.3155, 6.6310 and 9.7887 m s -1 , respectively.
The mathematical simulation considered the injection of an inert tracer into the molten bath to determine its behavior dependent time.The results of these measurements were the F curves shown in Figure 7 for the injection rates considered.The mixing time was determined in 0.062, 0.064 and 0.066 s for the injection rate of 3.3155, 6.6310 and 9.7887 m s -1 , respectively.

Lance depth
The effect of the lance depth on the residence and mixing times for the gas flow rate 3.3155 m s -1 is observed in Figure 8.The lance depth was varied at 30, 50, 70 and 90 % of the kettle height.The results show the velocity vectors and flow lines which indicate the flow trajectory and direction of fluid velocity.Figure 9 shows the results of the mixing and residence time for the three gas flow injections that have been analyzed.

Lance and kettle dimensions
The variation effect of lance and kettle dimensions on the residence and mixing times for the gas flow rate 3.3155 m s -1 is observed in Figure 10.The lance and kettle measurements reported in Table 1 were affected by a scaling factor to increase and decrease their original dimensions by 50 %.Figure 10 shows the residence and mixing determination times for different lance (Figure 10a) and kettle (Figure 10b) sizes for a flow rate of 3.3155 m s -1 .the gas flow rate is increased, a continuous jet is promoted, inclusively, a swirl motion on the liquid surface is reached 31 .It was observed the formation of a gas channel during the gas injection instead of the formation of discrete bubbles under all of the simulated conditions.In this work, the effect of the change of pressure in the gas bubble formation was considered in the nitrogen density as is reported in Table 1.

Residence and mixing time
The residence time represents the time required for the particles leaving the lance tip until they reach the liquid surface supported with the conveying gas.If the residence time decreases, the particles will have less opportunity to react with the liquid bath decreasing the refining metal processes efficiency.Figure 5 shows that when the primary bubble is eliminated by their dissolution on the liquid surface, entrainment of the particle into the plume zone produced by the insoluble gas, resulted in their rapid convection to the upper surface of the bath, where they remained unless the recirculation flow was sufficient to pull them back into the bulk liquid.The vortex formation zone keeps the particle in circulation in the bath when the gas flow rate was increased, as is observed in Figures 5b and 5c.Due to the substantial density difference between the nitrogen gas and the molten lead, the buoyancy force will be dominant

Gas flow rate
Figures 2, 3 and 4 show that on exit from the lance, at the time of 0.01 s, the nitrogen jet expands rapidly and penetrates only a short distance into the molten lead before rising vertically, at the time of 0.1 s.At this time, the gas plume is almost symmetrical; however, for the times of 0.2 and 0.3 s, the gas plume is not symmetrical, this behavior is due to the internal pressure of the bubble changes as the bubble rises to the surface, as was reported by Kelvin's equation 30 .P g = P l + 2 σ T /r , where P g is the pressure inside the gas bubble, P l is the pressure in the liquid at the level of the bubble, σ T is the surface tension and r is the radius of curvature.If the super-saturation pressure with respect to the transferred component exceeds P g , the bubble will grow until it reaches a size where the buoyancy forces exceed the surface tension forces and detachment occurs.It has been reported a transitional range based on Reynolds number 30 and modified Froude number 31 , where the bubbles form a continuous jet or discrete bubbles according with the gas flow rate.Generally, the inertial force of gas injected in from the top lance in the low gas flow rate is not strong, therefore there is a forming of discrete bubbles and when enough to assume that the initial velocity of the bubble will be close to final, leaving the particle on the liquid surface.Figure 6 shows the validation of the residence time obtained with Oguchi and phase field models.It must be stressed that the modelling results show ideal conditions of the particular lance-kettle system.Otherwise, Ohguchi's model takes into account an efficiency related with experimental factors such as impurities (Ag, Sb, As, Cu, etc) contained in lead and the temperature control, which directly influences the viscosity of the molten metal and hence the hydrodynamic behavior.However, the differences between both results were considered acceptable for this simulation.The mixing time was obtained to have a better comprehension of the reactor efficiency and hence the process efficiency, considering that the mixing time is an indirect measure of the transfer of momentum from the gas to the liquid 32 .Figure 7 shows the mixing time determination, which represent the time required to attain a desired level of homogeneity in the bath, this means that in a given time, the concentration in the metal bath is constant.

Lance depth
Figure 8 shows that when the lance depth was increased, the volume fraction of nitrogen was increased as well as the molten metal recirculation flow.A metal recirculation flow produced by a vortex formation zone was observed when the lance depth was 70 and 90 %, while a large dead zone was located at the bottom of the kettle to a lance depth of 30 and 50%.Therefore, an increment in lance depth increased, the liquid recirculation velocity and the particle residence time, despite the ensuing increment of recirculation velocity due to the vortex formation.The results of the mixing time show an increase when lance depth was increased; however, an opposite behavior was observed for the gas flow rate of 3.3155 and 6.631 m s -1 with a lance depth of 30% of the kettle height.This behavior is due to the faster dissipation of potential energy for the cases of lance depths of 30 and 50%, which allows a lower penetration of the jet, as observed in Figures 8a and 8b.A similar behavior was observed when the gas flow rate was increased for the different lance heights in Figure 9.

Lance and kettle dimensions
The results in Figure 10 show that the residence and mixing time were decreased when the lance diameter was increased.Considering that the variation of an average bubble diameter depends on the gas flow and the orifice's diameter, for this case, the gas flow remained constant and the orifice's diameter was modified.Hence, a large area in the lance tip generates a bigger bubble formation.The circulation of gas in a bigger bubble results in the interface movement of the liquid and in a reduction of the drag force on the bubble.Therefore, the buoyancy force is increased, which allows the bubble to rise faster, and   to decrease the residence and mixing time.To transfer an adequate momentum and to stir the molten metal in a bigger kettle, an increment in the lance depth 15 or gas flow injection is necessary.A bigger kettle shows an increase in the mixing and residence time for a constant flow rate.A better comprehension of the hydrodynamic behavior of molten metal-gas injection systems will be useful for further studies on the mass transfer of impurity removal on metal refining.The obtained results demonstrate that the injection rate of 3.3155 m s -1 and the increase of the lance height (deeper position) generate the higher residence and mixing time.This behavior is necessary to obtain an efficient disengagement of the particle from the carrier gas within molten lead.In addition, the kinetic behavior of the particle-metal interaction must be considered as well as operating parameters of the process in order to obtain the highest efficiency in the lead refining process.

Conclusions
A mathematical simulation of the operating parameters (lance height, lance and kettle dimensions) was carried out to establish the hydrodynamic behavior in the lance-kettle system for the lead refining process.The residence and mixing time were suitably simulated with the phase-field method and the Cahn Hilliard equation.The effect of gas flow rate, lance depth, and different kettle and lance dimensions were evaluated to obtain the residence and mixing time.The residence and mixing time were decreased when the injection rate and the lance diameter were increased.Therefore, the particle will have less opportunity to react with the liquid bath decreasing the refining metal processes efficiency.When the lance height and kettle dimensions were increased, the residence and mixing time were also increased.The mathematical modelling was validated with the determination of the residence time based on experimental data and Ohguchi's model.The results obtained with both models were similar for the gas flow rate of 3.3155 m s -1 .In order to reach the highest efficiency in lead refining, the residence and mixing time must be considered based on the operating parameters of the process.

Figure 1 .
Figure 1.Schematic diagram of the lance-kettle system and boundary conditions.

Figure 6 .
Figure 6.Residence time determination with the phase field and the Ohguchi models for different flow rates.

Figure 5 .
Figure 5. Residence time of 0.312, 0.28 and 0.262 s of the particles for a flow rate injection of a) 3.3155, b) 6.6310 and c) 9.7887 m s -1 , respectively.

Figure 7 .
Figure 7. Mixing time of 0.062, 0.064 and 0.066 s determined with F-curves for a flow rate injection of a) 3.3155, b) 6.6310 and c) 9.7887 m s -1 , respectively.

Figure 8 .
Figure 8. Residence time of 0.160, 0.174, 0.312 and 0.345 s of the particles for a flow rate injection of 3.3155 m s -1 at lance depth of 30, 50, 70 and 90 % of the kettle height, respectively.

Figure 10 .
Figure 10.Residence and mixing times for different lance and kettle sizes for a flow rate injection of 3.3155 m s -1 .

Figure 9 .
Figure 9. Residence and mixing times for different lance depth and flow rate injection of a) 3.3155, b) 6.6310 and c) 9.7887 m s -1 , respectively.

Table 1 .
Parameters considered in the residence time calculation with Ohguchi model.

Table 2 .
Experimental conditions and physical properties.