Résultats


Adaptation du programme au nouveau maillage


Pour les raisons évoquées dans la rubrique précédente nous avons récupéré le maillage de l'étude menée par M.Konan sur le huitième de la géométrie simplifiée 3D du cubilot. Cela nous a permis de lancer nos premiers calculs sous Neptune afin de vérifier la transportabilité du programme sur les ordinateurs de l'ENSEEIHT.

Cependant, nous nous sommes rapidemment heurtés à un problème de compatitilité entre les processeurs du CINES utilisés par M.Konan, et ceux de l'ENSEEIHT. Cet obstacle a beaucoup retardé l'avancement du projet puisque cela arrétait systématiquement les calculs juste après l'initialisation avant même la première itération, ce qui empêchait toute exploitation des simulations et rendait plus difficile la compréhension des modules que nous ne pouvions réellement tester. l'erreur provenait de l'interprétation du "zéro" qui est différente entre les ordinateurs du CINES et de l'ENSEEIHT, ce dont nous nous sommes rendus compte après discussion avec M.Konan et de multiples tentatives-test notamment sur les valeurs des fractions massiques initiales. Nous obtenions un message d'arrêt immédiat de Neptune dû à la mention "NaN" au lieu des valeurs calculées des paramètres dans le fichier "listing" donnant les résultats des convergences du calcul de Neptune, cela pouvant correspondre à des valeurs qui sont infinies ou divisées par zéro. Une modification par M.Konan du code fortran "ustssp" calculant les termes sources des scalaires passifs que nous suivons, a dans un premier temps permis de visualiser quelques itérations mais sans pour autant rendre le calcul robuste pour des jeux de paramètres très différents ou un nombre d'itération élevé. Ainsi, dans un second temps nous avons réalisé une importante opération de "débugage" des modules en regard des interprétations qu'en faisaient nos processeurs. Cela a aboutit à la mise en place de plusieurs boucles de détrompages et initialisation systématique de paramètres dont nous n'étions pas sûr. Ces ajustements sont disponibles dans les annexes.

Nous avons réalisé deux maillages de même géométrie et de même dimension, mais avec un nombre de mailles différent. Cela nous a permis de lancer des calculs sur le maillage grossier afin de vérifier rapidement la stabilité de la résolution, puis de relancer le calcul sur plus de temps afin d'obtenir une résolution plus précise avec le maillage fin. Ces simulations nous ont confortés dans l'idée que nos modifications avaient été concluantes, les calculs étant alors stables quelque soit les paramètres numériques ou thermodynamiques choisis. Les temps de calculs étant très importants pour tenter de simuler les 8 heures physiques de combustion, au vu du délai imparti pour ce BEI, les résultats proviennent uniquement de l'exploitation du maillage le plus lâche. Des simulations sur plusieurs milliers d'itérations ont été lancées avec succès sur le maillage fin mais cela nécessitera beaucoup plus de temps lors d'une éventuelle reprise ultérieure de ce projet, mais la stabilité est assuré par nos tests.


Inhomogénéité spatiale de la réduction des particules de coke


Les résultats que nous obtenons permettent d'observer une inhomogénéité spatiale de la réduction des particules de coke. En effet, lors de la visualisation des résultats sous Paraview, nous pouvons avoir accès à l'évolution temporelle de la réduction des diamètres des particules de coke dans le cube, qui est l'un des sclaires passifs créés pour le cas particulier de la combustion. Les particules de coke ont un diamètre de 0.005 mètres au début de l'expérience, puis se réduisent plus fortement sur les bords inférieurs du trapèze pour atteindre un diamètre d'environ 0.0038 mètres au bout de quelques instants (moins d'une minute en moyenne pour les points les plus excentrés des parties latérales du trapèze). Nous observons que les particules de coke diminuent de volume plus intensément sur les coins de la base du trapèze, mais de facon inhomogène. C'est à dire que la réduction de volume ne se fait pas symmétriquement des deux coté du trapèze. Ceci est conforme aux attentes de l'industriel qui s'appuyait sur une étude CFD mentionnée plus avant.

Les deux graphiques suivants montrent l'évolution du nombre de particules de coke par unité de masse (chi coke) et l'évolution du diamètre des particules (d coke) en fonction du temps. Nous observons que le nombre de particules de coke par unité de masse augmente au cours du temps, et ce très rapidement au début, puisqu'au bout de 10 secondes il y a deux fois plus de particules par unité de masse qu'au temps initial. L'évolution est ensuite plus lente à partir de 20 secondes, temps auquel chi vaut 2.4.
L'évolution du diamètre des particules de coke suit la même tendense que l'évolution du chi mais à la baisse. C'est à dire que le diamètre des particules de coke diminue très rapidement en début de simulation (le diamètre des particules de coke chute à 3.8mm alors qu'au début la valeur est de 5mm).
Les deux graphiques sont à rapprocher. En effet lorsque les particules de coke réduisent leur volume, elles sont plus nombreuses pour une même masse donnée. Cela revient à observer que les particules de coke perdent de la masse lors de la réduction.



        
Le graphique suivant a été obtenu avec le programme Paraview qui permet de visualiser les résultats des simulations numériques obtenues avec Neptune.
Nous observons qu'au bout de 5 secondes, le nombre de particules de coke par unité de masse (chi coke) n'est pas réparti de facon homogène sur le tas de braises. En effet, l'echelle nous informe que la densité de coke est plus importante sur la partie driote de la base du trapèze. Cela montre l'évolution inhomogène de la coke.






Nous avons cherché à tracer l'évolution spatiale du paramètre chi au même temps (5 secondes) au niveau de la forte densité de coke. Nous observons une répartition qui est elle aussi inhomogène à l'intérieur même de l'accumulation de coke.



        



Dans le graphique suivant, nous observons la répartition spatiale du taux de réaction dans le four au temps t=1.635 secondes. Le taux de réaction correspond au nombre de réactions entre le coke et l'air par unité de volume et de temps. Dans l'enceinte du four, le taux de réaction est très faible (de l'ordre de 10e-8). Au niveau de l'accumulation de coke sur la partie droite de la base du trapèze, on observe une augmentation du gamma. Cette observation est à relier à la répartition du paramètre chi, dans le sens où la réaction du coke avec l'air va diminuer la taille de la coke et ainsi augmenter la densité de celle-ci. de plus, nous obserons le caractère inhomogène de la répartition du taux de réaction.





Nous avons également tracé l'évolution de la vitesse moyenne du mélange gazeu suivant l'axe vertical "y" à un instant t=4.74 secondes. Le mélange gazeu a une vitesse importante au niveau du tas de braise et au dessus (0.27 m/s). Il y a des zones de recirculation sur les extrémités du four à bois (zones bleu foncé de valeur négative suivant l'axe ascendant "y").



Immobilité de la phase de coke


Lors de la mise en place de la simulation de la combustion, l'initialisation des fractions volumiques a été fait de façon très précise (forme trapézoïdale comme décrit dans la rubrique "travaux") et le maillage a été construit de manière à être plus raffiné dans cette zone du tas de braise. Ceci pour suivre cette variable primordiale, renseignant à la fois sur l'efficacité de la combustion mais aussi le validité du choix des termes sources de forçage ultérieur pour imposer la vitesse de chute du tas de particules ainsi que la durée totale de la combustion, sensée pouvoir se prolonger toute la nuit. Pour autant cette donnée nous a posé beaucoup de problèmes et apporté finalement assez peu d'informations. En effet lors de la récupération des modules de Neptune nous avons dans un premier temps mis l'accent sur la compréhension du rôle de chacun d'entre eux tout en l'adaptant à notre maillage, ce qui passait notamment par la supression artificielle de la troisième dimension (i.e. non utilisation et non appel de plusieurs routines ou boucles internes comme "interpolation.F"). Ce que nous n'avons pas fait, par commodité en premier lieu,c'est supprimer la présence d'une troisième phase (noté phase 2 dans les fichiers en annexe) correspondant à de la fonte supposée créer du coke en se consumant. Nous avons bien entendu limité l'impact de cette phase en supprimant la plupart des implications de sa présence sur les scalaires et variables users ainsi qu'en initialisant sa fraction volumique au minimum acceptable par les calculateurs. Pour autant, même de manière infime, la fraction volumique de fonte n'est pas nulle au sein du tas de braise, ce qui impacte sur les résultats de la fraction volumique de coke qui n'évolue quasiment pas étant donné qu'une partie du coke est recréé par combustion de la fonte.

Nous avons commencé la réécriture de tous les modules sans cette phase, mais cette réécriture doit être quasi complète en raison de l'emboitements des modules entre eux. Le temps manquant, imputable en grande partie à l'interprétation du "zéro" initialement mentionné, pour mené à bien cette opération, que nous évaluons à plus d'une semaine pour assurer un "nettoyage" totale, ne nous permet donc pas d'utiliser cette variable. Dans l'eventualité d'une reprise de ce travail, cette étape nous apparait comme la plus importante pour le groupe concerné.

Echanges radiatifs


Les modules fournis par M.Konan régissant les échanges termiques au sein du domaine de calcul tiennent uniquement compte du mécanisme de diffusion/convection entre la phase gazeuse et les particules de coke. En outre il est possible de définir dans Edamox, en tant qu'option, les transferts radiatifs inter-particulaire. Cette option est interne à Neptune. Pour autant les transferts radiatifs entre le gaz et les différentes parois (latérale, vitre et conduit d'évacuation) ne sont pas pris en compte alors qu'ils sont pertinents pour se rapprocher des résultats de la CFD obtenus par EDF. Le modèle correspondant le plus pertinemment avec notre étude est le cas d'échange radiatif entre deux surface l'une étant un plan considéré comme infini et l'autre étant une rangée de cylindres de diamètre d_coke et dont les centres sont écartés d'une distance s (on suppose une compacité élevé de l'ordre de 0.55). Dans la littérature fourni par DeWitt et Bergame, on trouve alors l'approximation suivante pour le facteur de receptivité ( fraction de radiation quittant la surface A et étant intercepté par la surface B):

                                      

La théorie des échanges thermiques radiatifs en surface close permet alors d'obtenir le flux net de chaleur par méchanisme radiatif entre le tas de braise et les différentes parois du four en fonction des propriétés physiques de celles-ci et de leur aire. soit la formulation suivante:

                                      

Ce flux va avoir pour effet d'augmenter la température en paroi, la radiation s'additionnant à l'échange convectif/diffusif du gaz avec la paroi. Or l'une des demandes spécifiques de l'industriel est une limitation de la température à l'intérieur de l'insert pour des raisons compréhensibles de résistance des matériaux utilisés, tant pour la vitre que pour les parois latérales et supérieures.

Ainsi il convient de modéliser aussi l'échange convectif entre la paroi et l'air ambiant extérieur qui tient lieu de refroidissement pour l'insert. Cet échange thermique est de type convection libre. Dans le cas particluier d'un echange avec une paroi plane ou incliné, des corrélations empiriques sont données pour calculer le nombre de Nusselt déterminant pour l'obtention du coefficient de transfert thermique moyen:

                                                                                   

avec

                                         et          


Nous avons ensuite ajouté cette modification d'importance dans le module des transferts thermiques "usth12" mais sans le rendre effectif (voir en annexe). En effet en 2D ces transferts n'ont pas de sens physiques réels étant donnés que les flux de chaleurs sont calculés pour des surfaces qui ne sont pas prtises en compte, avec un maillage 2D, par le solveur. Ainsi l'activation de ces méchanismes supplémentaires d'échanges thermiques pourront être utilisés lors de la tridimentionalisation de l'étude, ce que le temps ne nous permet pas de mettre en place.

Conclusion


Les résultats


BEI : Bureau d'Etude Industriel : Modélistation d'un four à bois (sujet 18)
Benoit Biton et Nicolas Gistau