Simulation Num´erique des cuves pour l’Electrolyse ^´ de l’Aluminium

Michel Flueck -EPFL Lausanne, CIRM 25 octobre 2006.

ASN/IACS/SB/EPF -Lausanne 25 octobre 2006

Collaboration : EPFL -Industrie de l’aluminium

´

Stabilit´e Lin´eaire des Cellules d’ Electrolyse de l’Aluminium

Michel Fl¨uck, ASN/ SB/ EPFL Lausanne.

Collaboration scientifique:

´

• EPFL, Ecole Polytechnique F´ed´erale de Lausanne, chaire d’Analyse et Simulation Num´eriques, prof. Jacques Rappaz.
• ALCAN -P´echiney (Alusuisse).

Financement:

´

• EPFL, Ecole Polytechnique F´ed´erale de Lausanne
• ALCAN -P´echiney (Alusuisse).
• Fonds nationaux Suisses.

Electrolyse de l’aluminium A

• alumine extrait du minerai bauxite: Al₂O₃
• r´eaction d’´electrolyse: 2Al₂O₃ + ´energie → 4Al +3O₂
• consommation des anodes: C + O2 → CO2
• m´elange ´electrolytique : cryolithe (960 degr´es C)

Electrolyse de l’aluminium B -Cuve numerique´

Electrolyse de l’aluminium BB -Vue de dessus

Electrolyse de l’aluminium C -cuves end to end

Electrolyse de l’aluminium D -cuves side by side

Electrolyse de l’aluminium E -en Chiffres

Donn´ees g´en´erales

• dimensions d’une cuve 10 × 3 × 0.5 [m]
• 5 Volts continu / 1-500’000 Amp`eres
• 2Al₂O₃(204) + 3C(36) → 4Al(108) + 3CO₂(132)
• 300 kg C / tonne Al, 1200 kg CO₂ / tonne Al
• 5640 kWh / tonne Al (th´eorique)
• 100 cuves par halle branch´ees en s´erie: 50MW
• 1 halle neuve � 1G Frs
• 25 M tonnes Al/an

´

Electrolyse de l’aluminium -X

Les faits importants -Mod´elisation

Les acteurs principaux sur le rendement

• Mouvement des fluides: forces de Lorentz ^�j ∧^�b -courant ^�j (6000A/m²): g´eom´etrie int´erieure de la cuve (´electrolyte) -induction ^�b: g´eom´etrie des barres de distribution du courant
• L’effet d’´ecran du container en acier
• Changements d’anodes

• Distance interpolaire (thermique) Observations -Mesures
• Variations du courant dans le ”tiges anodiques”
• Fr´equences ”gravitationnelles”
• Effets anodiques
• Court circuit Buts
• Maximiser le courant I_tot traversant la cuve et donc l’aluminium produit par une cuve
• Retrofit
• Cellule le plus stable possible (rendement) Simulations Num´eriques
• Solution stationnaire et Stabilit´e lin´eaire de cette solution
• Evolutif
• Solidification

Resultats Num´

eriques typiques: Solution Stationnaire

potentiel ´electrique stationnaire

Interface stationnaire

Vitesse stationnaire au niveau de l’interface

Resultats Num´e

eriques typiques: Diagramme de Stabilit´

valeurs propres MHD

Le mod`ele physique -A

Domaines

Schematic cross section of a cell

• Ω = domaine des fluides
• Λ = la cellule compl`ete
• Γ = interface entre aluminium et electrolyte
• S = autres fronti`eres inter-mat´eriaux
• extra bus bars

Le modele physique `-B

Hypoth`eses physiqyes

Le mod`ele physique:

• Aluminium et bain: Fluides Newtonniens, visqueux, incompressibles
• Immiscibilit´e des fluides `a leur interface
• Situation ´electromagn´etique statique

Hypoth`eses:

• Temp´erature constante dans la cellule
• Tension superficielle `a l’interface n´egligeable
• R´eactions chimiques n´eglig´ees
• Effet des bulles de gaz n´eglig´e

Lemod`-C

ele physique

Champs inconnus

Champs hydrodynamiques (limit´es `a l’aluminium et au bain):

• �u(x, y, z, t) : vitesse d´efinie dans Ω
• p(x, y, z, t) : pression d´efinie dans Ω
• z − h(x, y, t) = 0: param´etrisation de l’interface
• z = h = 0: position de l’interface au repos

Champs ´electromagn´etiques(dans tout l’espace E):

• �e : champ ´electrique
• ^�b : champ d’induction magn´etique
• ^�j : courant ´electrique

Propri´et´es physiques et param`etres (constants par mat´eriau):

• σ : conductivit´e ´electrique
• ρ : densit´e des fluides
• µ : viscosit´e des fluides
• I_tot : courant total traversant la cellule

Le mod`ele physique -D

Le probl`eme hydrodynamique

´

Equations pour les fluides: Navier-Stokes incompressible

ρ(∂_t�u +(�u.�)�u) − div τ(�u, p)= ^�j ∧^�b − ρ�g dans Ω div �u = 0 dans Ω

´

Equations pour la position de l’interface:

∂_th − �u.�(z − h)=0 surΓ(h)

Conditions aux bords:

�u = 0 sur ∂Ω

Conditions `a l’Interface: ([•]_Γ = •_bain −•_alu, �n normale unit´e) [�u]_Γ(h) = 0 sur Γ(h) [τ_i,j(�u, p)n_j]_Γ(h) =0 i =1, 2, 3, sur Γ(h)

Mod`ele de turbulence anisotrope τ_i,j(�u, p)= −pδ_i,j + µ_i,j(∂_iu_j + ∂_ju_i) , i, j =1, 2, 3

Le mod`ele physique -E

Le probl`eme ´electromagn´etique

´

Equations de Maxwell Statique: rot ^�b = µ₀^�j dans E rot �e = 0 dans E div ^�b = 0 dans E

´

Courant Electrique :

e = −^�

��ϕ dans Λ

^�j = σ(−^�u ∧^�dans Λ

�ϕ + �b)

j fix´e dans barres div ^�j = 0 dans Λ

Conditions aux bords et `a l’interface:

b, ϕ, j.�n continu sur Γ(h) ∪ S

− j.�n dσ = j.�n dσ = I_tot

Γin Γout

j.�n =0 Γ_isole

o`u I_tot = courant total traversant une cuve (constant, fix´e)

Question: solution stationnaire ? -stabilit´e ?

Stabilit´eaire

e Lin´

• probl`eme ´evolutif Y^�(t)= F(Y(t)), Y(0) = Y₀ avec condition initiale Y₀
• → ´evolution de cette condition initiale Y₀
• solution stationnaire 0 = F(y), Y(0) = y
• → y = condition initiale qui reste fixe
• Question: Si on perturbe cette condition initiale stationnaire?

• ”petites” perturbations autour de la solution stationnaire: Y(t)= y + Y (t)

iωt

• d´ependence en temps: eY (t) := Ye^iωt = Ye^−Im(ω)t {cos Re(ω)t + i sin Re(ω)t}
• lin´earisation du probl`eme ´evolutif: iωY = �F(y).Y ,

→ probl`eme aux valeurs propres pour (ω, Y )

• mode(ω, Y ) Im(ω) ≤ 0 la perturbation Y croˆıt, mode INstable Im(ω) > 0 la perturbation Y d´ecroˆıt, mode stable.

Viscosite´

Hypoth`eses Suppl´ementaires

1^o) Solution stationnaire: Navier-Stokes

• �u = 0 sur ∂Ω,
• [�u]_Γ =0 sur Γ, �u.�n =0 sur Γ

2^o) Perturbation: on suppose µ =0 equations d’Euler au lieu de Navier-Stokes

U.�n = 0 sur ∂Ω

�		�
•	�U.�n		Γ	= 0 sur Γ,

• p _Γ =0, [p]_Γ =0
• prise en compte a posteriori de la viscosit´e: terme de Moreau-Evans λ�u en place de µΔ�u

La solution stationnaire

L’Algorithme It´eratif

It´erations pour obtenir la surface libre

• G´eom´etrie (interface) fix´ee
• Calcul des forces de Lorentz
• Potentiel, Courant
• Induction
• Calcul des vitesses
• avec immiscibilit´e des fluides
• sans l’´equilibre des forces `
• Correction de l’interface

a l’interface

La solution stationnaire

L’Algorithme It´eratif -D´etails 1

00

Position initiale �u =0,h ⁰ =0, ^�b = 0, Poser k = 0:

• calcul du potentiel ϕ ^k+1 ∈ V_ϕ tq pour tout w ∈ V_ϕ:

k

_�ϕ k+1 k

_σ^k ���

.�w = σ^k(�u ∧^�b ).�w + gwds

ΛΛ Γin∪Γout k+1 _�ϕk+1

• calcul du courent ^�j = σ^k(−^�+ �u^k ∧^�b^k)
• calcul de l’induction (Biot-Savart)

� k+1

k+1 j (�y) ∧ (�x − �y) 0

^�b (�x)= d�y +^�b (�x)

| �x − �y |³

Λ k+1 k+1 k+1

• calcul des forces f^�= ^�j ∧^�b − ρ�g

La solution stationnaire
L’Algorithme It´eratif -D´etails 2
• Navier-Stokes �u u, k+1 ∈ V�	p k+1 ∈ Vp,		ψ k+1 ∈ Vψ tq:
• 2 µij�ij(�u )�ij(�v) − p k+1div �v + � � Ω Ω k k+1� � Ω Ω ρ k(���k+1 k+1k+1 + u .�)�u .�v = f .�v,				� Γ k ψ k+1 �kv.�n ds	∀�v ∈ V�u
• q div �u k+1 = 0, � Ω					∀q ∈ Vp
� Γ k • η �u .�n ds = 0, k+1k					∀η ∈ Vψ
• correction de l’ interface hk+1 = h k +		� � τi,j(�u , p i n k+1 k+1)n k k j Γ k [fz]

Stabilit´eaire A

e Lin´

Formulation en complexes standard

Solution stationnaire:

• h(x, y),
• �u(x, y, z),...

Syst`eme d’´equation lin´earis´ees pour

• H(x, y, t)

• U^�(x, y, z, t),... avec conditions initiales

• h(x, y)+ H(x, y, 0)
• �u(x, y, z)+ U^�(x, y, z, 0),...

D´ependance en temps

• H(x, y, t)=[e^αt(a cos βt + b sin βt)] H(x, y)
• U^�(x, y, z, t)=[e^αt(a cos βt + b sin βt)] U^�(x, y, z),...

Formulation complexe

• H(x, y, t)= Ree_� ^iωtH(x, y) _� iωt �
• U^�(x, y, z, t)= ReeU(x, y, z) ,...

Stabilit´eaire B

e Lin´

Lin´earisation -probl`eme de l’interface

Solution stationnaire:
• h(x, y), • p1(x, y, z), • p2(x, y, z),	Ω1(h), �u1(x, y, z), �u2(x, y, z),	Ω2(h),	Γ(h)) (x, y, z) ∈ Ω1(h), (x, y, z) ∈ Ω2(h),

Stabilit´eaire B

e Lin´

Lin´earisation -probl`eme de l’interface

Solution perturb´ee :
• Γ(h + H)),	Ω1(h + H),	Ω2(h + H)
• h(x, y) + H(x, y, t)
• p1(x, y, z) + P1(x, y, z, t),		(x, y, z) ∈ Ω1(h + H)
• p2(x, y, z) + P2(x, y, z, t),		(x, y, z) ∈ Ω2(h + H)

Stabilit´eaire B

e Lin´

Lin´earisation -probl`eme de l’interface

• extension de p₁ `
• aΩ<sub>1</sub>(h + H)
• extension de p<sub>2</sub> `
• aΩ₂(h + H)
• [p]_Γ = 0 sur Γ(h + H) devient: p₂(x, y, h(x, y)) + ∂_zp₂(x, y, h(x, y)) + P₂(x, y, (h + H)(x, y))

−p₁(x, y, h(x, y)) − ∂_zp₁(x, y, h(x, y)) − P₁(x, y, (h + H)(x, y)) = 0 finalement:

p + H∂_zp +[P ]_Γ(h+H) =0

Γ(h) Γ(h)

• stationnaire: p =0
• stabilit´e: H∂zp +[P ]_Γ(h+H) =0

Γ(h)

Γ(h)

Stabilit´eaire C

e Lin´

Le probl`eme aux modes propres (hydro Euler)

Probl`eme Hydro-dynamique :
iωρ�U + ��P				=	�F	dans Ω
div �U = 0						dans Ω
Equation de l’interface:
iωH − �U.�n = G						sur Γ
Conditions `a l’Interface et de bords:
�U.�n = 0						sur ∂Ω
�		�
�U.�n			= 0			sur Γ
Γ
[P ]Γ = [ρ]Γ gH						sur Γ
Termes de couplage:
�F	= �j ∧ �B + �J ∧�b					dans Ω
− ρ(�U.��) �u − ρ(�u. ��) �U
G = ∂z�u.�(z − h)H − �u.�H						sur Γ

Stabilit´eaire D

e Lin´

Le probl`eme aux modes propres (magneto)

´

Probl`eme Electromagn´etique rot B^�= µ₀J^�dans E div B^�= 0 dans E

⎧ � ��� ⎨�₋�_{U ∧}�_{u ∧} ��

J = σ �Φ+ ^�b + �B +Hj δ_Γ dans Λ

Γ

⎩�

J = 0 dans E − Λ div J^�= 0 dans E

+ ... Conditions `a l’Interface et de bords

´

Etant donn´e (^�esoudre le probl`electromagn´etique

U,H),onpeutr´eme´

pour Φ,J and B^�Importance des courans stationnaires horizontaux

Modes MHD -A

Le probl`eme pour les modes MHD

´

Toutes Contributions Electromagn´etiques contenues dans: _F�₍�_J(�_{b +} �_{j ∧ B}�₍�

U,H)= ^�U,H) ∧^�U,H) − ρ(�u.^�U − ρ(^��) �u.

�) ^�U.^�

F^�et G sont lin´eaires en (^�

U,H). Le probl`eme pour les modes MHD (limit´e a Ω`):

•	iωρ�U + ��P = �F (�U, H)	dans Ω
•	div �U = 0	dans Ω
•	iωH − �U.�n = G(�U, H)	sur Γ
•	[P ]Γ = [ρ]Γ gH	sur Γ

Modes MHD -B

Le probl`eme pour les modes MHD -R´esultats Th´eoriques

Modes Gravitationnels( F^�=0,G = 0)

• Pour I_tot = 0: ∃ ω_N < ...

0

00

0102

=0 <ω

<ω

< ... < ω

les valeurs propres sont r´eelles avec un espace propre de dimension finie. Modes MHD

• Pour de ”petites” valeurs de I_tot: _ωItot _,ωItot _,ωItot

210

, ... ω^Itot , ...

N

le spectre est ”peu” chang´e, mais peut devenir complexe. On s’int´eresse aux N premi`eres valeurs propres

Un cas particulier

La celule sans courant ´electrique

La solution stationnaire sans courant ´electrique I_tot = 0:

• Force ´electromagn´etique nulle
• b =0,j =0,ϕ =0
• Fluides soumis `a la gravit´e

�u =0,p = −ρgz + C, h =0

Stabilit´e lin´eaire de cette solution:

• Tous les ”modes gravitationnels” sont p´eriodiques (→ stable)
• Si Ω est un parall´elipip`ede rectangle:

Solution Analytique pour tous les modes propres

Les ”premiers” modes propres pour l’interface H

MHD Modes -C

M´ethodologie de R´esolution par Puissance Inverse avec Shift

1^o) Calcul de la solution stationnaire

2^o) Calcul des N premiers modes parall´elipip`ediques

3^o) Calcul des N premiers modes gravitationnels `ees parall´ediques ω, H

a partir des donn´elipip`

4^o) Continuation sur le courant total I_tot = γI_nominal, with 0 <γ< 1: calcul des modes MHD `es avec γ pr´edent

a partir des modes calcul´ec´

5^o) ”Diagramme de Stabilit´e”: trajectoires des valeurs propres dans le plan complexe ω en fonction du courant total I_tot traversant la cellule.

MHD Modes -E La formulation (�U, H) Rappel du probl`eme: • iωρ�U + ��P = �F (�U, H) • div �U = 0 • iωH − �U.�n = G(�U, H) • [P ]Γ = [ρ]Γ gH	dans Ω dans Ω sur Γ sur Γ
La formulation (�U, H) (Z, vitesses `a divergence nulle): Trouver (�U, H) ∈ Z, ω ∈ C tq ∀(�V , K) ∈ Z: ω �� Ω ρ�U.�V + � Γ [ρ]Γ gHKdσ � = �� � � −i [ρ]Γ g �U.�nK − [ρ]Γ g�V .�nH Γ Γ�� � � −i �F (�U, H). �V + [ρ]Γ gG(�U, H).K

ΩΓ

Re-formulation du probl`eme dans l’espaceZ :

. Trouver z ∈ Z, ω ∈ C tq ∀ζ ∈ Z ωr(z, ζ)= s(z, ζ) .

MHD Modes -E

La formulation (P, ^�

U,H) + Puissance Inverse avec Shift

Rappel du probl`eme:

U + ^�F (^�

• iωρ^��P = ^�U,H) dans Ω
• div U^�= 0 dans Ω
• iωH − ^�n = G(^�sur Γ

U.�U,H)

• [P ]_Γ =[ρ]_Γ gH sur Γ

Puissance Inverse avec Shift:

.
´ Etant donn´e (Pn, �Un, Hn) et ωn trouver (Pn+1, �Un+1, Hn+1) et ωn+1					= ωn + δωn			tq:
•		iωn �Un+1 + i(δωn)�Un +		1 ρ ��Pn+1		=	1 ρ �F (�Un+1, Hn+1)		dans Ω
• div �Un+1 = 0									dans Ω
•			iωnHn+1 + i(δωn)Hn − �Un+1.�n = G(�Un+1, Hn+1)						sur Γ
• �			[Pn+1]Γ = [ρ]Γ gHn+1						sur Γ
•	Γ Pn+1Hndσ = q.
.

Vector de d´epart: (P₀,U₀,H₀) and ω₀

(1)

”´eliminer” P_n+1 et δω_n comme fonctions de (U^�_n+1,H_n+1)

(2)

r´esoudre le probl`eme lin´eaire pour (U^�_n+1,H_n+1) avec GMRES

MHD Modes -F

La formulation (U, H^�) + P´enalisation de ’incompressibilit´e

. Trouver z ∈ Z^�,ω ∈ C tq ∀ζ ∈ Z^� ωr(z, ζ)= s(z, ζ)+ λp(z, ζ) .

• Z^� n’a plus de contrainte d’incompressibilit´e
• pollution du spectre: contrˆolable
• utilisation de solveur aux valeurs propres ARPACK

Approche Num´

erique

Solution stationnaire

Approche Num´

erique

Solution stationnaire

Approximation ´el´ements finis

• interpolation P₁ pour �u, p, h, ϕ, b
• interpolation P₀ for ^�j, f^�
• stabilization de type GLS pour Navier-Stokes

Solveurs

• iterative solvers :GMRES/ILU(0)
• Biot-Savart pour l’induction

Resultats Num´

eriques typiques: Solution Stationnaire

potentiel ´electrique stationnaire

Interface stationnaire

Vitesse stationnaire au niveau de l’interface

Approche Num´

erique

Stabilit´e Lin´eaire

Maillage hexa´edrique

Interpolation des champs stationnaires

Approximation ´el´ements finis

• interpolation Q₁ pour P
• interpolation Q₀/8 pour H, U^�

Solveurs

• Directs
• GMRES pour les modes MHD

Resultats Num´e

eriques typiques: Diagramme de Stabilit´

valeurs propres MHD

Mod`

ele stationnaire plus complet: Thermique A

Talus de cryolithe solidifi´ee

• trop mince: ne prot`ege pas assez la structure de la cuve
• trop ´epais: d´eroute le courant ´electrique
• isolation thermique

→ Mod`ele incluant la thermique (solidification)

Mod`

ele stationnaire plus complet: Thermique B

Mod`ele temp´erature (θ) -enthalpie (H):

∂H j²

•− div (k�θ)+ ρC_p �u.�θ = dans Λ

∂t σ

• θ = β(H) relation constitutive

∂θ

• k = α(θ_ambiant − θ) sur le bord ∂Λ

∂n

• H = H₀ condition initiale

• effets Joule dans l’´electrolyte
• sch´ema de Chernoff -P1/P1 (pseudo-temps)
• sch´ema it´eratif combinant le stationnaire existant (^�j, �u)
• terme de p´enalisation Karman-Cozeny (vitesse nulle dans le talus)
• interface ?

Vitesse modifi´

ee

Temp´

erature

Talus

Mod`e: Ferromagn´

ele stationnaire raffin´etisme A

Dessins: A. Janka

• la cuve compl`ete est construite dans un caisson d’acier (ferromagn´etique) de 2cm ´epaisseur

→ effet d’´ecran contre l’induction due aux cond. ext´erieurs

Mod`e: Ferromagn´

ele stationnaire raffin´etisme B

Magn´etostatique sans ferro(ici �j est donn´e):
rot �H0 = �j,	dans E
div �B0 = �0,	dans E
�B0 = µ0 �H0,	dans E
| �B0(�x) |= O(| �x |−2)	lorsque | �x |→ ∞
Magn´etostatique avec ferro:
rot �H = �j,	dans E
div �B = �0,	dans E
�B = µ(| �H |) �H, | �B(�x) |= O(| �x |−2)	dans E lorsque | �x |→ ∞
Formulation potentiel scalaire: �H − �H0 = −�ψ dans E.
Formulation potentiel vecteur: �B − �B0 = rot �A dans E.
Dans l’algorithme stationnaire:

• apr`es chaque calcul de B^�₀, on ajoute l’”effet” du caisson: B^�− B^�₀.

Mod`e: Ferromagn´

ele stationnaire raffin´etisme C

Formulation potentiel en scalaire

Syst`eme dans tout l’espace pour ψ:

− div µ_r�ψ = −div (µ_r − 1)H^�₀, dans Λ

[µ_r�ψ.�n]= −(µ_r − 1)H^�₀.�n, sur ∂Λ

| ψ(�x) |= O(| �x |⁻¹), lorsque | �x |→ ∞ Δψ =0, dans E − Λ

• Methode 1: on remplace le probl`eme ext´erieur par une relation int´egrale entre

dψ

ψ et sur ∂Λ → matrice PLEINE

d�n

• Methode 2: 2-boules (A. Janka) : -on plonge Λ dans 2 boules concentriques -on r´esout dans la grande boule avec des cond. de bord donn´ees -on corrige ces valeurs au moyen de la Formule de Poisson appliqu´ee entre la petite

sph`ere et la grande sph`ere

Mod`e: Ferromagn´

ele stationnaire raffin´etisme D

Induction sans -avec effet du ferromagn´etisme

Le mod`evolutif ´

ele

Motivations

• calcul de transient
• contrˆole de modes in/stables

M´ethode

• maillage se d´eforme -m´ethode ALE (P1-P1, GLS, GMRES)
• vitesse nulle au bord
• m´ethode ”level set” pour l’interface (AMG)
• calcul de l’induction (B.-S.) pas trop souvent

Validation

• converge vers la mˆeme solution stationnaire
• modes gravitationnels retrouv´es

Conclusions

Mod`ele physique

• effets MHD
• effets ferromagn´etiques
• solidification du bain → ledge

Calculs dans des configurations de cellules industrielles

• solution stationnaire
• stabilit´e lin´eaire
• ´evolutif

Perspectives

• effets des bulles de gaz
• viscosit´e mieux prise en compte pour la stabilit´e lin´eaire