Hydrodynamical Simulation of Perspective Installations for Electrometallurgy of Aluminium

Electrometallurgy is based on the electrolysis of fused chemical compounds of the metals. The process of electrolysis is not a simple chemical reaction of mix and match. As electricity is involved in this process, care is taken to understand and set up the apparatus as required. In view of this, basic requirements for theoretical, experimental understanding and following set up procedures are to be studied. The process inherent features, which define the hydrodynamics of the melt are:  high-strength electric current, which generates a large amount of Joule heat; spatial distribution of the volumetric power may be nonuniform  chemical reactions, which change the compound and generate a large amount of the gas bubbles; this may result in gradients of the melt density and current density  bubble rising, which strongly defines the velocities at the melt. The experimental explorations in the electrometallurgy are usually difficult or impossible because of the high temperature and hostile environment. For this reason electrometallurgical processes are the application field of CFD. The aim of this chapter  to demonstrate such application based on well developed and widely used methods of numerical analysis. The object of the development is an electrometallurgy of aluminium. The presented mathematical model of vertical electrode cell takes into account the electric field and current density, chemical composition, heating, natural convection, bubble flows. The integrated model of such kind may be used for the detailed 3D numerical simulations of the new installations.


Introduction
Electrometallurgy is based on the electrolysis of fused chemical compounds of the metals.The process of electrolysis is not a simple chemical reaction of mix and match.As electricity is involved in this process, care is taken to understand and set up the apparatus as required.In view of this, basic requirements for theoretical, experimental understanding and following set up procedures are to be studied.The process inherent features, which define the hydrodynamics of the melt are:  high-strength electric current, which generates a large amount of Joule heat; spatial distribution of the volumetric power may be nonuniform  chemical reactions, which change the compound and generate a large amount of the gas bubbles; this may result in gradients of the melt density and current density  bubble rising, which strongly defines the velocities at the melt.The experimental explorations in the electrometallurgy are usually difficult or impossible because of the high temperature and hostile environment.For this reason electrometallurgical processes are the application field of CFD.The aim of this chapter  to demonstrate such application based on well developed and widely used methods of numerical analysis.The object of the development is an electrometallurgy of aluminium.The presented mathematical model of vertical electrode cell takes into account the electric field and current density, chemical composition, heating, natural convection, bubble flows.The integrated model of such kind may be used for the detailed 3D numerical simulations of the new installations.

Aluminium electrolysis as the task for electrochemical hydrodynamics
It should be noted that for today's aluminium electrolysis cells (AE) CFD modeling unlikely gives essentially new results.One obvious reason is that modern aluminium production technology exists already several decades.To explain the other reasons and the problem statement let's consider the typical multielectrode electrolytic cell for aluminium production.The advantages and disadvantages of the modern AE will be evident from this consideration.The earth's crust contains about 9 percents of aluminium.But because of high chemical activity of this metal alumina (Al 2 O 3 ) is highly stable and so cannot be reduced by conventional reducing agents like coke, carbon monoxide or hydrogen.To detach the metal from the oxygen the sodium aluminum fluoride is used, which is named cryolite (Na 3 AlF 6 ).Actually the industrial electrolyte contains also the additions of aluminum, calcium, magnesium fluorides (AlF 3 , CaF 2 , MgF 2 ), which allow to decrease the melting temperature (Borisoglebsky et al., 1999).For the sake of simplicity this composition in the electrolyte will also be referred as cryolite.Molten aluminium is deposited under a cryolite solution with 3-5% alumina.The metal arises during the well known Hall-Heroult process (Grjotheim et al., 1982 ).The main products of the reactions are the aluminium and oxygen.The latter is generated at anode.Anodes of modern aluminium reduction cells consist mainly of the coal, which is inert with respect to cryolite.But oxygen joints the carbon, and carbon dioxide with small amount of carbon monoxide are the main gaseous reaction products.This results in the ecological problem involved with atmospheric pollution.The second technological problem is involved with the necessity of periodical changing of the spent electrodes because the coal anodes are consumed very quickly.While anode diminishes distance between the cathode and anode plates should be maintained constant.The only way to keep this distance is to make the working space horizontal as it is depicted at left in figure 1. Usually the aspect ratio of the horizontal working space has the order of 1/30.This horizontality causes the next two features.One is in the nonoptimal use of the area of the aluminium plant.The next feature is involved with so called anode effect: the rising of voltage, which can be caused by unsufficient rate of diffusion of alumina to the anode surfaces and also by accumulation of the generated gas in horizontal gap.This effect temporarily blocks the current and may result in stopping of the electrolysis in several cells.The mathematical modeling of modern horizontal AEs is used mainly for the optimization of their thermal conditions and electricity consuming.This is done by solution of heat conduction equation and Laplace's equation for electric potential.In general CFD may be also used for the flow prediction and optimization of horizontal AE functioning.The flows in vertical gaps between the electrodes can be modeled successfully.But in the working space, the high gas content in combination with low velocities results in the conditions, which are beyond the application area of the commonly used multiphase models.The dynamics of the growth and coalescence of the large bubbles in the narrow horizontal interelectrode gap too much depends on the local conditions, which are not well known.The flow in working space is found too complex, and "physical" accuracy of the numerical model is a priory insufficient for the detailed numerical study.Nevertheless, that doesn't mean that we cannot to model anything in the horizontal electrolytic cells.There exists many works on simulation in the electrometallurgy of Aluminium (see, for example, Purdie et al., 1992, Laszlo et al. 2005).One of the earlier works on numerical simulation by Fluent code was issued at 1993 (Purdie et al., 1992).But present-day AEs with horizontal working space are less predictable then those with vertical anode cells.

Perspective designs. Possibilities and aims of their numerical simulation
During the latest decades metallurgists tried to develop permanent anode for the AEs instead of the graphitic.It is also referenced as inert anode (La Camera et al., 1995, Dawless et al., 1999).The possibility of stable geometry of the working space allows vertical configuration of the anode and cathode plates as it is shown on right picture of fig. 1.The above mentioned disadvantages of the existing AEs will be reduced to minimum.We don't concern here with the problem of chemical stability of inert anode, which is very complex, and consider the vertical AE as an object of application of computational fluid dynamics.In assumption of the anode's chemical stability the tasks of the numerical simulations of VAE conditions by the methods of the theory of continuous medium will be the following:  Estimation and optimization of electric current distribution for effective use and uniform spending (through erosion) of the anodes.


Estimation of spatial distribution of alumina concentration for effective functioning of VAE.The problem of alumina feeding of AE exists: alumina is more dense than the melt and tends to sink and be accumulated on the bottom; this requires the modeling of solution of alumina and its consumption in working space.


Removal of volumetric Joule heating from the melt, which would be much more for vertical configuration.This heat is removed out through the installation's walls.The spatial temperature distribution and its extremes are the result of convectional heat removal to the installation's sidewall, and stationary temperature may be too large for electrochemical process.These tasks demand the conjugated solution of the problems from different topics:  Calculation of spatial distribution of electric potential and current density.


Simulation of multicomponent chemistry in the melt (species transport and reactions) or development of well grounded simplifications about the melt composition and its properties.


Heat transfer in solid structures and liquid electrolyte by heat conductivity and convective motion.


Modeling of twophase bubble flow: bubbles affect the spatial distribution of current density and strongly define the flow pattern and velocity.The space and time requirements for industrial applications are:  3D model in realistic geometry.


Modeling of steady states and transition regimes; modeling of unsteady physical processes (like alumina solution in a moving electrolyte).In contrast to horizontal AEs the vertical electrolytic cells weren't widely used for aluminium production, there is no practical experience for them.Although the vertical configuration of the electrode plates is commonly used structure in electrolytic cells, the extreme parameters of aluminium production in VAEs demands intensive study.Due to complexities of the experimenting with new apparatus dealing with hostile environment at high temperatures such apparatus are the object of numerical investigation.The example of mathematical model of such kind, which can be realized in commercial code by means of user's defined functions is presented below.

Mathematical model 2.1 Identification of the phenomena to be modeled and choice of the approaches
Let's briefly outline the physics to be described (or in some cases  to be neglected).The elementary electrolytic cell is sketched in Fig. 2. Due to applied electric field the current passes through the electrodes and cryolitealumina solution.Because of finite conductivity and flatness of the electrodes the electric potential distribution and the current density may be nonuniform (terminal effect, Tobias and Wijsman, 1953).Electric current in the electrolyte is the ion's motion that results in the deposition of aluminium at the cathode and oxygen  at the anode.The oxygen forms bubbles that rise to upper surface where they eliminate from the electrolyte.The bubbles affect the electrolyte conductivity  greater the gas concentration (upper region), less  the electrolyte conductivity.The size of bubbles increases as them float to the surface (due to decreasing of hydrostatic pressure) and due to their coalescence.The bubbleliquid interaction results in volumetric forces (mainly  drag forces) that strongly affect the liquid motion.The volumetric Joule heating (mainly at the electrolyte) results in heat expansion, volumetric buoyancy forces and natural convection of the liquid.Tangent forces from wall friction and bubble motion make the flow turbulent.Since the alumina is spent at the electrodes and is injected at some distant place (feed point) the concentration of the solution is nonuniform that may affect the current density near electrode and, hence  bubble generation.We see that even schematic representation of the alumina process requires conjugated modeling of several different phenomena, which occur in the flow.Each of them is described by its own equation(s).The mathematical model is based on the flow equations for multicomponent and multiphase medium with the specified particular terms, coefficients, boundary conditions, and model assumptions/simplifications.The electrolyte flow is partially a natural convection, and partially is a forced convection due to bubble driven forces.The flow is practically always turbulent.Hence, the equations describing heat transport, electric field, concentrations of electrolyte components, and turbulence are also necessary in the model.The energy equations and electricity equations are solved both in the melt and in a solid domains.Note about electromagnetic forces and possible magnetohydrodynamics wave effects.These effects can really be observed in horizontal AEs at the free interface boundary of liquid aluminium and cryolite.In such conditions (see fig. 1) the conducting medium having large sizes in two dimensions can be waved in electromagnetic field like shallow water (Dupuis & Bojarevics 2005).In the case of vertical AE the accumulated liquid aluminium is practically outside the electric field, and due to electroneutrality of the electrolyte (except the molecular double layer near electrode) the only possible electromagnetic force may be the ponderomotive force in magnetic field of direct current, which may be shown to be sufficiently small and may be taken into account like a source term in the momentum equation.The efficient computation may be achieved using the Eulerian twofluid model with RANS turbulence model for simulation of the two-phase gas-liquid flow.In what follows kε model is used for the turbulence.Electrochemistry effects in the hydrodynamic model are taken into account by addition the convective diffusion equation(s) for scalar function(s) describing the composition of the melt, and source/sink terms for these scalar(s).The effects of magnetic field aren't taken into account in following formulations because they don't dominate in the flow dynamics.The material equations deal with the definite material models.The choice of them postulates definite properties and actually is one of the assumptions of the whole numerical model.

Media: Material models
The following types of materials with corresponding properties should be taken into consideration in simulation of AE:  Solid structure materials  that of cathode, anode, walls and other solid structures of the facility.These materials don't move -the flow equations aren't solved for them.The material properties, which should be specified for solid material are: density, heat capacity, heat conductivity, electroconductivity.  Liquid material  electrolyte, which solidifies at temperature about 970 o C or lower depending on composition.It consists of several constituents, which can be schematized as two basic components: alumina and criolyte.The liquid is referred as a primary carrier phase in two-phase model, and gas -as a dispersed secondary phase, for which the bubble parameters are to be defined.The bubbles are assumed spherical.Both fluid and gas are modeled as incompressible liquid because the flow velocities are not greater than ~1m/s (the speed of sound u s at the gas should be close to that of the air, u sg ~330m/s, u s for the dense electrolyte should be of the order of that of water i.e. u sf ~1000m/s, hence Mach number for carrier flow is of the order of 10  3 10  2 ).Both liquids may be modeled as Newtonian and their material properties are: density, heat capacity, heat conductivity, viscosity, diffusivity (which actually is phenomenological coefficient for considered averaged mixtures), electroconductivity.  (optional) Melting/Solidification.It is usually described on the base of porous media model by introduction of effective sink terms for momentum and enthalpy within a solidusliquidus temperature interval (Fluent Inc., 2005).The alternative way is in introducing the effective capacity and viscosity within a solidusliquidus interval of a phase transfer.In such approaches the solidified material is modeled as liquid, i.e. the flow equations are solved for such material.

Twophase flow equations
The continuity equation for gas phase: where the source term Sg describes the electrochemical gas generation, g  -gas density.
Momentum equation: Here i g -components of gravity vector, , are a lift force, wall force preventing from bubbles reattachment, virtual mass force, and other body forces, respectively.The continuity equation for fluid: (1 ) 0 Momentum equation for fluid: (1 where ff f   , vector gf i R models drag force, O,f i F corresponds to other body forces.The main contribution to the viscosity coefficient is due to turbulent viscosity. The equations ( 1)-( 4) are closed by the definition of drag force There exists a number of the closure equations (Loth, 2008), built for different conditions.The spherical particle (bubble) having the diameter p d and moving relative to carrier flow with velocity gf u , is characterized by the particle Reynolds number and particle relaxation time: The drag force is written through the relative velocity: where interaction coefficient gf K is defined through the drag coefficient D C : The common form of the drag coefficient of the deformable bubble may be taken as linear interpolation between the extreme cases as, in particular at (Filippov, Drobyshevsky et al., 2010): where

C
 are defined as for a solid spherical particle in accordance with eqs.( 9) and ( 10).
The introduction of the expressions for the other forces ( ) acting on the rising bubble is more complex because the influence of these forces is more complex.In particular, the lift force can depend strongly on bubble position and be oppositely directed in close points.All these forces require comprehensive experimental data for particular cases.Because of the natural convection and bubble rising the flow in working space of VAE should be directed upward and spatial distributions of gas volume fraction and vertical component of velocity may be close to that of rising bubbles injected from bottom into isothermal upward water flow in a vertical tube (see Figs. 3-4).The near wall maximum of void fraction is the result of opposite action of lift force and wall force.Essential feature in the case of bubbles in VAE is the absence of the bottom gas injection and generation of gas in the vertical wall.That should result in not such deep minimum of volume fraction near the wall and its smaller values at the opposite side (i.e.near the cathode).Wang, 1987), line  simulation (Filippov, Drobyshevsky et al., 2010).
Note that the attempts to introduce the lift force and wall force in the full numerical model of VAE were based on the methodology outlined in (Filippov et al., 2010).As a whole these attempts didn't give a stable and robust algorithm in full model of VAE, although the separate tests on bubble motion with such model implemented in OpenFOAM  code were successfully performed (Filippov, Alipchenkov et al., 2010).Since the model of VAE is seemed to be overloaded by another physics, the simple phenomenology with effective horizontal force was introduced, which should integrally simulate the effect of . The expression for such effective force and important question about the bubble size and polydispersity of secondary phase are discussed lower.
The boundary conditions for carrier flow are the common wall conditions for channel flows.
No free surfaces are considered (that actually corresponds to crust formation at the upper boundary of electrolyte).The conditions for the secondary gas phase are stated as: the flux condition at the anode surface, in which gas is generated in accordance with the Faraday's low; the elimination conditions at the upper boundary:  =0; -slip conditions at the other boundaries.Note that in the case of crust of solid electrolyte situated atop the moving liquid the gas removal from the upper boundary of vertical working space meets with difficulties that may result in more wide area of bubbles expansion in the flow.In that case the simple condition of gas elimination  =0 may be substituted by more appropriate in accordance with the available empirical data.

Turbulence equations
For simulation of turbulent flows in multiple narrow gaps of working space of VAE only RANS models can be effectively used.The terms describing the generation of turbulence by buoyancy and bubble motion are better developed and tested for twoequation k  model which is used by authors in all simulations.The equations for TKE and its dissipation rate (Pope,2000, Wilcox, 1994) in common form taken from Fluent's theory guide as follows: The turbulent viscosity where C  is model constant.The density in these equations can be the density of mixture, the density of a carrier flow (dispersed turbulent model) or the density of a particular phase (if they are simulated as symmetric), dependently on the kind of the taken multiphase turbulent model.Additional terms are introduced in modifications of k   model for multiphase flows, which describe the effect of bubble presence, in particular, pseudo-turbulence induced by bubble wakes (van Wijngaarden, 1976).The question about effect of a turbulence even in the dispersed multiphase flows is even more complex matter than the above mentioned problem of modeling the body bubble forces.In absence of reliable experimental data the default data of the used commercial code are taken as the starting point.

Enthalpy equations
The enthalpy equations for primary and secondary phases include the balance of transport, diffusion, heat sources and (optionally) viscous heating: The enthalpy s h of s-th phase (where s=f, g) is: where T s is the temperature, f V is viscous heating, fg Q is interphase heat exchange, f S  other heat sources, cT () is constant pressure heat capacity.Heat conductivity coefficient s  takes into account turbulent heat transfer.The interphase heat exchange term fg Q depends on the phase temperature difference.It is written through the heat transfer coefficient fg H : which for dispersed flows usually is calculated with the aid of the well known Ranz-Marshall correlation.
The temperature of the generated gas in AE should be equal to that of liquid and there is no reasons for these temperatures to be different in an incompressible liquid.Since there is no significant evaporation at the electrolyte, there is no significant interphase heat-and mass exchange between the primary and secondary phases.Therefore, interfacial heat-mass transfer doesn't play significant role in AEs.Practically all generated gas erases from the system.It is interesting to compare the heat removal due to gas outlet with the total heat generation.The mass generated from 1 m 2 of the electrode plate is given by Faraday' lows: where and j is normal component of current density.The temperature of the generated gas is equal to that of its environment, hence, the rate of gas enthalpy generation can be estimated as: .We see that for rough simulations in many cases the heat removal by the gas may be neglected and energy equation for the secondary phase -be omitted.Boundary conditions for the enthalpy equations, which for incompressible flows are often solved for the temperatures, are the usual BCs for heat conduction equation.These BCs are set the same for both phases.On a coupled boundary between the solid and liquid materials the conjugation condition of continuous normal heat flux is stated.The same condition is stated on a boundary of different materials for electric potential equation considered below.

Electric field and current density
In conducting material the distribution of electric potential  is related with current density by Ohm's low: Due to electro-neutrality of the all media div =0 j or This equation defines the spatial distribution of electric potential.The current density is used for calculation of the source term of volumetric Joule heating in the energy equation Eq. ( 15) for liquid phase.
The gas bubbles occupy the significant part of the volume that reduces correspondingly the cross-section of the interelectrode space, i.e. partially screens electric current.In the case of relatively small gas volume fractions (approximately  <0.2) the decrease of conductivity  due to decreasing of effective conducting section depends on the bubble concentration through the Bruggeman's equation (Bruggeman, 1935): where 0  is the conductivity of homogeneous electrolyte.
When gas volume fraction is close to 0.5 or larger the flow is poorly predictable.In such a case the equation ( 25) or similar to it may be used like phenomenology that gives right result in extreme case  =1.The boundary conditions for potential equation ( 24) are specified in the places of current input/output, i.e. in a poles of the electric circuits.The values of current density at one pole P  and electric potential at other pole P + are specified: On the other boundaries  insulation condition (default BC): The overpotential effect (Newman & Thomas-Alyea, 2004) in the case of small potential drop along the electrodes may be taken uniform for all anode boundary and thus may be easily modeled.In general, the overpotential have to be determined for every definite installation with use of the detailed data on the electric resistance.
In calculation sequence of electrical part of full model Eq. ( 24) is solved first, then the current is defined from Eq. ( 23) and source term (15) for the enthalpy equation as well.Since the current density depends on spatial distribution x ()  of secondary phase, the common integration scheme for the full system with sequential solution of flow equations, energy equation and equations for field, concentration etc. will be explicit.But owing to slow change of x ()  obtaining of stable numerical solution is possible.

Electrolyte composition
The effect of nonuniform distribution of alumina concentration is of importance for consideration of the influence of the local deficiency of alumina on the current density and aluminium/gas generation.If the local concentration is greater than the definite minimum value (about 2%) the process is ruled by the electric potential and conductivity.The variations of electrolyte composition weakly influence the flow in comparison with the bubble motion (see below).But variations of the concentration of alumina can violate the electrolysis (the value of alumina concentration  <2% results in anode effect) and thus maintaining of  is of primary importance.From such point of view the modeling of electrolyte composition is important but somewhat separate problem, which may be solved with different levels of detailing.The zero level corresponds to neglect of electrolyte composition.A such model was outlined above.The following level of detailing is twocomponent model of composition and we will try to show that the possibility of such simplification really exists.
Actually the cryolite and alumina dissociate during their solution and form several charged ions.But this microscopic structure isn't exhibited macroscopically.The average distribution of concentration of the ions smoothly varies in space, except the very narrow double layer at anode surface, which cannot be modeled macroscopically (Borisoglebsky et al., 1999).In some extent the local stoichiometric composition depends on the dissolved alumina concentration only.It is important that having at the input of Hall-Heroult process one mole of Al 2 O 3 we obtain two moles of Al and one and half moles of oxygen O 2 (or 1.5CO 2 ), i.e. no new chemical compounds arises in the solution during the process of aluminium reduction.
All that allows to consider the alumina electrolyte as the neutral idealized compound of alumina and cryolite only.The properties of such system (in particular, the density) will depend on its composition  the empirical data should be used.Note that electrolyte of industrial AE can contain up to 10% of additions (Borisoglebsky et al., 1999).The generalizations of the idealized case of twocomponent mixture that is considered below depend on particular application of the numerical model.For example, we don't concern here the problem of alumina dissolution since this complex process depends on the details of technology and apparatus of alumina feeding.
The two components of the electrolyte are described by the alumina mass concentration  .Correspondingly, the criolyte concentration will be 1   .The transport equation for  is the equation of convective diffusion (Newman & Thomas-Alyea, 2004).In more or less general form this equation can be written as follows Actually this is a balance equation for the alumina/cryolite masses, which are together form the primary phase.Since the gas density is negligible the density  in this equation may be the mean density of gas-liquid flow, which carries the concentration parameter  .The relation of actual density of electrolyte (primary phase) f  with  should include the temperature and is taken from the empirical data.Note that here the diffusion coefficient  Abramovitz, 1981).The places of sources and sinks of alumina mass i.e. the sources and sinks of concentration for Eq. ( 29) are the alumina feed places and vertical working space of AE respectively.If we don't model the chemistry in electrolyte, the place of the sink of the averaged component (alumina) is not exactly positioned in the working space, and it is naturally to set this sink at the electrode surface.Both cathode and anode areas may be taken as the boundaries of concentration sinks, since due to intensive bubble motion the electrolyte composition should be effectively smoothed.The mass flux boundary condition describing the sink of alumina during the electrolysis at the boundary out  is ruled by Faraday's lows, analogously to that of gas generation (18)(19).The generation/elimination of xt (,)  at every from the two electrode surfaces is defined as: where function out mx t (,)  is the given spatial distribution of alumina feed sources or is the sink at the electrode, which is calculated from the Faraday's lows using the calculated normal component of current density: The source term of alumina feed may be specified both through the volumetric source S A and the analogous flux boundary condition in a some part of the AE walls.
The velocity field u i (x) in the convective term is known from the flow equation.The transport of species in the liquid should be referred to the velocity f i ux () () of a primary phase, which carries the species.Since the gas density is small, the velocity f i ux () () is close to the centre of mass velocity in a mixture.

Further simplification, adaptation and verification of the integrated model
The equations presented above are well known and widely used.For real application of them to VAE the source/sink terms, coefficients, and boundary conditions should be specified.For the sake of effective solution of these conjugative equations some simplifications should be done.Some of them are discussed below.The possibilities of the simplifications are to be checked by verification and validation of the numerical models.

Some features of heatmass transfer in horizontal and vertical AEs
The thermal regime of AE installations can be characterized by the averaged volumetric heat generation per unit volume fAA p Ap SWV /  , where Ap W is the total heat generation in the cell, Ap V is the estimated heated volume.As it was mentioned, vertical AEs may have the parameter fA S some times greater than that of horizontal AEs.The reason is in their geometry as it can be seen from Fig. 1: considerable volume of the horizontal AE is occupied by the electrodes and electrolyte lying out of working space.The heat scattering in horizontal AEs is successfully calculated with use of common heat conductivity equation, the areas of the electrolyte flow are taken into account in such a simulations through effective heat conductivity coefficients.Since the working space in horizontal AE is situated downward, the alumina granules being introduced in the electrolyte dissolute near the bottom and reach the working space without considerable problems.
In vertical AEs the fA S parameter is larger since the working space occupies relatively larger volume.The dominant physical process in heat removal is the convection.Let's compare the mass flow rates resulting from the natural convection and from bubble driven forces.The cause of the motion in both cases may be characterized as an averaged density gradient in mixture.In case of natural convection in a single-phase liquid the gradient arises due to thermal expansion and nonuniform composition.The latter factor gives the relative density variation  less than 2 10  (taking into account that variations of alumina concentration in working space are normally not much greater than 2%).Assumed temperature drops in working space shouldn't be greater than T~100K  , heat expansion The averaged buoyancy force involved with presence of the gas bubbles is G fg 0 ~ .
Usually gas volume fraction  has an order of 0.1, i.e. in that case we have GT ff  .Hence, in working space bubble driven (drag) force should dominate and mainly defines the flow at the installation.This allows some roughness in determination of the density through the temperature and composition, but requires to pay larger attention to modeling of interfacial momentum exchange, namely to drag force.

Simplifications in multiphase model
The gas flow at the lowest part of vertical electrolytic cell consists of separate bubbles, gas volume fraction through the width of the gap is relatively small.The problem of description of a dispersed bubble flows in the vertical gap is satisfactory studied to this day for monodisperse bubble flow with 0.1 0.2  .There exists several models for polydisperse flows (see, for example (Krepper et al., 2007), (Silva & Lage, 2011), (Guan Heng Yeoh & Jiyuan Tu, 2010).Since the gas volume fraction in VAE may amount to the values of some tens percents, the bubble flow in the vertical working space may be really polydisperse due to various departure diameters and coalescence as it is sketched in Fig. 5.The implementation of the model of polydispersity in the whole numerical model of VAE will result in the increasing of the number of equations, because several additional equations will be added, which describe the introduced size ranked groups of bubbles (usually  34 or more).We should add the model of coalescence, probability density function (PDF) for bubble diameters.These equations and models require the experimental data, which for bubbles in cryolite are rather poor.The result of the such sophistication will be in an improvement of the model of dispersed bubble flow for the range of gas volume fraction of the order of 0.1-0.2.For higher volume fraction the applied MCFD methods will be essentially phenomenological, and here we have the common dilemma: before the implementation of a new improved approach the researcher should decide whether the improvement be really reached and what is the calculation cost of the work.Because of complexity of the full numerical model of VAEs and the lack of the experimental data the models of polydispersity were not considered for them b y a u t h o r s .I n a l l simulations the bubble size was taken as an averaged value, which may be varied during serial quality assurance calculations.The effects of the variable bubble size in a working space may be investigated numerically as a separate effects and be taken into account through the model coefficients.Another important question to MCFD model concerns with the above mentioned body bubble forces.They influence the distribution of  across the vertical gap, that defines the cross profile of vertical component of a velocity and in somewhat may effect the mass flow rate in the interelectrode gap.There exists several forms of equations for the lift force (see Antal et al., 1991, Legendre & Magnaudet, 1998), which were validated and used by many authors.The testing and identification of the models of the lift force and wall force (Antal et al., 1991, Lopez de Bertodano, 1994, Troshko & Hassan, 2001), which prevents from the bubbles reattachment, is usually performed in the bubble column with downward bubble injection.Airwater flows are commonly considered in such experiments.The analogous experiments for wall generation of the bubbles in fused cryolite are unknown to authors.The experience shows that implementation of new correlations in multiphase models always requires some fitting of the coefficients during the model verification.Since the using of lift force model in common form often reduces the stability of numerical calculations, more simple form of horizontal component force was taken for simulation of bubble motion across the gap.The robustness of simulations is an argument to develop and use the simplified "engineering level" models, which have to be tested and studied for prototypical conditions.Another argument to use the simplified form for simulation of the bubble body force is in the above mentioned relatively small established range of the validity of disperse models.The experiments with bubble flows are conducted usually for airwater mixtures in conditions of clean and even walls of bubble columns and low gas volume fractions.In conditions of the aluminium reduction cell the wall roughness, the bubble size and gas volume fraction are usually relatively high, and the model validated on the water experiments may be incorrect in the conditions of VAE working space.In the experiments with vertical wall gas generations (see, for example, Cheung & Epstein, 1987, Ziegler & Evans, 1986) the average thickness of the near-wall bubble layer increases with the height.The proposed phenomenology takes into account only this fact.An effective body force is introduced for the gas phase, it acts in the secondary phase normally to wall, in which gas is generated (Fig. 6).The simplified phenomenological forms of expressions for effective body force, acting on the bubbles in a narrow vertical channel and causing their horizontal drift, were introduced by different authors in different forms.This approach was previously used to model vertical electrolysers in a 2D formulation, in particular, using discrete particle model of Fluent code (Mandin et al., 2006, Bech et al., 2003), where the horizontal drift of bubbles rising along the vertical wall was caused by additional effective force.Bech et al., (2003) considered such horizontal force as being a function of the distance to the vertical wall and expressed it via the force potential as where A and n are constants and gg rd 2  is the bubble radius.The effective body force acting on the gas phase is as follows: Mandin et al., (2006) assumed the horizontal force to be independent on the distance across the gap and to be such that its ratio to the drag force was equal to the assumed slope of the bubble trajectory.A significant disadvantage of expressions ( 32) and ( 33) for the effective horizontal force is in their dependence on the domain geometry and the independence of the local gas concentration that can restrict applicability of this model for high gas content and leads to nonphysical effects.The effective force is better to be expressed only through local characteristics of a two-phase flow.During introducing this force it is natural to assume that the interaction between bubbles and turbulent flow, which is caused partially by other bubbles depends mainly on bubble's spacing and their average size.A simple assumption is that this force depends explicitly only on the gas volume fraction, while the dependence on the liquid velocity, viscosity, etc. can later be introduced via the corresponding multipliers.Then, it is also natural to introduce a gradient of  into the expression for the effective force, since, as it may be assumed, the gradient should be involved with the bubble repulsion direction in a case of high gas content.A simple expression for such an effective force is as follows: where A and n are some constants.This formula takes into account the direction of the gradient of the gas volume fraction  .In the region of homogeneous  (which usually corresponds to the uniform distribution of bubble velocities), the effective force vanishes.At the edge of the bubbly region, the effective force (34) acting on the gas phase decreases with the concentration.The particle relaxation time ( p τ ) and corresponding length ( p l ) are small.
For a bubble diameter of d = 3 mm, we obtain Therefore, when a bubble leaves the two-phase zone, its motion related to the effective body force is rapidly terminated.In the region of small  , the motion of a gas phase is determined only by their buoyancy and the carrier flow.The vertical component of gradient of  may be or may be not taken into account.

Numerical implementation and solver options
The modern CFD codes at parallel computers can solve all the equations presented above, on the meshes of the order of ten million cells even for transient problems.The minimal set of the equations includes usually the flow equations (1-4) with equations of turbulence model and energy equation (14-15) for all phases.The additional equations of convective diffusion type describing definite physical scalar values (such as concentration equation ( 29) or electric potential equation ( 24)) may be added together with their coefficients (user's defined scalar, UDS in Fluent code).Handling with coefficients of all equations is possible due to user's defined functions (UDF).
The particular attention should be paid to choice of solver's options since the numerical diffusion may distort the results.The problem of non-uniqueness of the obtained numerical solutions also exists.These issues lead authors to the following preferences in CFD code options which seem to be more or less general to all complex simulations of VAES:  steady problems is commonly solved as transient by relaxation to steady state  pressure-velocity coupling procedure uses SIMPLE (Versteeg & Malalasekera, 1995) algorithm for steady problems and PISO or Coupled (simultaneous solution of four flow equations for pressure and velocity components) solvers for really transient statements.The Coupled solver is preferable because of its better numerical stability  second order spatial approximation in all solved PDEs  realizable or RNG versions of k   turbulence model.

Verification
The verification of the model was done for the separate phenomena, the modeling of which isn't included in standard possibilities of the used CFD code.For the models in consideration such verification set should include the tests concerning with:  electroconductivity problems  natural convection of heat generating liquid  rising of spherical bubbles generated at the vertical walls  convective diffusion.It is essential that the verification of twophase procedures should be done in prototypical conditions, close to that of AEs as it was discussed above.To this day the following problems were passed:  electroconductivity problems: a) simple simulative problem; b) the problem of terminal effect (Filippov, Korotkin, Urazov et al. 2008)  natural convection of molten salt with solidification at cooled boundary (Filippov et al. 2009)  conjugated simulation of effect of rising bubbles on the electrical conductivity (Filippov, Korotkin, Kanaev et al. 2008).The verification and identification of the coefficients of the described model of effective body forces is not yet satisfactory performed because of luck of data on bubbles motion in vertical generating channels.One performed test concerned with CFD calculations is briefly outlined in what follows.

Simulation of elementary electrolytic cell
The elementary electrolytic cell (see Fig. 7) was considered in 2D geometry: cross section by vertical plane, which is normal to the electrode plates was considered.Thermal boundary condition of third kind   Figure 8 shows the vertical profiles of the current density, which were obtained numerically for problems (a) and (b) in comparison with the corresponding results of estimations using formulae presented in (Filippov, Korotkin, Urazov et al. 2008).In analytical solution the plates occupy all the height of the cell.The potential difference in analytical calculation was taken at a half-height of the electrode plates.Note that the partial screening of potential in the upper part of the working space promotes leveling of the current distribution along the vertical axis.The effect of a gas phase in the flow is illustrated by Fig. 9, which shows the path lines for cases (a) and (b).Gas phase distribution is shown in Figs. 10 and 11.The introduction of the gas phase increases the maximum flow velocity (see Table 1) and changes the flow pattern.Table 1 presents some integral characteristics of the flow regime in the two cases under consideration.Here, the average gas volume fraction m  was determined as the ratio of the volume occupied by the gas phase to the total volume of the calculation domain.The data of Table 1 show that both the maximum temperature and the temperature drop through the domain in case (b) are smaller that in case (a).This is related to a more intense motion of liquid and more intense heat exchange in the presence of the gas phase.Note that the temperature variation through the calculation domain is relatively small.


The physical state of the media at the aluminium electrolysis may be described in terms of two-phase hydrodynamics, heat transfer, electric current, and convective diffusion with definite assumptions on chemical compositions.


The uncertainties of models, material properties, boundary conditions and others limit the multiphase model's accuracy, therefore, in the integrated model of VAE including simulation of the bubble motion in cryolite, only the "engineering" level accuracy can be achieved.This means on the one hand that we shouldn't try to achieve the "excellent agreement" with the results of the precise airwater experiments.On the other hand this allows some simplifications in the integrated model for the sake of its robustness and effectivity.From such point of view - The integrated mathematical model of the processes in VAE was developed and realized on the base of Fluent code.Some reasons and the ways of simplifications were discussed in the chapter.


The model was verified for the separate phenomena in simple geometrical configurations.The set of the tests includes heat transfer, natural convection, electric current distribution.The verification and identification of the introduced simple correlations for body forces defining the bubble's motion is required.


Because of the simplifications made during the development, the validation of the full model on the integrated experiments with AE installations is required.


The simulations carried out with the built integrated model demonstrated the robustness and efficiency of the calculations.
coefficient of the spherical particle having almost zero Weber number, bubbles in contaminated liquid.The value of We is limited by We<12 that corresponds to the fragmentation of the bubbles.One of the most commonly used relation for between the solution for spherical segment and Stokes law for a sphere.For a contaminated liquid (that is the case for the electrolyte) coefficients

MD
has the dimensionality of dynamical viscosity, Pas.It is equal to common diffusion coefficient D, multiplied by mean density.It includes the laminar and turbulent diffusion coefficients:
the top boundaries of the electrodes.Here the heat exchange coefficient HW m K 200 /( )  corresponded to the estimated heat removal.The bottom and vertical lateral boundaries were adiabatic that should be close to the condition in the elementary cell in the middle section of VAE.The potential difference was equal to 1V.The results of two calculations are considered below, which were performed for the cases: (a) the absence of gas generation, and (b) the presence of gas generation.

Fig. 10 .
Fig. 10.Gas volume fraction The source terms Tb GG , describe the production of k and  due to velocity gradient and buoyancy, S k , S e stand for other possible production processes, 20)where reference temperature T ref corresponds to room temperature, T ref =300K.Let the distance between the electrode plates is b.The Joule heat generation in working space per 1m 2 of electrode plates is www.intechopen.com