機(jī)械設(shè)計(jì)外文翻譯-耦合系統(tǒng)熱問題的數(shù)值模擬在鋁電池中具有磁流體動(dòng)力效應(yīng)【中文6270字】【PDF+中文WORD】,中文6270字,PDF+中文WORD,機(jī)械設(shè)計(jì),外文,翻譯,耦合,系統(tǒng),問題,數(shù)值,模擬,電池,具有,流體,動(dòng)力,效應(yīng),中文,6270,PDF,WORD Numerical simulation of thermal problems coupledwith magnetohydrodynamic effects in aluminium cellqY.Safa*,M.Flueck,J.RappazInstitute of Analysis and Scientific Computing,Ecole Polytechnique Fe de rale de Lausanne,Station 8,1015 Lausanne,SwitzerlandReceived 27 December 2006;received in revised form 4 February 2008;accepted 8 February 2008Available online 29 February 2008AbstractA system of partial differential equations describing the thermal behavior of aluminium cell coupled with magnetohy-drodynamic effects is numerically solved.The thermal model is considered as a two-phases Stefan problem which consistsof a non-linear convectiondiffusion heat equation with Joule effect as a source.The magnetohydrodynamic fields are gov-erned by NavierStokes and by static Maxwell equations.A pseudo-evolutionary scheme(Chernoff)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.?2008 Elsevier Inc.All rights reserved.Keywords:Aluminium electrolysis;Chernoffscheme;Heat equation;Magnetohydrodynamics;Ledge;Solidification1.IntroductionA phase changing problem motivated by the modelling of thermal problem coupled with magnetohydro-dynamic effects in a reduction cell is studied.In a smelting cell operating with HallHe roult 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 effect due to the electric resistivityof the bath.0307-904X/$-see front matter?2008 Elsevier Inc.All rights reserved.doi:10.1016/j.apm.2008.02.011qSponsors: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) 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 efficiency.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 differential equations consisting of theheat equation with Joule effect 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 presentedin 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:?X X1 X2:fluids and solid ledge,?N N1 N2:electrodes,?K X N:domain representing the celland we define 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.AluminiumElectrolyteCathode LiningAnode BlocksFrozen ledgeFrozen ledgeFig.1.Transverse cross section of aluminium reduction cell.1480Y.Safa et al./Applied Mathematical Modelling 33(2009)14791492The unknown physical fields with which we shall deal are listed as follows:Hydrodynamic fields:?u:velocity field in Xi;i 1;2;(u 0 in solid ledges),?p:pressure.Electromagnetic fields:?b:magnetic induction field,?e:electric field,?j:electric current density.Thermal fields:?H:enthalpy,?h:temperature.The material properties are defined as?q:mass density,?rband r:electrical conductivity in and,respectively,outside the bath,?g:viscosity of the fluids,?l0:magnetic permeability of the void,?k:thermal conductivity,?Cp:specific heat,?: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 effect due to the current crossing the cell.9.Effects of chemical reactions 7,Marangoni effect 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 formC?h x;y;z:z?hx;y;x;y 2 D?,where D is usually a rectangle corresponding to the parametrizationof aluminiumcathode interface.We denote the dependence of X1;X2and C with respect to?h by usingXi Xi?h;i 1;2;C C?h:Y.Safa et al./Applied Mathematical Modelling 33(2009)147914921481From assumption(vii)we get the following relation:ZD?hx;ydxdy V1;whereV1is the volume of aluminium:The unit normal to C?h pointing into X2?h is given byn 1krz?hkrz?h:We consider the following standard set of equations for hydrodynamic fields:qu;ru?div2lDu?p qgzI j bin X1?h X2?h;1divu 0in X1?h X2?h;2u;rz?h 0on C?h;3withDu 12ru 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 completethose equations by introducing the conditions on the boundaries of the domains X1?h and X2?h containingthe fluids.For any field w,w?C?hdenotes the jump of w across C?h,i.e.w?C?h wbath?waluminium.For thefields u and p we haveu 0on oX;4u?C?h 0;5?pI 2lDun?C?h 0:6The fluid part of Xi?h i=1,2 is only a subdomain of the domain Xi?h delimited by the front of solidification.In order to solve the hydrodynamic problem in a fixed domain Xi,we use the method of fictitious domain”involving a penalization tool.The velocity and the pressure will then be defined in both liquids and solids.Weadd to NavierStokes equation the term Kfsu;fsis the solid fraction which is a function of temperature.Thefunction K is given by Carman Kozeny”law:Kfs lCf2sP21?fs3;where P is the mean pore size and C is a constant obtained experimentally(see 10).Eq.(1)may then be mod-ified toqu;ru?div2lDu?p qgzI Ku j bin X1?h X2?h:7If 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 mimicsthe 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)14791492We finally obtain the hydrodynamic problem PHD:for given j;b and h,find u,p and?h such thatqu;ru?div2lDu?p qgzI Ku j bin X1?h X2?h;8divu 0in X1?h X2?h;9u;n 0on C?h;10u 0on oX;11u?C?h 0;12?pI 2lDun?C?h 0;13ZD?hx;ydxdy V1:142.4.The electromagnetic problemWe 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?r/,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:rot b l0j and Ohms law:j r?r/u b in K,we then obtain the electric conservation law given bydiv?rr/u b 0in K:We denote byoonthe operator n;r,here n is the outer unit normal on oK.We introduce the following boundary conditions concerning electric potential/:?ro/on 0on oK n R;?ro/on j0on R2;/0on R1;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:bx l04pZKjy x?ykx?yk3dy b0 x8x 2 K;where 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 and?h,find/;b and j suchthatdivr?r/u b 0in K;15?ro/on 0on oK n R;16?ro/on j0on R2;17/0on R1;18j r?r/u bin K;19bx l04pZKjy x?ykx?yk3dy b0 x8x 2 K:202.5.The thermal problemWe consider as known the hydrodynamic field u and the electromagnetic field j.The steady solution we arelooking for will be here obtained as the limiting case of a time dependent heat equation.Y.Safa et al./Applied Mathematical Modelling 33(2009)147914921483In this subsection,we thus introduce the evolutionary thermal model.In our convectiondiffusion problemthe location of the front of solidification(interface separating ledge and liquid bath)is not known a priori andso needs to be determined as part of the solution.Such problems,widely referred as Stefan problems”,arehighly non-linear.In order to overcome difficulty related to the non-linearity of the Stefan interface condition,an enthalpy function is defined,it represents the total heat content per unit volume of material.The enthalpycan be expressed in terms of the temperature,the latent heat and of the solid fraction fs,namely:Hh Zh0qCpsds 1?fsh:21Since the enthalpy Hh is a monotonic function we can introduce the function b defined by the relationbHh h:22The 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 formoHot?divkhrh qCpu;rh S;23h bH;24which is a non-linear convectiondiffusion system.The term u;rh denotes the scalar product of u withrh,S is the heat source provided by Joule effect only.It takes the formS rkr/k2:25The advantage of this temperatureenthalpy formulation,taken in distributional sense,is that the necessityto carefully track the location of solidliquid interface is removed and standard numerical technique can beemployed to solve our phase change problem.The temperature h is subject to the Robin boundary condition:kohon aha?hon oK;26whereohonis the derivative in the direction of the outward unit normal on oK;a is the coefficient of thermaltransfer,which may depends on both space and temperature,and hais the temperature outside K.The heattransfer is due to convection and radiation.The radiation is implicitly taken into account by using:a ah c1 c2h?c3W=m2?Cwhere 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,find h and H such thatoHot?divkhrh qCpu;rh Sin QT;27h bHin QT;28khohon ahha?hon RT;29H H0in K;for t 0:301484Y.Safa et al./Applied Mathematical Modelling 33(2009)147914922.6.The full problemWe have just described the hydrodynamic,the electromagnetic and the thermal problems.In each of thosewe have assumed that the other fields were given.The problem we want to solve is to find the velocity u,the pressure p,the electrical potential/,the enthalpyH and the temperature h satisfying the three problems above;the functions bH;ah;Cph;b0 x;H0 x andfs fsh are given and the constants q;l0;Cp;ha;V1;r and l are known.3.The numerical approachThe 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.Aglobal pseudo-evolutive”algorithm involving a space discretization by finite element method is applied forthe solving of the three coupled problems.3.1.Computation of the hydrodynamic fieldsThe hydrodynamic problem is iteratively solved.In each solving step,we first solve the problem in a fixedgeometry 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 twofollowing steps:?Step 1:we solve the hydrodynamic problem for the given geometry C?h and by taking into account theinterface 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 Cstej 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 ZDd?hx;ydxdy 0:An iterative scheme is used to compute um1;pm1;Cstem1;wm1and?hm1as functions of the valuesobtained from the previous iteration m.We set w?pI 2lDun;n?C?hand fz j b;ez?qg and we define the solving step SmHDbyqum;rum1?div2lDum1?pm1 qgzI Kum1 jm bmin X1?hm X2?hm;31divum1 0in X1?hm X2?hm;32um1;n 0on C?hm;33?pm1I 2lDum1n;t?C?hm 08 t tangent on C?hm;34wm1?pm1I 2lDum1n;n?C?hm;35?hm1?hm?wm1 Cstem1fmz?C?hm;36D?hm1?hmdxdy 0:37Y.Safa et al./Applied Mathematical Modelling 33(2009)147914921485Noting that the stop condition for this algorithm is based on H1norm estimation of um1?um,which has tobe smaller than a tolerance?.3.2.Computation of the electromagnetic fieldsThe magnetic induction b depends on electrical current j and implicitly on potential field/,we have to com-pute the values of these electromagnetic fields for a known velocity field um1.To find/we apply an iterativescheme in which,at the solving step m,we use the value of bmto compute successively/m1by using(15)andthe boundary conditions(16)(18)and then jm1by using(19).Subsequently,we apply BiotSavart law to findthe value of bm1as function of jm1.We carry out the solving step SmEMbydiv?rr/m1 rum1 bm 0in K;38?ro/m1on 0on oK n R;39?ro/m1on j0on R2;40/m1 0on R1;41jm1 r?r/m1 um1 bmin K;42bm1x l04pZR3jm1y x?ykx?yk3dy b0 x8 x 2 K:43The stop condition is based on L2norm estimation of/m1?/mwhich has to be smaller than a tolerance?.3.3.Computation of the thermal fieldsAs already mentioned,we use a pseudo-evolutive description as a mathematical mean which convergestoward the steady solution of thermal problem(27)(30).Making use of semi-implicit discretization in time of(27)(30)we obtainHm1?Hms?r?khmrhm1 qCpum;rhm1 Sm;44where Hm;hmand Smare the values of H,h and S at time tm mtau and s is the time step discretization.Inorder to close the system(44),we make use of the Chernoffscheme,namelyHm1 Hm chm1?bHm;45where c is a positive relaxation parameter.By replacing(45)in(44),we obtain a scheme in order to computethe temperature at time tn1,i.e.chm1?bHms?r?khmrhm1 qCpum1;rhm1 Sm:46It is shown in 11 that the scheme is stable as long as c satisfies the following condition:0 0;0if kum1k 0(54with dKisdKPeKif PeK 1;1if PeKP 1;?55and hKis the size of K,the term PeKis the local Peclet number.We denote by km khmh and am ahmh,respectively,the thermal conductivity and the heat transfer coefficient at time step m.Remark.It is obvious that the enthalpy H has a jump on the solidification front which is not a priori known,but we approach H by Hm1hwhich is in the subspaces of continuous functions.Thus the approximationpresents a strong gradient in a narrow region where the exact enthalpy makes a jump.In spit the fact that ourdiscrete problem is well posed,we note that the convergence of the enthalpy approximation toward thesolution H is only true in L2K norm.4.Numerical resultsWe use GMRES to solve the matrix system resulting from the hydrodynamic problem SmHD and from thethermal problem SmTh.In the other hand,since the matrix system related to the computation of electricalfields SmEM is symmetric positive definite,we use the Algebraic Multi Grid method AMG or the ConjugatedGradient method CG to solve this problem.The MHD-Thermal calculation is carried out by using(PC Pentium(R)4,CPU 2.80 GHz,2 GB RAM)andthe convergence of the global algorithm is obtained in 10 hours.The results relative to the computation of theelectric potential are presented in Fig.3.Fig.4 shows the distribution of the temperature through the cell.Theshape of the ledge is given in Fig.5.We observe clearly the effect of numerical diffusion related to SUPG sta-bilization on the narrow region of the solidifications front.It may be worth pointing out the agreement of thispicture with the one representing the velocity field(Fig.6)and to notice in particular that the locations wherethis field is small correspond to the locations where the ledge is large.It is easy to observe that the numerical results of the velocity field computations match Darcy law on thepart of domain where the liquid fraction is less than 1.This shows the efficiency of using the method of aux-iliary domain”involving a penalization of the velocity field in the hydrodynamic model.The fluid layer along their depth is discretized with several elements to have more accurate representation ofhydrodynamic field.Further,as mentioned before,the nodes at interface aluminiumbath are allowed to move1488Y.Safa et al./Applied Mathematical Modelling 33(2009)14791492along vertical line to insure the vertical force equilibrium(Step 2 in the Section 3.1).The main value of errorwith respect to the liquid depth is 6%.At the level of aluminiumbath interface we have the highest convectioneffects and consequently the finest thickness of the ledge see Fig.5It is important to insure that the total heat dissipation matches the internal heat generation.It is critical thatthis is achieved in order to consider the results are corresponding to cell steady-state conditions.The total heatdissipation is obtained from the sum of Joule heat produced in each part Kj,j 1;.;N,where N is the num-ber of the parts of the cell traversed by the electrical current:total heat dissipation XNj1ZKjrkr/k2dx 408:50 kW:The total heat lost corresponds to convective dissipation on the cells boundary oK:ZoKaha?hdC 398:3 kW:Fig.3.Electric potential results in the cell.Fig.4.Temperature results in the cell.Y.Safa et al./Applied Mathematical Modelling 33(2009)147914921489Fig.6.Velocity field shown before solidification(above)and after solidification(below)wtih mean value=0.8 cm/s.Fig.5.Liquid fraction showing the ledge.1490Y.Safa et al./Applied Mathematical Modelling 33(2009)14791492The numerical error relative to total Joule dissip