外文原文.pdf

AbstractqSponsors: Alcan-Pechiney Company and Swiss National Science Foundation; Grant No. 200020-101391.*Corresponding author. Tel.: +41 22 379 23 66; fax: +41 22 379 22 05.E-mail addresses: yasser.safaepfl.ch, yasser.safaobs.unige.ch (Y. Safa).Available online at Applied Mathematical Modelling 33 (2009) 14791492/locate/apm0307-904X/$ - see front matter C211 2008 Elsevier Inc. All rights reserved.A phase changing problem motivated by the modelling of thermal problem coupled with magnetohydro-dynamic eects in a reduction cell is studied. In a smelting cell operating with HallHeroult process, the metalpart is produced by the electrolysis of aluminium oxide dissolved in a bath based on molten cryolite 1. Var-ious phenomena take place in such a cell for which a transverse section is schematically pictured in Fig. 1.Running from the anodes through liquid aluminium and collector bars, the steady electric current spreadsin the electrolytic bath. The important magnetic field generated by the currents carried to the alignment ofcells, coupled with the currents running through the cells themselves gives rise to a field of Laplace forceswhich maintains a motion within these two conducting liquids. A magnetohydrodynamic interaction takesplace in the cell. In the other hand a heating source is produced by the Joule eect due to the electric resistivityof the bath.A system of partial dierential equations describing the thermal behavior of aluminium cell coupled with magnetohy-drodynamic eects is numerically solved. The thermal model is considered as a two-phases Stefan problem which consistsof a non-linear convectiondiusion heat equation with Joule eect as a source. The magnetohydrodynamic fields are gov-erned by NavierStokes and by static Maxwell equations. A pseudo-evolutionary scheme (Cherno) is used to obtain thestationary solution giving the temperature and the frozen layer profile for the simulation of the ledges in the cell. A numer-ical approximation using a finite element method is formulated to obtain the fluid velocity, electrical potential, magneticinduction and temperature. An iterative algorithm and 3-D numerical results are presented.C211 2008 Elsevier Inc. All rights reserved.Keywords: Aluminium electrolysis; Cherno scheme; Heat equation; Magnetohydrodynamics; Ledge; Solidification1. IntroductionNumerical simulation of thermal problems coupledwith magnetohydrodynamic eects in aluminium cellqY. Safa*, M. Flueck, J. RappazInstitute of Analysis and Scientific Computing, Ecole Polytechnique Federale de Lausanne, Station 8, 1015 Lausanne, SwitzerlandReceived 27 December 2006; received in revised form 4 February 2008; accepted 8 February 2008Available online 29 February 2008doi:10.1016/j.apm.2008.02.011Electrolyte Anode BlocksFig. 1. Transverse cross section of aluminium reduction cell.1480 Y. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492On the wall of the cell, a solidified bath layer, the so-called ledge is created. These ledges protect the cellsidewall from corrosive electrolytic bath and reduce the heat loss from the cell (see 2 page 23). Moreover,its profile strongly influences the magnetohydrodynamic stability causing oscillations of the aluminiumbathinterface which could decrease the current eciency. Consequently an optimal ledge profile is one of the objec-tives of cell sidewall design.The thermal solidification problem in smelting cell has been treated by several authors 35. As far as weare aware, this problem has never been considered when coupled with the magnetohydrodynamic fields. Theaim of this paper is to deal with such fields interaction. Let us mention that the details on this problem can befound in Safas thesis 6.Mathematically, the problem is to solve a coupled system of partial dierential equations consisting of theheat equation with Joule eect as a source, Maxwell law equations with electrical conductivity as a function oftemperature and NavierStokes equations. The interface between aluminium and bath is an unknown. Theledge is considered as electrical insulator, the thermal model is a stationary two-phases Stefan problem.The outline of this paper is as follow: in Section 2 we introduce the physical model, the algorithm is presentedAluminium Cathode LiningFrozen ledge Frozen ledgein Section 3 and we give the numerical results in Section 4.2. The modelIn order to introduce the model we first describe some geometrical and physical quantities.2.1. General descriptionsThe geometry is schematically defined by Fig. 1. We introduce the following notations:C15 X X1X2: fluids and solid ledge,C15 N N1N2: electrodes,C15 K XN: domain representing the celland we define the interfaces:C15 C oX1 oX2: free interface between aluminium and bath, which is an unknown,C15 Ri oK oNi; i 1;2,C15 R R1R2: outer boundary of the electrodes.Y. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492 1481C15 Cp: specific heat,C15 : latent heat.2.2. Physical assumptionsThe model leans on the following basic hypotheses:1. The fluids are immiscible, incompressible and Newtonian.2. In each domain Xi, i = 1, 2, the fluids are governed by the stationary NavierStokes equations.3. The electromagnetic fields satisfy the stationary Maxwells equations, Ohms law is moreover supposed tobe valid in all the cell K.4. The electrical current density outside the cell is given (current in the collector bars).5. The electrical conductivity r is function of temperature h in the fluids and electrodes parts.6. The viscosity g, the density q and the specific heat Cpare temperature independent.7. The volumes of the domains X1and X2have given values (mass conservation).8. The only heat source is produced by the Joule eect due to the current crossing the cell.9. Eects of chemical reactions 7, Marangoni eect 8,9, surface tension as well as the presence of gas floware neglected.2.3. The hydrodynamic problemIn this part we consider the temperature field h and the electromagnetic fields j and b as known. We chooseto represent the unknown interface between aluminium and bath by a parametrization of the formCC22hx;y;z : z C22hx;y;x;y2DC138, where D is usually a rectangle corresponding to the parametrizationof aluminiumcathode interface. We denote the dependence of X1;X2and C with respect toC22h by usingThe unknown physical fields with which we shall deal are listed as follows:Hydrodynamic fields:C15 u: velocity field in Xi;i 1;2; (u 0 in solid ledges),C15 p: pressure.Electromagnetic fields:C15 b: magnetic induction field,C15 e: electric field,C15 j: electric current density.Thermal fields:C15 H: enthalpy,C15 h: temperature.The material properties are defined asC15 q: mass density,C15 rband r: electrical conductivity in and, respectively, outside the bath,C15 g: viscosity of the fluids,C15 l0: magnetic permeability of the void,C15 k: thermal conductivity,Xi XiC22h; i 1;2; C CC22h:hx;ydxdy V1; where V1is the volume of aluminium:1C22 C22Here3thosethe fluids.fieldsThe fluidC22 C22In orderinvolving a penalization tool. The velocity and the pressure will then be defined in both liquids and solids. Wefunction K is given by Carman Kozeny” law:the Darcy1482 Y. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492When fs! 1, we get Kfs!1and then u 0 in the solid zone.law:rp qgzC0Ku j b:If only liquid phase is present we have K 0 and the above equation reduces to the usual NavierStokes equa-tion. Inside the mushy zone K may be very large, compared to the other terms, and the above equation mimicsqu;ru C0div2lDuC0p qgzIKu j b in X1C22hX2C22h: 7where P is the mean pore size and C is a constant obtained experimentally (see 10). Eq. (1) may then be mod-ified toKfslCf2sP21C0 fs3;add to NavierStokes equation the term Kfsu;fsis the solid fraction which is a function of temperature. Thepart of Xih i = 1, 2 is only a subdomain of the domain Xih delimited by the front of solidification.to solve the hydrodynamic problem in a fixed domain Xi, we use the method of fictitious domain”u 0onoX; 4uC138CC22h 0; 5C0pI2lDunC138CC22h 0: 6(.,.) is the usual scalar product on R . Eqs. (1)(3) correspond to 1st and 2nd assumptions. We completeequations by introducing the conditions on the boundaries of the domains X1C22h and X2C22h containingFor any field w, wC138CC22hdenotes the jump of w across CC22h, i.e. wC138CC22h wbathC0 waluminium. For theu and p we havewithDu12ru ruT; I dij i;j 1;2;3:qu;ru C0div2lDuC0p qgzIj b in X1hX2h; 1divu 0inX1C22hX2C22h; 2u;rz C0C22h 0onCC22h; 3We consider the following standard set of equations for hydrodynamic fields:n krz C0C22hkrz C0C22h:DThe unit normal to CC22h pointing into X2C22h is given byZC22From assumption (vii) we get the following relation:u;n0onCC22h; 102.4. Therot bZC0r 0onoKnR; 162.5. TheY. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492 1483We consider as known the hydrodynamic field u and the electromagnetic field j. The steady solution we arelookingC0ro/on j0on R2; 17/ 0onR1; 18j rC0r/ u b in K; 19bxl04pZKjyx C0 ykx C0 yk3dy b0x8x 2 K: 20thermal problemon4pK kx C0 yk3where b0is some magnetic induction field due to the electric currents which flows outside the cell.The electromagnetic problem PEMis then formulated as following: for given u andC22h, find /;b and j suchthatdivrC0r/ u b 0inK; 15o/bxl0jyx C0 ydy b0x8x 2 K;where j0is the given current density on the outer boundary of the anode R2. Notice that magnetic inductionb is obtained as a function of electrical current j by using BiotSavart relation:C0ro/on 0onoKnR;C0ro/on j0on R2;/ 0onR1;We denote byoonthe operator n;r, here n is the outer unit normal on oK.We introduce the following boundary conditions concerning electric potential /: l0j and Ohms law: j rC0r/ u b in K, we then obtain the electric conservation law given bydivC0rr/ u b 0inK:We consider the velocity field u as well as the temperature h are known. From the Faradays law we have rote 0, the electric field is then given by e C0r/, where / is the electric potential field computed in K. We stilldenote by u the continuous extension of the velocity by zero in K, taking into account Amperes law:u 0onoX; 11uC138CC22h 0; 12C0pI2lDunC138CC22h 0; 13ZDC22hx;ydxdy V1: 14electromagnetic problemWe finally obtain the hydrodynamic problem PHD: for given j;b and h, find u, p andC22h such thatqu;ru C0div2lDuC0p qgzIKu j b in X1C22hX2C22h; 8divu 0inX1C22hX2C22h; 9for will be here obtained as the limiting case of a time dependent heat equation.Inthe locatso needshighlyan encan beHh qCpsds 1C0 fsh: 21Therespondingwhich is a non-linear convectiondiusion system. The term u;rh denotes the scalar product of u withto carefemployed to solve our phase change problem.whertransfer,wherAnThe1484 Y. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492oHotC0divkhrhqCpu;rhS in QT; 27h bH in QT; 28khohon ahhaC0h on RT; 29QT KC2C1380;T and RT oKC2C1380;T:thermal problem PThtakes the form: for given u;C22h and j, find h and H such thatFor a given scalar value T, which will represent the integration time, we denote:e c1;c2and c3are positive values provided by experimental estimation.initial condition on enthalpy Hx;0H0on K is assumed.transfer is due to convection and radiation. The radiation is implicitly taken into account by using:a ahc1c2hC0 c3 W=m2 C14Ckon ahaC0h on oK; 26eohonis the derivative in the direction of the outward unit normal on oK;a is the coecient of thermalwhich may depends on both space and temperature, and hais the temperature outside K. The heatThe temperature h is subject to the Robin boundary condition:ohully track the location of solidliquid interface is removed and standard numerical technique can berh, S is the heat source provided by Joule eect only. It takes the formS rkr/k2: 25The advantage of this temperatureenthalpy formulation, taken in distributional sense, is that the necessityh bH; 24lem as a Stefan problem in temperature and enthalpy under the formoHotC0divkhrhqCpu;rhS; 23function bH is computed by mathematical processing (interpolation) in the list of h;H values cor-to the inverse relation H bC01h given in Eq. (21). With this relation we can formulate the prob-Since the enthalpy Hh is a monotonic function we can introduce the function b defined by the relationbHh h: 220Zhthis subsection, we thus introduce the evolutionary thermal model. In our convectiondiusion problemion of the front of solidification (interface separating ledge and liquid bath) is not known a priori andto be determined as part of the solution. Such problems, widely referred as Stefan problems”, arenon-linear. In order to overcome diculty related to the non-linearity of the Stefan interface condition,thalpy function is defined, it represents the total heat content per unit volume of material. The enthalpyexpressed in terms of the temperature, the latent heat and of the solid fraction fs, namely:H H0in K;for t 0: 302.6. TheWewe haveTheH anfs f3. TheThegeometry without normal force equilibrium condition on the interface and then we update the interface posi-Y. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492 1485u;n0; on Ch;C0pI2lDun;tC138CC22h 0; 8 t tangential vector on CC22h;the problem is then easily formulated on a weak formulation.C15 Step 2: we update the position of interface in order to verify the equilibrium of normal forces on the inter-face C, we chooseC22h :C22h dC22h with:dC22h C0C0pI2lDun;nC138CC22h Cstej b;ezC0qgC138CC22h:Here we denote by ezthe unit vector of Oz axis and by Cste the constant obtained from the condition:ZZDdC22hx;ydxdy 0:An iterative scheme is used to compute um1;pm1;Cstem1;wm1andC22hm1as functions of the valuesobtained from the previous iteration m.We set w C0pI2lDun;nC138CC22hand fzj b;ezC0qg and we define the solving step SmHDbyqum;rum1C0div2lDum1C0pm1qgzIKum1 jm bmin X1C22hmX2C22hm; 31divum10inX1C22hmX2C22hm; 32um1;n0onCC22hm; 33C0pm1I2lDum1n;tC138CC22hm 0 8 t tangent on CC22hm; 34wm1 C0pm1I2lDum1n;nC138CC22hm; 35C22hm1C22hmC0wm1 Cstem1fmzC138CC22hm; 36C22m1C22m1: we solve the hydrodynamic problem for the given geometry CC22h and by taking into account theconditions:C22tion by using the non-equilibrium normal force. The solving deals with the alternative application of the twofollowing steps:C15 Stepinterfacehydrodynamic problem is iteratively solved. In each solving step, we first solve the problem in a fixed3.1. Computation of the hydrodynamic fieldspseudo-evolutive” algorithm involving a space discretization by finite element method is applied forlving of the three coupled problems.The numerical solution of the mathematical above problems is based on an iterative procedure in which wecarry out alternatively the computation of the three types of unknowns: hydrodynamic HD, electromagneticEM and thermic Th. In this section we present the iterative schemes for the problems PHD;PEMand PTh.Aglobalthe sofull problemhave just described the hydrodynamic, the electromagnetic and the thermal problems. In each of thoseassumed that the other fields were blem we want to solve is to find the velocity u, the pressure p, the electrical potential /, the enthalpyd the temperature h satisfying the three problems above; the functions bH;ah;Cph;b0x;H0x andsh are given and the constants q;l0;Cp;ha;V1;r and l are known.numerical approachDh C0 h dxdy 0: 37be smal3.2. Computschemethe boundary conditions (16)(18) and then j by using (19). Subsequently, we apply BiotSavart law to findo/AstowardMwherorder0 0;(Remark.but wediscrete problem is well posed, we note that the convergence of the enthalpy approximation toward thesolutionWthermalfieldsGradientThethe conelectricshapebilizationpicturethis fieldItpart ofiliaryThehydrodynami1488 Y. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492e use GMRES to solve the matrix system resulting from the hydrodynamic problem SmHD and from theproblem SmTh. In the other hand, since the matrix system related to the computation of electricalSmEM is symmetric positive definite, we use the Algebraic Multi Grid method AMG or the Conjugatedmethod CG to solve this problem.MHD-Thermal calculation is carried out by using (PC Pentium(R) 4, CPU 2.80 GHz, 2 GB RAM) andvergence of the global algorithm is obtained in 10 hours. The results relative to the computation of thepotential are presented in Fig. 3. Fig. 4 shows the distribution of the temperature through the cell. Theof the ledge is given in Fig. 5. We observe clearly the eect of numerical diusion related to SUPG sta-on the narrow region of the solidifications front. It may be worth pointing out the agreement of thiswith the one representing the velocity field (Fig. 6) and to notice in particular that the locations whereis small correspond to the locations where the ledge is large.is easy to observe that the numerical results of the velocity field computations match Darcy law on thedomain where the liquid fraction is less than 1. This shows the eciency of using the method of aux-domain” involving a penalization of the velocity field in the hydrodynamic model.fluid layer along their depth is discretized with several elements to have more accurate representation of4. Numerical resultsH is only true in L2K norm.presents a strong gradient in a narrow regio