Article pubs.acs.org/IECR
Modeling the Dynamics of Precipitation and Agglomeration of Oxide Inclusions in Liquid Steel N. Rimbert,*,† L. Claudotte,‡ P. Gardin,‡ and J. Lehmann‡ †
LEMTA, Université de Lorraine, CNRS, ESSTIN, 54500 Vandœuvre-lès-Nancy, France ArcelorMittal Global R&D, Maizières Process, 57283 Maizières-lès-Metz, Cedex, France
‡
ABSTRACT: Deoxidation and refining of liquid steel is the site of many complex phenomena, such as liquid−gas−solid multiphase flow, turbulence, precipitation, agglomeration, and extraction of multicomponent inclusions. The steel industry is very interested in developing a framework to integrate these different mechanisms. This is based on three pillars: chemical thermodynamics, population balance, and transport phenomena. Accordingly, a homebrewed thermo-kinetics code is interfaced with a commercial Computer Fluid Dynamic code through a population balance module. The population balance module resorts to the classic Quadrature Method of Moments extended to multicomponent inclusions. This paper is focused on a comparison between three-dimensional numerical results obtained with or without taking into account precipitation kinetics. This distinction is important since under an assumption of local equilibrium, harmful inclusions are completely dissolved before the continuous casting occurs, whereas when using precipitation kinetics, harmful inclusions can be present in the final product. This is consistent with industrial samples.
1.1. Argon Bubbling of a 3D Liquid Steel Ladle. During steel deoxidation, the ladle can be stirred for instance, by the injection of argon bubbles at its base. This is the technique that has been modeled in this study and the parameters relevant to the ladle geometry, porous plug size, position and flow rate are given in Table 1. The bubble nozzle inlet is located in the
1. INTRODUCTION Inclusion cleanliness is one of the key points of high-grade steel elaboration, such as interstitial free (IF) steel widely used in automotive industry.1 More details on this important industrial topic can be found in a recent review by Zhang and Thomas.2 Previous work mostly focused on the evolution of one component inclusions (cf., ref 3, for instance, which also contains a small review). Switching to multicomponent inclusions has been the incentive of a new kind of integrated modeling named multi-QMOM, based on software that couples the solving of chemical kinetics, fluid motion, and population balance equation (PBE) of inclusions.4,5 Though the code can theoretically handle any steel composition, for the sake of both simplicity and computational time, only combinations of alumina and titanium oxides have been considered so far.6 One of the motivations behind this choice is that these compounds are known to increase clogging of the nozzle during continuous casting.1 The aim of the present paper is to present new results stemming from the three-dimensional multi-QMOM model, which clearly show that the combination of fluid transport and of the kinetic aspects of nucleation, growth, and aggregation of particles are essential for assessing the occurrence of out-ofequilibrium harmful inclusions, such as alumino-titanate oxides. This Article starts with a short study of the hydrodynamics of the ladle chemical reactor under consideration. Former results of the coupling of a chemical equilibrium code (CEQCSI) with the computer fluid dynamic (CFD) commercial code that showed spurious redissolution of AlTiOx are reiterated. Next, in a second part the multi-QMOM framework is introduced. Finally, in the third and last part, results of alumino-titanate oxides’ precipitation, including chemical kinetics, are detailed and it is shown that this can account for observed occurrences of AlTiOx. © 2014 American Chemical Society
Table 1. Dimensions and Reference Values for the Steel Ladle ladle diameter (m) ladle height H (m) position of the porous plug diameter of the porous plug (m) argon injection temperature TN (K) liquid steel temperature T (K) liquid steel density ρsteel (kg/m3) liquid steel dynamic viscosity (Pa·s) Ar flow rate Qg,N (NL/min) Ar density ρAr at T = 1873 K (kg/m3) stirring power (W/t) bubble size DB (m)
4 3.7 centered 0.014 300 1873 7 000 0.0062 300 0.28 ∼23 0.002
middle of the ladle base. The computation was done using Ansys/Fluent discrete phase model so that the bubbles are tracked as lagrangian particles. One of the main reasons to use this lagrangian approach over competing eulerian modeling is that, provided gas hold-up is low, the same mesh can be used, independently of the injection location. The forces acting on the bubble are the added mass, the drag from the surrounding liquid steel, the lift force and the buoyancy force induced by Received: Revised: Accepted: Published: 8630
November 25, 2013 April 2, 2014 April 5, 2014 April 5, 2014 dx.doi.org/10.1021/ie403991e | Ind. Eng. Chem. Res. 2014, 53, 8630−8639
Industrial & Engineering Chemistry Research
Article
where Cb is an adjustable parameter fixed as 0.05 in this study and α is the bubble volume fraction. For ε, the extra term reads ε Sε = Cε Sk (7) k
gravity. Variation of bubble density with ferro-static pressure is also taken into account. The fluid phase is treated as a continuum by solving steady Reynolds-Averaged Navier−Stokes equations. The dispersed phase is solved by tracking a large number of bubbles through the calculated flow field. The bubbles’ trajectories are computed individually at specified intervals during the fluid phase calculation, and can exchange momentum with the carrier phase (two-way coupling). Continuity and momentum equations for the continuous phase are written, respectively, as
∂ ρUj = 0 ∂xj ρ Uj
with Cε = 1.2. The drag coefficient given in eq 7 was initially proposed by Tomiyama,11 who determined a set of equations for different interface situations: pure, slightly contaminated and contaminated systems. The contaminated system was considered because this is always the case in water and probably also in liquid steel, and the following expression was chosen and implemented:
(1)
∂ ∂P ∂ Ui = − + ρg i − ρ⟨u′iu′j⟩ + Fi ∂xj ∂xi ∂xj
⎛ 24 8 Eo ⎞ ⎟ C D = max⎜ (1 + 0.15Re 0.687), ⎝ Re 3 Eo + 4 ⎠
(2)
where ρ is the fluid density, P the liquid pressure, U the average liquid steel velocity, and u′ the liquid fluctuating velocity. Closure for the Reynolds stress ⟨u′i u′j ⟩ is obtained using a realizable k−ε turbulence model.7,8 The momentum transfer F from the continuous phase to the dispersed phase is computed by examining the change in momentum of a bubble as it passes through each control volume. It is given by FD =
With the Eötvos number defined as Eo =
(3)
m=N
F=
∑
FD(v − U)ṁ BΔt
∂ ∂ ∂ ⎛ ∂YC ⎞ (YC) + (UY ⎜Γ ⎟ + YĊ i C) = ∂t ∂xi ∂xi ⎝ ∂xi ⎠
(4)
m=1
where ρB, DB, ṁ B, v, and N are, respectively, the following bubble properties: density, diameter, mass flow rate, velocity, and number per control volume. μ is the liquid steel dynamic viscosity, CD the bubble drag coefficient, and Re the Reynolds number based on bubble diameter, slip velocity, and liquid kinematic viscosity. The bubble trajectory is predicted by integrating the equation of motion, written below
ρ dU ρB dt
(9)
(10)
with Γ = (νt/SCt) and Schmidt number SCt = 1, Ẏ C is the chemical source term of species C. Figure 1 shows the average flow field obtained after the computation has converged toward a steady state. Convective fluxes are computed using a QUICK second order scheme.4,5 However, it is difficult to validate the result of such a computation since measurements are difficult to carry out in situ. To do so, we chose to compare the computed stirring power to Sundberg’s formula12 given by eq 11.
ρ −ρ 1 ρ d dv (U − v) = FD(U − v) + g B + dt 2 ρB dt ρB +
ΔρgDB 2 σB
where σB is bubble surface tension and Δρ the density difference. Finally, bubble growth with ferro-static pressure evolution between the gas injection point and free surface is taken into account (ideal gas law) thanks to the User Defined Function developed in Fluent. For the chemical species mass fraction YC evolution, the conservation equation expressed as follows is solved
3μC DRe 4ρB DB 2
(8)
(5)
where g is gravity acceleration. In this equation, the most significant forces (drag, buoyancy, added mass, and pressure gradient) are taken into account.9 The dispersion of bubbles because of turbulence in the fluid phase is predicted using a “random walk” model,7 which includes the effect of instantaneous turbulent velocity fluctuations on the bubble trajectories through the use of stochastic methods (commonly called the “eddy lifetime” model). However, this model assumes that the bubbles have no direct impact on the generation or dissipation of turbulence in the carrier phase; therefore some modeling corrections are necessary. Accordingly, a correction has been implemented both in k and ε equations, as was initially proposed by Boisson and Malin,10 taking into account the production Sk of turbulent kinetic energy because of drag forces Sk = 0.75C bC Dρα(1 − α)
|| U − v ||3 DB
Figure 1. Argon stirring of the steel ladle. Average flow field on a vertical plane. Porous plug injecting argon bubbles is located on the axis of symmetry.
(6) 8631
dx.doi.org/10.1021/ie403991e | Ind. Eng. Chem. Res. 2014, 53, 8630−8639
Industrial & Engineering Chemistry Research ε=
RTQ g,N ⎛ ⎛ ρ gH ⎞ ⎛ T ⎞⎞ ⎜⎜ln⎜1 + steel ⎟ + ⎜1 − N ⎟⎟⎟ msteel Vm,g ⎝ ⎝ Patm ⎠ ⎝ T ⎠⎠
Article
Table 2) are then progressively added over 10 s. They are considered to be discrete particles, which dissolve after 1 s in (11)
Table 2. Mass Fraction of the Different Elements Added to the CFD Computation
In this equation, g is the gravity acceleration, R is the gas constant, Vm,g is the molar volume of the injected gas, msteel is the total mass of steel. TN is the temperature of the injected gas, and Patm is the atmospheric pressure at the top of the ladle. This formula takes into account the gas expansion, and it can be seen that a vacuum at the top of the ladle, as during vacuum treatment or in Ruhrstahl−Heraeus (RH) technology, would greatly enhance the stirring power. The application of formula 11 to the present case leads to a value for ε of 22.4 W/t, while the CFD results predict ε ≈ 23 W/t. Therefore, it will be supposed in the remaining part of this paper that this corresponds to a realistic flow field in a steel ladle. Note that fluid dynamic equations are computed at the beginning of the simulation and that once a steady state has been reached the flow field is frozen. This will be used to transport the nucleated particles, which are considered as tracers since they are small and quite dilute. Moreover, the whole ladle will be assumed to remain at the fixed temperature of 1873 K. Figure 2 shows the concentration of a passive scalar (i.e., following eq 10 with no source term) injected locally into the
YO_dissolved YAl_added YTi_added
350 ppm 700 ppm 400 ppm
the liquid steel13 with a melting rate given by a dedicated User’s Defined Function of Fluent. Once melting is complete, the added element will react with dissolved molecular oxygen and nucleation will begin. The location of particle injection is at the top right of the ladle. Table 2 shows the mass fractions that were added to the ladle, as well as the initial value of dissolved oxygen. The value of the physical properties of the elements relevant to their transport as passive scalar in the ladle is shown in Table 3. Note that these scalars are passive as far as hydrodynamics is Table 3. Physical Properties of Liquid Steel and Solute Elements steel ρ (kg/m ) μ (kg/m/s) D (m2/s) 3
7000 6.2 × 10−3
O
Al
12.58 × 10−9
Ti
2700
4500
3.02 × 10−9
8.3 × 10−9
concerned, but are quite chemically active. For instance, precipitation of alumina is both highly exothermic and fast; so sink terms, for both oxygen and added elements, and source terms, for oxides, must be summed accordingly. Table 4 Table 4. Chemical and Physical Properties of the Oxides Al2O3 ρ (kg/m3) ΔG0 (J)
3990 −1 199 869 + 393.212T 2.5554 × 10−5
Vm (m3/mol) σ (J/m2) 0.1
Figure 2. Upper (lower) black line is the maximum (minimum) of the local concentration ratio to the average concentration. Dotted lines are values of the same ratio obtained for 11 virtual captors located throughout the ladle.
TiO2
Al2TiO5
4250 −645 874 + 224.346T 1.8799 × 10−5
3700 −1 976 384 + 663.581T 4.9151 × 10−5
0.1
0.1
summarizes the Gibbs energy variations associated with the different chemical reactions considered
flow to evaluate the mixing time of the ladle. The concentration is nondimensionalized by dividing by the average concentration. Values of the scalars at 11 points scattered throughout the ladle (i.e., “virtual sensors” indicated by the 11 dashed lines in Figure 2) as well as the maximum and minimum values (solid lines) for the whole ladle are tracked. From Figure 2, the 95% mixing time can be estimated to be around 250 s. This is an important sequel to the ladle’s hydrodynamics since chemical species are considered as passive scalars in the present study. The mixing time represents the time needed to relax to a uniform concentration and therefore, for an instantaneous reaction, to the time needed to reach a global chemical equilibrium. 1.2. Precipitation of Titano-aluminate Oxides. As stated above, once the steady state is reached the flow field computation is stopped and Ti and Al addition elements (cf.,
2Al + 3O → Al 2O3(s)
(12)
Ti + 2O → TiO2(s)
(13)
2Al + Ti + 5O → Al 2TiO5(s)
(14)
This table also shows the relevant physical properties of the solid precipitates. In this study, two test cases have been analyzed. Case A corresponds to the (quasi) simultaneous introduction of aluminum and titanium which is known to be likely1 to lead to the production of AlTiOx. Actually, there is a 15 s delay between the two injections, titanium being added after aluminum. Case B corresponds to a long delay between the 2 additions (180 s). From industrial experience, Case B is unlikely to produce AlTiOx. 1.3. Preliminary Results Using a Local Chemical Equilibrium Hypothesis. First attempts14 to model the 8632
dx.doi.org/10.1021/ie403991e | Ind. Eng. Chem. Res. 2014, 53, 8630−8639
Industrial & Engineering Chemistry Research
Article
reactor (around 110s). Figure 5 shows the same evolution for case B. Case B is unlikely to produce AlTiOx and therefore, in
precipitation of AlTiOx were made using a different reactor technology (an RH reactor, whose flow field is depicted in Figure 3). The same modeling was used for bubble stirring as
Figure 5. Case B: 3D coupling between Fluent and Ceqcsi in an RH reactor. Time evolution of mass fractions is shown. Figure 3. Contour of the velocity magnitude (m/s) in the middle plane of the reactor. Air bubbles are injected into the right arm, creating circulation in a counterclockwise direction.
this case, the equilibrium computation agrees with the wellknown industrial fact. It can be seen that case A creates transiently much more AlTiOx than case B, which is encouraging and close to the industrial trends but, as expected, the AlTiOx species that appear are never stable and eventually dissolve in a rapid return to equilibrium. To conclude concerning these preliminary computations, the local chemical equilibrium hypothesis does not seem to be able to cope with the observed long-term occurrence of AlTiOx when Al and Ti are injected almost simultaneously. This was the incentive to develop a method based on out-of-equilibrium kinetics and on a population balance equation for the precipitates. This new method should allow the tracking of the growth of agglomerates both by diffusion and by collision, as well as their dissolution.
described in section 1.1 and was based on a local chemical equilibrium hypothesis, that is, each mesh volume is assumed to produce the equilibrium quantity of product deduced using CEQCSI software, that is, by minimizing Gibbs energy cell by cell, from the value given in Table 4. The precipitates that appear are considered to be new scalars passively transported by the flow (cf., eq 10). Figure 4 shows the time evolution of the different mass fractions of all chemical species for case A. Whereas case A is known to be likely to produce AlTiOx, some AlTiOx are indeed created from t = 25s on, but dissolve after a lapse of time following the addition of Ti close to the mixing time of the RH
2. MULTI-QMOM FRAMEWORK Multi-QMOM4 is a framework designed to turn one PBE with several inner degrees of freedom (such as size, composition etc.) into several PBEs with a single degree of freedom (hence the prefix multi). The PBEs are then solved using the Quadrature Method Of Moments developed by McGraw15 for aerosol applications and widely extended to other applications in chemical engineering.16 The advantage of this modeling over other classical methods such as the class method is that it is quick and efficient; however it cannot fully take into account particle size segregation.17 Figure 6 sketches the gist of this modeling; two kinds of particles are introduced: clusters and elementary particles. The population balance is, therefore, split into several parts, one set devoted to tracking the different nuclei and their evolution (or elementary particles), the other set devoted to tracking the aggregates (or clusters): ∂n p(L ; x , t ) ∂t
Figure 4. Case A: 3D coupling between Fluent and CEQCSI in an RH reactor. The time evolution of different species mass fractions is depicted.
+
∂Un i p(L ; x , t ) ∂xi
= Jp (L ; x , t ) + Gp(n p(L ; x , t )) 8633
(15)
dx.doi.org/10.1021/ie403991e | Ind. Eng. Chem. Res. 2014, 53, 8630−8639
Industrial & Engineering Chemistry Research
Article
Figure 6. The multi-QMOM framework developed to model multi component agglomerates. Geometry and composition are decoupled so that in this binary example the global population balance equation is split into three distinct population balances: one for agglomerates (open circles) and two for elementary components (gray and black disks).
where βij = β(Li, Lj) and is the length based aggregation kernel between two spherical particles of size Li and Lj; D/Dt stands for the total derivative using the fluid velocity. Np is the number of types of precipitate, Nq is the number of QMOM nodes, wi and Li are, respectively, weights and abscissas of QMOM and are classically obtained from the moments mk using the product-difference algorithm.4,15,16 Note that only the integral transform has been applied to the left side of eqs 15 and 16, while both integral transform and quadrature approximation have been applied to the source terms on the right side of eqs 15 and 16. To give a more physical interpretation to eqs 18 and 19, we emphasize that m0 is the number of particles per unit volume, m1/m0 is their average size, m2/m0 is their average surface, m3 is the total volume of particles per unit volume and is conserved when no precipitation or dissolution occurs. Moreover, term a stands for the nucleation rate of each individual species, while b governs the growth rate of nuclei that have already precipitated. Term c is the precipitation rate of the clusters. Likewise, term d is the growth rate of aggregates due to the diffusion growth of their elementary particles. Lastly, term e stands for the growth of the aggregates by aggregation. Finally let us point out that growth terms in eqs 18 and 19 are not used directly as they may create numerical instabilities related to the fact that the set of moments may become unrealizable i.e. no longer be a set of moments corresponding to a real distribution.18 Appendix A and ref 4 gives the exact method used. The main idea is to use the distribution representation as a set of Dirac peaks, move the peaks according to the growth or dissolution of the particles they correspond to and then compute the moments of the shifted distribution, thus ensuring that the new set of moment is realizable. 2.1. Nucleation. The nucleation rate J and the diffusion growth rate G are computed by the Mipphasolacido thermokinetic code, which is based, according to the classical nucleation theory (CNT)19 on the minimization of the Gibbs energy π ΔG(L) = σπL2 + ΔGv L3 (20) 6
∂nc(L ; x , t ) ∂Un i c (L ; x , t ) + ∂t ∂xi = Jc (L ; x , t ) + Gc(nc(L ; x , t )) + Ac(nc(L ; x , t ), nc(L ; x , t ))
(16)
where U is the fluid velocity, np(L;x,t) is the number density function (or NDF, that is, number of particles per size length per unit volume; sometimes also referred to as the particle size distribution or PSD) of elementary particle of chemical species p, of size L, located in a neighborhood of position x at time t. nc is the cluster NDF. Jp/c are nucleation source terms concerning either elementary particles p or clusters c and Gp/c are related growth terms. The independence hypothesis between size and composition during aggregation translates into an aggregation operator Ac that is applied only to clusters and not to elementary particles. The two sets of population eqs 15 and 16 are not independent and the cluster equation is linked to elementary particle equations through several constraints: (i) The nucleation term Jc is built from the different nucleation terms Jp as each new isolated nuclei is considered to be a single particle cluster. (ii) Similarly, the growth term Gc is built from the different Gp terms as the growth of each elementary particle of the cluster is assumed to contribute to the whole cluster’s growth. (iii) The total mass of the cluster is equal to the total mass of the elementary particles. However, eqs 15 and 16 are not directly solved in the present work. Let us consider the NDF moments and simultaneously introduce the Gaussian quadrature approximation that is used in the computation: mk =
∫0
+∞
The competition between the volume nucleation driving force and the surface creation energy cost leads to a minimum given by
Nq
Lk n(L) dL ≈
k ∑ wL i i i=1
(17)
where Li and wi are, respectively, the Nq Gaussian quadrature points and their corresponding weights. Applying this integral transform and quadrature to eqs 15 and 16 results in the following set of equations:4
ΔG* =
16πσ 3 3ΔGv2
(21)
where values of the interfacial energy σ can be found in Table 4, and ΔGv is obtained from 8634
dx.doi.org/10.1021/ie403991e | Ind. Eng. Chem. Res. 2014, 53, 8630−8639
Industrial & Engineering Chemistry Research
Article
where a0 is the interatomic distance, Dj the diffusion coefficient of element j in liquid steel, and xj its atomic fraction in the metal. 2.2. Growth. When interfacial kinetics is not a limitation factor for growth rate, the diffusion flux at the surface for an inclusion p of size Li is obtained as follows:
ΔGv = I(∑ CjiXi)/Vm (22)
i,j
where I denotes the driving force of precipitation. Cij is the number of atoms of component j in the oxide i, Xi is the mole fraction of oxide i, and Vm is the oxide mole volume. For instance in the present case, the oxides are respectively Al2O3, TiO2 and Al2TiO5 so the matrix C reads, in the (Al, O, Ti) coordinates: ⎛2 3 0⎞ ⎜ ⎟ C = ⎜0 2 1⎟ ⎜ ⎟ ⎝2 5 1⎠
Φp(Li) = k(Li)(Cp(l) − Cp(i))
C(l) p
where and are, respectively, the concentrations in mol· m−3 in liquid steel far from the precipitate and the interfacial concentration. When the growth is controlled by diffusion, C(i) p is equal to the equilibrium concentration. k(Li) is the mass transfer coefficient. In the present case, because of the size of the precipitates, the mass transfer coefficient is assumed to be
(23)
For a given stoichiometric oxide compound, the driving force I obeys the expression (for a binary oxide; we leave it to the reader to make the straightforward extension for ternary oxides) 1/ x + y ⎫ ⎧⎛ ⎪ (aE)x (aO) y ⎞ ⎪ ⎜ ⎟ ⎬ I = RT ln⎨⎜ ⎟ ⎪⎝ K ExOy ·a ExOy ⎠ ⎪ ⎩ ⎭
k(Li) =
Dp Li /2
(24)
Gp(Li) = Φp(Li)Voxide
(25)
where Z is the Zeldovitch factor, βn a frequency factor, Ns the number of nucleation sites per unit volume (here taken to be equal to the Avogadro number divided by the molar volume of liquid steel, as only heterogeneous nucleation is considered), and kB the Boltzmann constant. In the thermo-kinetic software,21 it is considered that all precipitates are spherical, with the diameter of new nuclei given by CNT:
βBr (Li , Lj) =
1 π (L*)2 2
Vm 5A
(26)
σ kBT
∑ Γi i
1/2 ⎡ ⎛ ((Li + Lj)/η)2 ⎞⎤ ⎟⎥ × ⎢fu , i + fu , j − 2fu , i fu , j ⎜⎜1 − ⎟ ⎢⎣ 601/2Reλ ⎠⎥⎦ ⎝
(34)
(27)
where u′ ≈ (2k/3)1/2 is the fluid turbulent (rms) velocity, Reλ ≡ (15u′4/εν)1/2 is the Reynolds number based on the Taylor length scale, ε the turbulent dissipation rate, ν the fluid kinematic viscosity and η = (ν3/ε)1/4 is the Kolmogorov length scale. The f u,i are the response coefficients24 of the particle i to the fluctuation of the fluid velocity. It is related to the correlation between the statistics of the carrier phase and the dispersed phase. The value of this coefficient is obtained23 through a biexponential approximation to the Lagrangian autocorrelation function of the fluid velocity. Appendix B gives details of the computation of these coefficients. By studying eq 34, it can be seen that when response coefficients are equal to one (i.e., the two colliding particles are perfectly following the carrier fluid flow), the Saffman and Turner formulation24 is recovered
(28)
where Γi stands for the number of oxide monomers i that can be added to the nucleus per unit of time and surface. Several models for Γi exist. Here, the model of LeGoues and Aaronson21 has been used as it is relevant for diffusion-limited growth. For a stoichiometric compound containing several elements, one gets Γi =
1 (∏ (Djxj)Cij )1/ ∑k Cki a04 i
(33)
βturb(Li , Lj) = (8π )1/2 (Li + Lj)2 u′
The frequency factor βn is given for a nonstoichiometric composition by βn = π (L*)2
2 2kBT (Li + Lj) 3μ LiLj
where Li and Lj are particle diameters and μ is the fluid viscosity. Turbulent Collisions. Zaichik et al.23 have developed a turbulent model to calculate the collision rates of bidispersed inertial particles in isotropic turbulence
The Zeldovitch factor in eq 25 is the number of oxide monomers that can be added or removed from the critical nucleus without changing the energy value by more than kBT19 Z=
(32)
where Voxide is the oxide mole volume in m3/mol. 2.3. Aggregation. The aggregation mechanism is classically assumed to be a second-order rate process, in which the rate of collision is proportional to the product of concentration of the two colliding species.22 In this work, three collision kernels and one collision efficiency are used, which will be presented now. The different kernels can be classified according to their range of influence. Brownian Collisions.
20
L* = −4σ /ΔGv
(31)
where Dp is the diffusion coefficient of the element p (given in Table 1) and Li the diameter of the precipitate. The diffusion growth rate is finally calculated using
where aE and aO are, respectively, the activities of element E and of oxygen referred to the 1% dilute solution in liquid iron, aExOy is the activity of ExOy in the oxide, KExOy is the solubility equilibrium constant of ExOy in liquid iron, and x and y are the appropriate stoichiometric coefficients. Finally the rate of homogeneous nucleation of inclusions is given by ⎛ ΔG* ⎞ J = ZβnNs exp⎜ − ⎟ ⎝ kBT ⎠
(30)
C(i) p
(29) 8635
dx.doi.org/10.1021/ie403991e | Ind. Eng. Chem. Res. 2014, 53, 8630−8639
Industrial & Engineering Chemistry Research βturb(Li , Lj) =
8π 15
ε (Li + Lj)3 ν
Article
(35)
In contrast, when response coefficients are close to zero, the turbulent collision kernel of Abrahamson26 is found (i.e., the particles are in the inertial part of the turbulence cascade and there is no longer any correlation between the flow and the particle collision) βturb(Li , Lj) = (8π )1/2 (Li + Lj)2 (vi′ 2 + v′j 2)1/2
(36)
where v′ is the fluctuating velocity of particles. Stokes Collisions. The Stokes collision kernel takes into account the settling difference, or more exactly in our case, the rising velocity difference between small and large inclusions. It reads: βSt (Li , Lj) =
2πg (ρ − ρc ) 9μ
|Li2 − Lj2|(Li + Lj)2
(37)
Figure 7. Comparison of the influence range of different aggregation kernels. The reference particle size is 1 μm when the reference temperature is 1873 K, and the reference turbulence dissipation rate is 2.2 × 10−3 m2/s3.
where ρ and ρc are, respectively, the fluid density and the estimated cluster density. Collision Efficiency. Higashitani et al.27 introduced a collision efficiency factor α(Li, Lj), defined as the probability that two agglomerates stick together when they collide. While this collision efficiency was initially developed to calculate turbulent kernel in viscous domain (eq 35), we similarly decided to apply the collision efficiency to generic turbulent kernel (eq 34). This coefficient is multiplied by the turbulent collision kernel to obtain a turbulent collision rate, which takes into account the van der Waals force and the lubrication forces needed to eliminate the liquid film between two colliding particles:28,29
model for oxide particles in liquid steel, thus the efficiency of Brownian collision is assumed to be equal to one. Uncertainty pertaining to this modeling is not a drawback as the influence of Brownian motion decreases as particle sizes increases and Higashitani’s formula applies. Moreover, there is a strong assumption that this efficiency is very high due to the possible existence of a gaseous bridge between the particles that sticks them together,28 hence no fragmentation of the aggregates has been considered.
⎛
αij = α(Li , Lj) = 0.8H
−0.18 ⎛ Li + Lj ⎞3⎞ ⎟ ⎟ ⎜36πμγ ⎜⎝̇ 2 ⎠ ⎟⎠ ⎝
0.18⎜
3. PRECIPITATION OF TITANIUM ALUMINATES IN THE MULTI-QMOM FRAMEWORK As can be seen in the preceding part of this paper, this new framework allows for the full kinetic simulation of the precipitation and aggregation of multicomponent oxides in a 3D steel ladle. Coupling with the thermo-kinetic code is carried out using Fluent User Defined Scalars (UDS) and User Defined Memory (UDM) and this allows for an out-ofequilibrium computation of the precipitation of the nonmetallic oxides. Figure 8 shows the results of this computation for case A, where the injection of titanium was delayed by 15 s relative to aluminum. In this case AlTiOx molecules are more likely to be found in the final product than when there is a long delay between additions of the elements. As can be seen from Figure 8, the AlTiOx aggregates are stable and do not dissolve during the mixing time of the ladle (250 s) after the injection has been realized. They seem, however, to slowly dissolve thereafter and to eventually reach a new kinetically driven equilibrium value from ∼350 s onward. In contrast, and as expected, the simulation of case B leads to the dissolution of the AlTiOx species (cf. Figure 9). In that case a value close to the ladle mixing time was used for the introduction of Titanium (i.e., the reactor was well-mixed, therefore decreasing the nucleation driving force) and the present out-of equilibrium computation gave results in agreement with a classical equilibrium assumption. The history of mixing is therefore an important parameter that interacts with kinetics of precipitation, growth and aggregation in an intricate way to explain the occurrence of
(38)
where H is the Hamaker constant. Its value is unknown for mixed oxides, but a value of 2.3 × 10−20 J is referenced30 for pure alumina in liquid steel at 1873 K. For turbulent flows, the strain rate is given by
γ turb ̇ =
ε ν
(39)
The collision efficiency factor (eq 38) has also been extended to be used with the Stokes kernel. In that case, the associated strain rate is related to the difference of terminal rising velocity by γ St ̇ =
2g |ρ − ρP ||Li − Lj| 9μL L
(40)
Finally, the aggregation kernel implemented is the sum of these three kernels: βij = βijBr + αijturbβijturb + αijStβijSt
(41)
Figure 7 shows the behavior and the relative intensity of theses kernels: Brownian motion governs the small scales, whereas turbulent collision is dominant in the intermediate range and floatation takes over at the larger scales. Higashitani’s formula gives very high collision efficiency for nuclei-sized particles (around 40 nm) and this efficiency decreases as the size of the agglomerate increases. To our knowledge there is no orthokinetic/Brownian agglomeration collision efficiency 8636
dx.doi.org/10.1021/ie403991e | Ind. Eng. Chem. Res. 2014, 53, 8630−8639
Industrial & Engineering Chemistry Research
Article
Figure 8. Case A: Results of the 3D coupling between Fluent, MultiQMOM and Mipphasolacido. A delay of 15 s occurred between the injection of aluminum and titanium.
Figure 10. Evolution of inclusion composition predicted by Fluent3D/MultiQMOM/Mipphasolacido for Case A.
oxides and liquid metal is, to our knowledge, unknown and varying (for instance with the oxygen content that acts as a surfactant). In present computations, it has been compulsory to maintain this at a low value (0.1 J·m−2) for nucleation to occur. This low value may be valid in the case of heterogeneous nucleation but in this case the number of nucleation sites should be reduced to reach the real number of heterogeneous nucleation sites. This is likely to somehow limit nucleation and therefore enhance growth by diffusion of already nucleated inclusions. Taking this into account will slightly increase the complexity of the modeling, as the number of heterogeneous nucleation sites is unknown, and a separate population balance equation (or at least a separate scalar equation) may eventually be needed. Another advantage of this approach would be to enable particles of detached slag/refractory to be taken into account, which could also serve as heterogeneous nucleation sites. Whether by dissolution of these slag/refractory inclusions in the liquid metal or by chemical evolution of the aggregates, this would allow computation of the evolution of the particles’ composition, in relation to the slag or refractory composition, as the two are known to evolve concurrently.31 However, this intricate problem is left for future work.
Figure 9. Case B: Results of the 3D coupling between Fluent, MultiQMOM, and Mipphasolacido. There was a delay of 180 s between the injection of aluminum and titanium.
these harmful out-of-equilibrium inclusions in real industrial cases. Lastly, Figure 10 shows the evolution of the composition of inclusions during the first 150 s of the simulation for case A. It can be seen that TiO2 and Al2TiO5 appears concomitantly but that TiO2 dissolves after 90 s (in agreement with Figure 8) and that most titanium is found in the Al2TiO5 form after this time. It can also be seen that variation of the mass fraction of Al2TiO5 is then related to variations of the mass fraction of both dissolved Ti and dissolved Al (cf. Figure 8). Many other results could be postprocessed, such as the location of the most harmful inclusions or the average size of the inclusion. However, using the present modeling, the size of the inclusion barely reaches the micrometer scale. This seems somewhat underestimated. One possible explanation could be that the largest inclusions in real ladles could be related to slag entrainment. But an alternative explanation may be linked to the fact that, for the time being, only homogeneous nucleation has been considered. Moreover, the surface tension between
4. CONCLUSIONS To conclude, using a recent framework involving MultiQMOM PBE resolution coupled to CFD and thermo-kinetics, this paper presents the first results showing that AlTiOx appearing in steel processing needs out-of-equilibrium (kinetic) chemistry to be computed. One of the major result of the paper is the rediscovery, well-known from an industrial standpoint, that waiting for the mixing time of the reactor between the injection of the addition elements (here Al and Ti) is a good way to get rid of these harmful inclusions, as they eventually dissolve. From a qualitative viewpoint, this is an interesting step in the modeling of the chemistry of a steel ladle. However, there still remain some improvements to be made, both in the modeling (for instance by introducing heterogeneous nucleation) and in the computational efficiency of the coupling. Not previously mentioned is the introduction of the fractal dimension of the 8637
dx.doi.org/10.1021/ie403991e | Ind. Eng. Chem. Res. 2014, 53, 8630−8639
Industrial & Engineering Chemistry Research
Article
oxide aggregate. This is another potentially important parameter32 as some of them are known to be nonspherical (Al2O3 for instance). This may be introduced, for the sake of simplicity as a fixed parameter or, more realistically as an inner parameter evolving through time. However, the latter possibility needs some development, since bivariate population balance equations are more difficult to solve. Recent advances like DQMOM,33 CQMOM,34 or multifluid35 approaches could be used. This will also allow taking into account particle inertia and velocity distribution when they are large enough (more than 1 μm as a rule of thumb). Software computational efficiency can also be improved, as profiling of the code has shown that a bottleneck can be found in the chemistry computation. As computation in the three-dimensional case can take several weeks or months, In Situ Adaptative Tabulation (ISAT36) could be used to enhance the efficiency of the thermo-kinetics code. This will be the goal of future work.
■
Li ,c(t + Δt ) = Li ,c(t )hf
with a factor independent of the abscissa 1/3 ⎛ Δm3,c ⎞ ⎟⎟ hf = ⎜⎜1 + m3,c ⎠ ⎝
■
APPENDIX B fu , i =
2Ωi + ζ 2 2
2Ωi + 2Ωi + ζ
2
, Ωi =
τp, i τL
,ζ=
τλ f TLf
(B.1)
with ⎛ ρp 1 ⎞ 8 R c, i τp, i = ⎜⎜ + ⎟⎟ 2 ⎠ 3 C D, iVR, i ⎝ ρF
Growth Treatment for Elementary Particles
Instead of directly using the diffusion growth rate Gp,i presented in eq 15 or 18, where the weights of the Dirac’s distributions are modified, the abscissas are shifted ensuring growth in a proper manner. Here, we describe more precisely this procedure. First, the oxide volume is connected to the third moment according to the relationship
τLf
τL̅ =
τη
=
(B.2)
2(Reλ + 32) 151/2 ·7
(B.3)
⎛ 2Re ⎞1/2 11 + 7Reλ = ⎜⎜ 1/2 λ ⎟⎟ , C0 = τλ̅ = τη 205 + Reλ ⎝ 15 C0 ⎠ τλ f
(A.1)
where kv is the volume shape factor (equals to π/6 for spherical particles). Then, we note that only the size of the particles (and not their number) is concerned by diffusion growth, so that it is useless to try to modify the weights. The volume variation for node i during a time-step Δt is obtained through
(B.4)
where τp,i, τLf, τλf, and τη are, respectively, the response time of the particle i to the flow, the longitudinal integral time scale, the longitudinal Taylor time scale, and the Kolmogorov time scale. If particles or aggregates are assumed to be in a Stokes flow, then one gets
ΔVi ,p = k v Δ(wi ,pLi ,p3) = k vwi ,p(Li ,p3(t + Δt ) − Li ,p3(t )) (A.2)
CD =
giving finally the new abscissa used to reconstruct the vector of moments for the next time step 1/3 ⎡ ΔVi ,pk v −1 ⎤ ⎥ Li ,p(t + Δt ) = Li ,p(t )⎢1 + ⎢⎣ wi ,pLi3,p(t ) ⎥⎦
(A.5)
Equations A.4 and A.5 supersede eq 16 or 19 in the actual implementation. Note that once the new positions of the Dirac peaks have been computed they are used to develop a new set of moments using quadrature formula 17.
APPENDIX A
Vi ,p = k vwi ,pLi ,p3
(A.4)
24 Rep
(B.5)
which is valid when Ga =
(A.3)
8R c, i 3g |ρF − ρp |ρF μ2
> 10
(B.6)
As a consequence, the product-difference algorithm is not used, because the new treatment of particle diffusion growth only shifts the node abscissas. Therefore, the set of moments obtained (and convected) thereafter is still the set of moments of a real distribution. It may be remarked that even if this new way of treating the particle growth may seem quite simple, it is an important modeling step as it enforces the numerical stability. It has been observed that if weights and abscissas were modified as presented in eq 18, the product-difference algorithm might fail.
where Ga is the Galilei number. This condition translates into
Growth Treatment for Clusters
*E-mail:
[email protected].
R c, i
