sdarticle (4).pdf

Numerical simulation of thermal problems coupled with magnetohydrodynamic eff ects in aluminium cell q Y. Safa*, M. Flueck, J. Rappaz Institute of Analysis and Scientifi c Computing, Ecole Polytechnique Fe de rale de Lausanne, Station 8, 1015 Lausanne, Switzerland Received 27 December 2006; received in revised form 4 February 2008; accepted 8 February 2008 Available online 29 February 2008 Abstract A system of partial diff erential equations describing the thermal behavior of aluminium cell coupled with magnetohy- drodynamic eff ects is numerically solved. The thermal model is considered as a two-phases Stefan problem which consists of a non-linear convectiondiff usion heat equation with Joule eff ect as a source. The magnetohydrodynamic fi elds are gov- erned by NavierStokes and by static Maxwell equations. A pseudo-evolutionary scheme (Chernoff ) is used to obtain the stationary solution giving the temperature and the frozen layer profi le for the simulation of the ledges in the cell. A numer- ical approximation using a fi nite element method is formulated to obtain the fl uid velocity, electrical potential, magnetic induction and temperature. An iterative algorithm and 3-D numerical results are presented. ? 2008 Elsevier Inc. All rights reserved. Keywords: Aluminium electrolysis; Chernoff scheme; Heat equation; Magnetohydrodynamics; Ledge; Solidifi cation 1. Introduction A phase changing problem motivated by the modelling of thermal problem coupled with magnetohydro- dynamic eff ects in a reduction cell is studied. In a smelting cell operating with HallHe roult process, the metal part 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 spreads in the electrolytic bath. The important magnetic fi eld generated by the currents carried to the alignment of cells, coupled with the currents running through the cells themselves gives rise to a fi eld of Laplace forces which maintains a motion within these two conducting liquids. A magnetohydrodynamic interaction takes place in the cell. In the other hand a heating source is produced by the Joule eff ect due to the electric resistivity of the bath. 0307-904X/$ - see front matter ? 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2008.02.011 q Sponsors: 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 On the wall of the cell, a solidifi ed bath layer, the so-called ledge is created. These ledges protect the cell sidewall from corrosive electrolytic bath and reduce the heat loss from the cell (see 2 page 23). Moreover, its profi le strongly infl uences the magnetohydrodynamic stability causing oscillations of the aluminiumbath interface which could decrease the current effi ciency. Consequently an optimal ledge profi le is one of the objec- tives of cell sidewall design. The thermal solidifi cation problem in smelting cell has been treated by several authors 35. As far as we are aware, this problem has never been considered when coupled with the magnetohydrodynamic fi elds. The aim of this paper is to deal with such fi elds interaction. Let us mention that the details on this problem can be found in Safas thesis 6. Mathematically, the problem is to solve a coupled system of partial diff erential equations consisting of the heat equation with Joule eff ect as a source, Maxwell law equations with electrical conductivity as a function of temperature and NavierStokes equations. The interface between aluminium and bath is an unknown. The ledge 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 presented in Section 3 and we give the numerical results in Section 4. 2. The model In order to introduce the model we fi rst describe some geometrical and physical quantities. 2.1. General descriptions The geometry is schematically defi ned by Fig. 1. We introduce the following notations: ? X X1 X2 : fl uids and solid ledge, ? N N1 N2: electrodes, ? K X N: domain representing the cell and we defi ne the interfaces: ? C oX1 oX2: free interface between aluminium and bath, which is an unknown, ? Ri oK oNi;i 1;2, ? R R1 R2: outer boundary of the electrodes. Aluminium Electrolyte Cathode Lining Anode Blocks Frozen ledgeFrozen ledge Fig. 1. Transverse cross section of aluminium reduction cell. 1480Y. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492 The unknown physical fi elds with which we shall deal are listed as follows: Hydrodynamic fi elds: ? u: velocity fi eld in Xi;i 1;2; (u 0 in solid ledges), ? p: pressure. Electromagnetic fi elds: ? b: magnetic induction fi eld, ? e: electric fi eld, ? j: electric current density. Thermal fi elds: ? H: enthalpy, ? h: temperature. The material properties are defi ned as ? q: mass density, ? rband r: electrical conductivity in and, respectively, outside the bath, ? g: viscosity of the fl uids, ? l0: magnetic permeability of the void, ? k: thermal conductivity, ? Cp : specifi c heat, ? : latent heat. 2.2. Physical assumptions The model leans on the following basic hypotheses: 1. The fl uids are immiscible, incompressible and Newtonian. 2. In each domain Xi , i = 1, 2, the fl uids are governed by the stationary NavierStokes equations. 3. The electromagnetic fi elds satisfy the stationary Maxwells equations, Ohms law is moreover supposed to be 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 fl uids and electrodes parts. 6. The viscosity g, the density q and the specifi c 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 eff ect due to the current crossing the cell. 9. Eff ects of chemical reactions 7, Marangoni eff ect 8,9, surface tension as well as the presence of gas fl ow are neglected. 2.3. The hydrodynamic problem In this part we consider the temperature fi eld h and the electromagnetic fi elds j and b as known. We choose to represent the unknown interface between aluminium and bath by a parametrization of the form C?h x;y;z : z ?hx;y;x;y 2 D?, where D is usually a rectangle corresponding to the parametrization of aluminiumcathode interface. We denote the dependence of X1;X2and C with respect to?h by using Xi Xi?h;i 1;2;C C?h: Y. Safa et al./Applied Mathematical Modelling 33 (2009) 147914921481 From assumption (vii) we get the following relation: Z D ? hx;ydxdy V1;whereV1is the volume of aluminium: The unit normal to C?h pointing into X2?h is given by n 1 krz ?hk rz ?h: We consider the following standard set of equations for hydrodynamic fi elds: qu;ru ? div2lDu ? p qgzI j bin X1?h X2?h;1 divu 0in X1?h X2?h;2 u;rz ?h 0on C?h;3 with Du 1 2 ru ruT;I diji;j 1;2;3: Here (.,.) is the usual scalar product on R3. Eqs. (1)(3) correspond to 1st and 2nd assumptions. We complete those equations by introducing the conditions on the boundaries of the domains X1?h and X2?h containing the fl uids. For any fi eld w, w?C? h denotes the jump of w across C?h, i.e. w?C? h wbath? waluminium. For the fi elds u and p we have u 0on oX;4 u?C? h 0;5 ?pI 2lDun?C? h 0:6 The fl uid part of Xi?h i = 1, 2 is only a subdomain of the domain Xi? h delimited by the front of solidifi cation. In order to solve the hydrodynamic problem in a fi xed domain Xi , we use the method of fi ctitious domain” involving a penalization tool. The velocity and the pressure will then be defi ned in both liquids and solids. We add to NavierStokes equation the term Kfsu;fsis the solid fraction which is a function of temperature. The function K is given by Carman Kozeny” law: Kfs lCf 2 s P21 ? fs3 ; where P is the mean pore size and C is a constant obtained experimentally (see 10). Eq. (1) may then be mod- ifi ed to qu;ru ? div2lDu ? p qgzI Ku j bin X1?h X2?h:7 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 mimics the Darcy law: rp qgz ?Ku j b: When fs! 1, we get Kfs ! 1 and then u 0 in the solid zone. 1482Y. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492 We fi nally obtain the hydrodynamic problem PHD : for given j;b and h, fi nd u, p and?h such that qu;ru ? div2lDu ? p qgzI Ku j bin X1?h X2?h;8 divu 0in X1?h X2?h;9 u;n 0on C?h;10 u 0on oX;11 u?C? h 0;12 ?pI 2lDun?C? h 0;13 Z D ? hx;ydxdy V1:14 2.4. The electromagnetic problem We consider the velocity fi eld u as well as the temperature h are known. From the Faradays law we have rot e 0, the electric fi eld is then given by e ?r/, where / is the electric potential fi eld computed in K. We still denote by u the continuous extension of the velocity by zero in K, taking into account Amperes law: rot b l0j and Ohms law: j r?r/ u b in K, we then obtain the electric conservation law given by div?rr/ u b 0in K: We denote by o on the operator n;r, here n is the outer unit normal on oK. We introduce the following boundary conditions concerning electric potential /: ? r o/ on 0on oK n R; ? r o/ on j0on R2; / 0on R1; where j0is the given current density on the outer boundary of the anode R2. Notice that magnetic induction b is obtained as a function of electrical current j by using BiotSavart relation: bx l0 4p Z K jy x ? y kx ? yk3 dy b0x8x 2 K; where b0 is some magnetic induction fi eld due to the electric currents which fl ows outside the cell. The electromagnetic problem PEMis then formulated as following: for given u and? h, fi nd /;b and j such that divr?r/ u b 0in K;15 ? r o/ on 0on oK n R;16 ? r o/ on j0on R2;17 / 0on R1;18 j r?r/ u bin K;19 bx l0 4p Z K jy x ? y kx ? yk3 dy b0x8x 2 K:20 2.5. The thermal problem We consider as known the hydrodynamic fi eld u and the electromagnetic fi eld j. The steady solution we are looking for will be here obtained as the limiting case of a time dependent heat equation. Y. Safa et al./Applied Mathematical Modelling 33 (2009) 147914921483 In this subsection, we thus introduce the evolutionary thermal model. In our convectiondiff usion problem the location of the front of solidifi cation (interface separating ledge and liquid bath) is not known a priori and so needs to be determined as part of the solution. Such problems, widely referred as Stefan problems”, are highly non-linear. In order to overcome diffi culty related to the non-linearity of the Stefan interface condition, an enthalpy function is defi ned, it represents the total heat content per unit volume of material. The enthalpy can be expressed in terms of the temperature, the latent heat and of the solid fraction fs, namely: Hh Z h 0 qCpsds 1 ? fsh:21 Since the enthalpy Hh is a monotonic function we can introduce the function b defi ned by the relation bHh h:22 The function bH is computed by mathematical processing (interpolation) in the list of h;H values cor- responding to the inverse relation H b?1h given in Eq. (21). With this relation we can formulate the prob- lem as a Stefan problem in temperature and enthalpy under the form oH ot ? divkhrh qCpu;rh S;23 h bH;24 which is a non-linear convectiondiff usion system. The term u;rh denotes the scalar product of u with rh, S is the heat source provided by Joule eff ect only. It takes the form S rkr/k2:25 The advantage of this temperatureenthalpy formulation, taken in distributional sense, is that the necessity to carefully track the location of solidliquid interface is removed and standard numerical technique can be employed to solve our phase change problem. The temperature h is subject to the Robin boundary condition: k oh on aha? hon oK;26 where oh on is the derivative in the direction of the outward unit normal on oK;a is the coeffi cient of thermal transfer, which may depends on both space and temperature, and hais the temperature outside K. The heat transfer is due to convection and radiation. The radiation is implicitly taken into account by using: a ah c1 c2h ? c3W=m2 ?C where c1;c2and c3are positive values provided by experimental estimation. An initial condition on enthalpy Hx;0 H0on K is assumed. For a given scalar value T, which will represent the integration time, we denote: QT K?0;TandRT oK?0;T: The thermal problem PThtakes the form: for given u;? h and j, fi nd h and H such that oH ot ? divkhrh qCpu;rh Sin QT;27 h bHin QT;28 kh oh on ahha? hon RT;29 H H0in K;for t 0:30 1484Y. Safa et al./Applied Mathematical Modelling 33 (2009) 14791492 2.6. The full problem We have just described the hydrodynamic, the electromagnetic and the thermal problems. In each of those we have assumed that the other fi elds were given. The problem we want to solve is to fi nd the velocity u, the pressure p, the electrical potential /, the enthalpy H and the temperature h satisfying the three problems above; the functions bH;ah;Cph;b0x;H0x and fs fsh are given and the constants q;l0;Cp;ha;V1;r and l are known. 3. The numerical approach The numerical solution of the mathematical above problems is based on an iterative procedure in which we carry out alternatively the computation of the three types of unknowns: hydrodynamic HD, electromagnetic EM and thermic Th. In this section we present the iterative schemes for the problems PHD;PEMand PTh. A global pseudo-evolutive” algorithm involving a space discretization by fi nite element method is applied for the solving of the three coupled problems. 3.1. Computation of the hydrodynamic fi elds The hydrodynamic problem is iteratively solved. In each solving step, we fi rst solve the problem in a fi xed geometry without normal force equilibrium condition on the interface and then we update the interface posi- tion by using the non-equilibrium normal force. The solving deals with the alternative application of the two following steps: ? Step 1: we solve the hydrodynamic problem for the given geometry C?h and by taking into account the interface conditions: u;n 0;on C?h; ?pI 2lDun;t?C? h 0;8 t tangential vector on C?h; the problem is then easily formulated on a weak formulation. ? Step 2: we update the position of interface in order to verify the equilibrium of normal forces on the inter- face C, we choose?h :?h d?h with: d?h ? ?pI 2lDun;n?C? h Cste j b;ez ? qg?C? h : Here we denote by ezthe unit vector of Oz axis and by Cste the constant obtained from the condition: Z Z D d?hx;ydxdy 0: An iterative scheme is used to compute um1;pm1;Cstem1;wm1and ? hm1as functions of the values obtained from the previous iteration m. We set w ?pI 2lDun;n?C? h and fz j b;ez ? qg and we defi ne the solving step Sm HD by qum;rum1? div2lDum1? pm1 qgzI Kum1 jm bmin X1?hm X2?hm;31 divum1 0in X1?hm X2?hm;32 um1;n 0on C?hm;33 ?pm1I 2lDum1n;t?C? hm 08 t tangent on C?hm;34 wm1 ?pm1I 2lDum1n;n?C? hm; 35 ? hm1?hm? wm1 Cstem1 fm z ?C? hm ;36 D?hm1? ? hmdxdy 0:37 Y. Safa et al./Applied Mathematical Modelling 33 (2009) 147914921485 Noting that the stop condition for this algorithm is based on H1norm estimation of um1? um, which has to be smaller than a tolerance ?. 3.2. Computation of the electromagnetic fi elds The magnetic induction b depends on electrical current j and implicitly on potential fi eld /, we have to com- pute the values of these electromagnetic fi elds for a known velocity fi eld um1 . To fi nd / we apply an iterative scheme in which, at the solving step m, we use the value of bmto compute successively /m1by using (15) and the boundary conditions (16)(18) and then jm1 by using (19). Subsequently, we apply BiotSavart law to fi nd the value of bm1as function of jm1. We carry out the solving step Sm EM by div?rr/m1 rum1 bm 0in K;38 ? r o/m1 on 0on oK n R;39 ? r o/m1 on j0on R2;40 /m1 0on R1;41 jm1 r?r/m1 um1 bmin K;42 bm1x l0 4p Z R3 jm1y x ? y kx ? yk3 dy b0x8 x 2 K:43 The stop condition is based on L2norm estimation of /m1? /mwhich has to be smaller than a tolerance ?. 3.3. Computation of the thermal fi elds As already mentioned, we use a pseudo-evolutive description as a mathematical mean which converges toward the steady solution of thermal problem (27)(30). Making use of semi-implicit discretization in time of (27)(30) we obtain Hm1? Hm s ? r ? khmrhm1 qCpum;rhm1 Sm;44 where Hm;hmand Smare the values of H, h and S at time tm mtau and s is the time step discretization. In order to close the system (44), we make use of the Chernoff scheme, namely Hm1 Hm chm1? bHm;45 where c is a positive relaxation parameter. By replacing (45) in (44), we obtain a scheme in order to compute the temperature at time tn1, i.e. c hm1? bHm s ? r ? khmrhm1 qCpum1;rhm1 Sm:46 It is shown in 11 that the