Le modèle thermo-hydraulique 3D dynamique a été implémenté dans le code commercial ANSYS-FLUENT qui utilise la méthode de résolution par volumes finis. Dans un souci de compréhension, on rappelle ici les hypothèses du modèle numérique développé [Chatelain16].
Transferts de chaleur III.1.2.1.
Dans l’objectif de pouvoir simuler le procédé de solidification dirigée avec des conditions thermiques les plus proches possibles de la configuration réelle, il est nécessaire de résoudre précisément les transferts thermiques entre les différents composants situés à l’intérieur de l’enceinte du four.
Si la convection dans la phase gazeuse (atmosphère d’argon) est négligeable devant les échanges par rayonnement, la convection dans le silicium liquide a quant à elle une influence forte sur la solidification dirigée puisqu’elle modifie les gradients de température à l’interface solide/liquide du silicium. Il est donc nécessaire de décrire la convection naturelle dans le bain de silicium liquide. Les équations de Navier-Stokes sont donc résolues dans le silicium avec une modélisation de la turbulence par un modèle Reynolds Averaged Navier-Stokes (RANS), et la convection naturelle est prise en compte par l’hypothèse de Boussinesq.
L’équation de la chaleur est résolue dans l’ensemble du domaine et le rayonnement est modélisé par l’approche Discret Ordinates proposée par FLUENT [ANSYS, Inc11]. L’équation de transport du rayonné est alors résolue pour un nombre fini d’angles solides (associés à différentes directions de l’espace). L’intensité du rayonnement émis et absorbé par les surfaces solides, supposées opaques, est définie par leur émissivité .
On rappelle que le modèle numérique est délimité par l’assemblage isolant du four (voir figure 3.1). La cuve métallique qui sert d’enceinte extérieure n’est pas modélisée. Des conditions aux limites sont donc appliquées sur les faces externes des isolants afin de retranscrire les pertes de chaleur vers le milieu ambiant. Le flux de chaleur vers l’extérieur comprend alors un terme lié au flux par convection et un terme lié au flux par rayonnement. On suppose que la température ambiante et la température de la cuve sur laquelle rayonnent les isolants est de 300K. Le coefficient d’échange convectif est fixé à 25 W.m-2
.K-1. Cette valeur a été déterminée préalablement par des simulations 2D de l’ensemble du four (incluant la cuve extérieure).
Chapitre III
75 On suppose également des contacts parfaits entre les différents solides du modèle, ce qui revient à négliger les résistances thermiques de contact.
Modèle enthalpique fusion/solidification du silicium III.1.2.2.
La solidification dirigée est un problème complexe faisant intervenir les transferts de chaleur dans l’enceinte du four, les mouvements convectifs dans la phase liquide et le changement de phase à l’interface solide/liquide. Il s’agit d’un problème à frontière libre puisque la forme et la position de l’interface solide/liquide sont inconnues à priori. La simulation des procédés de solidification nécessite une procédure de suivi de l’interface afin d’y imposer les conditions aux limites sur le champ de vitesse et les flux de chaleur et de masse. Il existe deux approches pour la simulation numérique des interfaces solide/liquide.
La première, appelée méthode enthalpique, consiste à définir l’intervalle de température dans lequel le changement de phase se produit et à calculer en chaque point une variable appelée fraction liquide. L’interface est donc définie implicitement par une isovaleur de la fraction liquide. Cette technique utilise une grille de calcul fixe mais nécessite un raffinement constant du maillage afin d’obtenir une discrétisation suffisante pour la transition solide/liquide.
La seconde méthode, généralement appelée méthode de maillage adaptatif (ou maillage déformable), consiste à définir explicitement l’interface solide/liquide. Pour cela, les nœuds du maillage de l’interface doivent être contraints à suivre son déplacement et sa déformation. Le déplacement de l’interface est calculé à chaque pas de temps par un bilan de flux et le maillage se déforme pour s’adapter à cette nouvelle position. Cette technique permet d’appliquer un raffinement local à l’interface, réduisant ainsi les temps de calculs par rapport à la méthode enthalpique (à précision égale). La définition des conditions aux limites est également simplifiée par la connaissance explicite de l’interface.
C’est la première méthode qui est proposée dans FLUENT [ANSYS, Inc11]. Le modèle mis en œuvre est celui proposé par Voller and Prakash [Voller87].
Afin de décrire le changement de phase solide/liquide, la chaleur latente doit être définie en fonction de la température. On utilise alors la fraction liquide locale fL pour exprimer ΔH(T) de la manière suivante : ( ) { Equ. 3. 18
Ts et Tl désignent respectivement les températures du solidus et du liquidus et L représente la chaleur latente de changement de phase, qui vaut 1.79 × 106 J·kg−1 pour le silicium [Mills00]. La fraction liquide locale fL varie entre 0 (matériau solide) et 1 (matériau liquide). La zone de changement de phase (0 < fL < 1) est appelée zone pâteuse.
La fermeture du problème de solidification nécessite de définir la fraction liquide locale en fonction de la température. Pour une solidification avec zone pâteuse (Ts≠ Tl) on suppose une évolution linéaire de fL, qui s’exprime par la relation suivante :
Chapitre III 76 { Equ. 3. 2
Cette formulation ne peut être appliquée dans le cas d’une solidification isotherme (Tl = Ts). La fraction liquide est alors calculée à l’aide d’une méthode alternative, proposée par Voller and Swaminathan [Voller91]. Le lecteur peut se référer à la publication originale pour une présentation détaillée du calcul. On notera que cette méthode conduit implicitement à la définition arbitraire d’un intervalle de température, faible mais non nul, pour le changement de phase.
Dans la zone pâteuse, le matériau est traité comme un milieu poreux. La vitesse du fluide est alors amortie et tend vers 0 lorsque la fraction liquide locale diminue.
Pour la solidification d’un corps pur, avec un front lisse, les températures du solidus et du liquidus sont donc égales (Ts = Tl) et la transition solide/liquide s’effectue de manière brutale, ce qui est équivalent à une zone pâteuse d’épaisseur nulle. On définit donc un coefficient d’amortissement élevé afin d’imposer des vitesses nulles dès que la fraction liquide devient inférieure à 1 et de conserver une zone pâteuse "numérique" la plus fine possible, dans le meilleure des cas correspondant à l’épaisseur d’une maille (5 mm dans notre cas).
Discrétisations spatiale et temporelle III.1.2.3.
Le volume du silicium a été discrétisé par un maillage caractéristique de 5 × 5 × 5 mm3. Un raffinement du maillage est appliqué aux interfaces creuset/silicium et silicium/argon afin d’améliorer le calcul des couches limites hydrodynamiques. Les mailles adjacentes aux parois ont une hauteur de 0.5 mm et un taux de croissance de 1,3 est imposé pour les mailles suivantes. Dans les autres parties du four, les mailles ont des dimensions caractéristiques comprises entre 5 et 15 mm. Le maillage complet comprend 1,6 million d’éléments. La figure 3.2 présente une vue en coupe du maillage utilisé.
Pour la simulation du procédé de cristallisation, l’amplitude des pas de temps est fixée à 1 minute. Cette discrétisation temporelle semble raisonnable compte tenu de la vitesse de déplacement de l’interface (de l’ordre de quelques micromètres par seconde) et de la durée totale du cycle. On utilise un solveur ségrégé et des schémas numériques au second ordre pour la discrétisation spatiale et temporelle.
Chapitre III
77 Figure 3. 2. Vue en coupe du maillage du four de solidification dirigée [Chatelain16].