Modélisation de la dégradation de la production de puissance d’une pile à combustible suite aux sollicitations mécaniques Tilda Akiki To cite this version: Tilda Akiki. Modélisation de la dégradation de la production de puissance d’une pile à combustible suite aux sollicitations mécaniques. Autre. Université de Technologie de Belfort-Montbeliard; Université Saint-Esprit (Kaslik, Liban), 2011. Français. ฀NNT : 2011BELF0156฀. ฀tel-00607268฀ HAL Id: tel-00607268 https://tel.archives-ouvertes.fr/tel-00607268 Submitted on 8 Jul 2011 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. N°d’ordre : 156 Année 2011 THESE pour obtenir le grade de DOCTEUR DE L’UNIVERSITE DE TECHNOLOGIE DE BELFORT-MONTBELIARD ECOLE DOCTORALE : SCIENCES POUR L’INGENIEUR ET MICROTECHNIQUES Spécialité : Mécanique et conception par Tilda AKIKI MODELISATION DE LA DEGRADATION DE LA PRODUCTION DE PUISSANCE D’UNE PAC SUITE AUX SOLLICITATIONS MECANIQUES Soutenance prévue le 3 mars 2011 devant le jury composé de : M. Abel CHEROUAT M. Olivier LOTTIN M. Jean-François BLACHOT M. Willy CHARON M. Naïm OUAINI M. Raed KOUTA Rapporteur Rapporteur Examinateur Directeur de thèse Co-directeur de thèse Encadrant 2 Remerciements Je vous présente mon manuscrit de thèse dont la lecture, je l’espère, vous sera agréable. Il est aussi pour moi une occasion pour remercier toutes les personnes qui m’ont soutenue et encouragée et pour montrer ma gratitude envers les établissements qui m’ont accueillie durant mon travail de recherche. Je remercie donc l’Université de Technologie de Belfort-Montbéliard et spécialement l’institut FCLAB pour m’avoir accueillie maintes fois chez eux. Je remercie aussi l’Université Saint Esprit de l’USEK, établissement au sein duquel je travaille, pour m’avoir facilité la tâche en m’assurant des conditions favorables pour la recherche. Je tiens à remercier très fortement mon directeur de thèse en France Professeur Willy Charon pour sa rigueur scientifique, sa très haute qualité d’encadrement et son sens de justesse. Je remercie aussi mon directeur de thèse au Liban Professeur Naïm Ouaini pour sa présence et sa volonté à résoudre les problèmes rencontrés. Son aide m’est très précieuse. De plus, je remercie mon encadrant Dr. Raed Kouta pour ses conseils et ses encouragements depuis le début de ma thèse ainsi que Dr. Gilbert Accary pour ses compétences scientifiques et son engagement inconditionnel qui m’est très cher. Je n’oublie surtout pas mon amie Marie-Christine pour son aide, son travail d’équipe et sa présence ainsi que Dana et Béatrice qui ont rendu mes séjours en France si agréables. Leur amitié me marquera toujours. Pour mes parents, ma sœur, ma famille et mes amis, mes remerciements ne seront jamais à la hauteur de mes sentiments, de ma gratitude et de mon amour. Pour mes enfants, pour mon mari. Pour tantie … Tilda KARKOUR AKIKI 3 Table des matières Introduction 0.1 Cadre de la recherche 0.2 Objectifs de la thèse 0.3 Présentation du mémoire 1 Mécanique de la PEMFC et couplages multiphysiques dans les publications 1.1 PEMFC, opération et performances 1.2 Equations générales de la PEMFC 1.3 Modèles de PEMFC 1.4 Techniques de caractérisation des piles à combustible 1.5 Phénomènes mécaniques influant la production de puissance d'une PEM 2 Les outils pour la modélisation 2.1 Simulateurs systèmes 1D 2.2 Simulation mécanique 2D ou 3D 2.3 Simulation multiphysique 2D ou 3D 2.4 Interface entre logiciels 3 Influence de la porosité et de la perméabilité locales de la couche de diffusion des gaz sur les performances des PEMFC 3.1 Approche générale 3.2 Les modélisations proposées 3.3 Les résultats de l'analyse numérique 3.4 Conclusion 4 Influence sur les performances des PEMFC de la pression locale à l'interface entre la couche de diffusion des gaz et la plaque bipolaire 4.1 Approche générale 4.2 Les modélisations proposées 4 4.3 Les résultats de l'analyse numérique 4.4 Conclusion 5 Influence sur les performances des PEMFC de la résistance de contact entre la couche de diffusion des gaz et la plaque bipolaire 5.1 Introduction 5.2 Le modèle multiphysique global 5.3 Les modélisations de la résistance de contact 5.4 Les résultats de l'analyse numérique 5.5 Conclusion 6 Etude de sensibilité et de variabilité 6.1 Objectifs du chapitre 6.2 Présentation des résultats dans un plan d'expériences 6.3 Analyse de la variance 6.4 Régression multilinéaire 6.5 Incertitudes sur les modèles étudiés 6.6 Conclusion Conclusion générale et perspectives Annexes A Nomenclature B Liste des figures C Liste des tables D Calcul du champ de porosité Références 5 Introduction 0.1 Cadre de la recherche 0.2 Objectifs de la thèse 0.3 Présentation du mémoire 6 0.1 Cadre de la recherche L’action de recherche décrite dans ce mémoire s’est inscrite dans le cadre de la coopération établie entre l'Université Saint-Esprit de KASLIK (USEK) au Liban et l'Université de Technologie de Belfort-Montbéliard (UTBM). Cette coopération comprend les dimensions pédagogie et recherche. La recherche poursuivie tout au long de cette thèse fait partie d’une part du projet Systèmes Mécaniques Adaptatifs (SMA) du laboratoire mécatronique M3M (Méthodes, Modèles et Métiers) de l’UTBM impliqué dans l’institut FCLAB de recherche sur les systèmes pile à combustible et d’autre part du projet de l’équipe de recherche « Modélisation Multiphysique », en cours de constitution, du département Sciences et Technologies de l’USEK. Ce partenariat vise la maîtrise du processus de conception des systèmes mécatroniques et des conditions de leur production industrielle en grandes séries. L’intérêt de FCLAB porte essentiellement sur les systèmes pile à combustible pour le transport en général qu’il soit terrestre, marin ou aérien. Dans le domaine du transport, seules les piles à basse température de type à membrane échangeuse de protons peuvent délivrer le niveau de puissance exigé. Le domaine de la pile à combustible et ses retombés attendues dans l’amélioration des conditions et des conséquences de la production d’énergie présente un intérêt considérable non seulement pour la France mais aussi pour le Liban. 0.2 Objectifs de la thèse La technologie des piles à combustible est en constante évolution. Elle doit encore cependant relever de nombreux défis scientifiques, techniques et économiques avant d’être commercialisée. Un des principaux défis scientifiques et techniques pour la maîtrise de la conception et la production d’une pile à combustible est de connaître son comportement dans son environnement réel d’usage. De par le contexte dans laquelle elle s’est exécutée, cette thèse s’intéresse aux piles à basse température de type à membrane échangeuse de protons (PEMFC est l’abréviation anglosaxonne que nous conservons dans ce mémoire). 7 Les PEMFC font l’objet de nombreuses recherches pour augmenter leurs performances et diminuer leur coût. Les solutions passent par le développement de nouvelles membranes polymères (électrolytes), par des améliorations sur les cœurs de pile (ensembles électrodesmembrane), par la mise au point de nouveaux catalyseurs, par de nouveaux concepts de plaques bipolaires, métalliques et composites, par un concept innovant de distribution des gaz, ainsi que par une optimisation des systèmes à pile à combustible. La plupart des études se concentrent sur les aspects physico-chimiques des piles à combustible. Par contre cette thèse se propose de mettre en évidence l’influence sur la production d’énergie des sollicitations mécaniques statiques, dynamiques voire thermiques (serrages, vibrations, frottements, …) comme phénomènes couplés relevant du domaine multiphysique (interactions fluide-structure, électrique …). En effet ces influences sont encore largement ignorées bien qu’elles soient de plus en plus ressenties comme déterminantes. Cette thèse s’est donc intéressée aux méthodes et modélisations pour l’étude de cette problématique. Dans la modélisation des problèmes relatifs aux piles à combustible, une distinction est souvent faite entre les modélisations physiques (mécaniques, électrochimiques) et les modélisations multiphysiques. Ainsi, les études s’appuyaient sur des approches de 3 types : • approche mécanicienne pour étudier les déformations possibles à l’intérieur d’une pile à combustible (PAC), • approche électrochimique pour incorporer la dégradation de la tension de la pile en ajoutant un paramètre de vieillissement au modèle électrochimique généralisé en régime permanent et en faire une Analyse des Modes de Défaillance, de leurs Effets et de leur Criticité (AMDEC) [Fowler 2002 et Kundu 2006], ou • approche multiphysique [P. Zhou 2006 et 2007, Y. Zhou 2009, Zhang 2006, Nitta 2007, Fekrazad 2007, Hottinen 2007, Sun 2005, Al-Baghdadi 2007 et 2009] qui offre une description complète de la dynamique, du transfert de chaleur et du courant généré. Dans ce contexte, nous savons qu’une pile est soumise à différentes sollicitations dynamiques dans son environnement réel d’usage d’où la nécessité de prendre en compte des couplages multiphysiques dans la modélisation de la pile. 8 L’étude présentée dans ce mémoire porte sur cette dernière catégorie de modélisations, pour lesquelles un couplage multiphysique est primordial dans le bilan de performance énergétique de la pile. En effet, sous une sollicitation mécanique comme la compression, des propriétés physiques de la pile comme la porosité de la couche de diffusion de gaz (GDL est l’abréviation anglo-saxonne que nous conservons aussi dans ce mémoire), la pression des gaz à l’interface entre la couche de diffusion de gaz et la plaque bipolaire (PB), et la résistance de contact entre la couche de diffusion de gaz et la plaque bipolaire présentent de fortes variations. Ces comportements spécifiques sont à l’origine de la modification des courbes de densité de puissance et donc de la performance de la pile. Au cours des quelques dernières années, quelques études expérimentales et numériques ont été consacrées à l’étude de l’interaction entre les comportements mécaniques et les couplages multiphysiques et ont contribué à une meilleure compréhension de la pile. Le travail présenté dans ce mémoire est la continuation des études numériques basées sur l’approche multiphysique qui consiste à résoudre les équations de Navier-Stokes, de Darcy, de Maxwell-Stefan et de la mécanique des fluides. La prise en compte des effets 3D permet aussi une description plus réaliste du comportement de la pile. Un des objets de ce travail a donc été de développer des approches pour évaluer les effets des sollicitations mécaniques sur la performance et le rendement des PEMFC. 0.3 Présentation du mémoire Le premier chapitre de ce mémoire présente l’état des connaissances des propriétés des piles à combustible ainsi qu’un inventaire des différentes études antérieures, permettant de mieux situer notre contribution. Dans le deuxième chapitre, les outils logiciels concernés par la thèse sont présentés en fonction des objectifs du sujet de thèse et de l’institut FCLAB. Le code de calcul 2D, pour la prise en compte des champs de porosité locale et de perméabilité locale de la GDL dans la modélisation de la production de puissance de la PEMFC, est décrit au troisième chapitre. 9 Dans le quatrième chapitre, un modèle 3D de mécanique des fluides permet d’étudier l’influence de la compression de la pile sur la pression locale d’interface entre la couche de diffusion des gaz et la plaque bipolaire. L’influence de cette pression locale d’interface sur le comportement de la PEMFC est déduite des courbes de densité de courant et de puissance. Le code de calcul 2D du chapitre 3 est ensuite complété pour prendre en considération un champ local de résistance de contact entre la GDL et la PB. Enfin, dans le sixième chapitre, les résultats des études menées précédemment sont rassemblés selon une approche issue des plans d’expériences. Cette approche permet d’effectuer une étude de sensibilité et une modélisation par surface de réponse. De plus, dans certains cas, une analyse est effectuée autour de l’influence des incertitudes des paramètres de conception d’une pile. Ces analyses permettent de se projeter dans des perspectives de recherche de compromis entre les coûts de production et les performances élevées. 10 Chapitre 1 Mécanique de la PEMFC et couplages multiphysiques dans les publications 1.1 PEMFC, opération et performances 1.2 Equations générales de la PEMFC 1.2.1 Eléments de thermodynamique 1.2.2 Cinétique des réactions 1.2.3 Transport des charges 1.2.4 Transport des espèces non chargées 1.3 Modèles de PEMFC 1.3.1 Modèle à une dimension 1.3.2 Différents modèles de la littérature 1.4 Techniques de caractérisation des piles à combustible 1.4.1 Techniques de caractérisation électrochimiques In Situ 1.4.2 Techniques de caractérisation électrochimiques Ex Situ 1.5 Phénomènes mécaniques influant la production de puissance d'une PEM 1.5.1 Composants de la pile dans les publications 1.5.2 Tableau comparatif de modèles relatifs au sujet de thèse 1.5.3 Le système pile Le sous-système PAC Le sous-système de gestion thermique Le sous-système de traitement et de fourniture du combustible Le sous-système de l’électronique de puissance Le sous-système de gestion de l’air 1.5.4 Synthèse et conclusion de l’étude bibliographique 11 L’idée de ce premier chapitre est de partir des modélisations issues de la littérature, d’identifier ensuite les paramètres qui sont soumis à la mécanique et à ses interactions multiphysiques et finalement d’introduire les phénomènes mécaniques couplés influençant les performances et pourquoi pas plus tard la durabilité. En ce sens ce chapitre présente une sorte de synthèse de l’état des connaissances sur les piles à combustible de type PEMFC. Les équations régissant les PEMFC sont d’abord présentées. Cela est nécessaire pour pouvoir identifier les équations qui pourraient être affectées par les phénomènes mécaniques couplés. Ensuite, une section est dédiée aux différents modèles de la littérature, là aussi pour pouvoir agir et intégrer les couplages mécaniques. Il ne s’agit pas d’une présentation détaillée des différents composants d’une PEMFC et de leurs fonctions. Pour cela, le lecteur pourra se reporter à la référence [Biyikoglu 2005]. Il est également intéressant de décrire les techniques de caractérisation des piles à combustible pour examiner si elles peuvent donner des renseignements sur l’influence des couplages mécaniques. Les modélisations sont présentées d’après l’excellent livre d’O’Hayre, Cha, Colella, et Prinz [O’Hayre 2006] dont il conserve la structure dans la mesure où les informations sont pertinentes pour la thèse. Dans ces conditions, les différentes sections sont éventuellement concentrées sur les PEMFC et complétées par des informations en provenance d’autres publications. La première partie du livre porte sur les principes des Piles à Combustible (PAC). Le chapitre 1 définit la pile à combustible, son mode d’opération, sa performance… Le chapitre 2 rappelle les principes de la thermodynamique : l’énergie interne notée U, l’enthalpie H, l’énergie libre de Helmholtz F, l’énergie libre de Gibbs G et tire l’équation de Nernst. Le chapitre 3 parle de la cinétique des réactions, de l’énergie d’activation et tire l’équation de Butler-Volmer pour enfin déduire les effets de la concentration des réactifs, de la barrière d’activation, de la température et des réactions de site dues à des électrodes rugueuses, sur la 12 performance. A partir de l’équation de Butler-Volmer, une équation simplifiée de la cinétique d’activation est tirée : c’est l’équation de Tafel. Le chapitre 4 explique le transport des charges à l’intérieur d’une pile à combustible, dans l’électrolyte principalement. Le chapitre 5 décrit le transport des espèces non chargées à l’intérieur d’une pile à combustible, dans les électrodes par diffusion et dans la structure de flux par convection. Le chapitre 6 propose différents modèles de PAC. Le chapitre 7 présente les techniques de caractérisation des piles à combustible. La seconde partie du livre explicite la technologie des piles à combustible. Le chapitre 8 présentant les différents types de PAC, ne nous concerne pas directement. Les chapitres 9 et 10 présentent les sous-systèmes inclus dans la conception d’une PAC : la PAC elle-même, le sous-système de gestion thermique, le sous-système de livraison/traitement du combustible, le sous-système de l’électronique. Finalement, le chapitre 11 quantifie l’impact des PAC sur l’environnement. 1.1 PEMFC, opération et performances Une pile à combustible transforme l’énergie chimique stockée dans un carburant (l’hydrogène) en énergie électrique. En produisant du courant, la PAC n’est pas consommée, c’est la différence primordiale entre une PAC et une batterie. La figure 1.1 illustre les éléments de modélisation d’une cellule de PAC et les étapes dans la génération électrochimique de l’électricité : Transport des réactifs, réactions électrochimiques, conduction ionique et électronique, évacuation des produits des réactions. 13 Fig. 1.1 – Eléments de modélisation d’une cellule [O’Hayre 2006 - 0]. Notons que les PEMFC fonctionnent à une température avoisinant 80°C. Du côté anode, un flux d’hydrogène pénètre dans le serpentin, diffuse à travers la GDL poreuse et subit une réaction d’oxydation sur le catalyseur qui sert à accélérer les réactions. H 2 → 2H + +2e- (1.1) Du côté cathode, un flux d’oxygène humidifié pénètre dans le serpentin des plaques bipolaires, diffuse à travers la GDL poreuse et atteint le catalyseur qui recouvre les électrodes où une réaction de réduction se produit. 1 O 2 +2H + +2e- → H 2 O 2 (1.2) L’eau produite par la réaction ci-dessous doit être évacué de la pile afin de ne pas l’inonder. Dans les électrodes, le transport de masse (réactifs et produits) se fait par diffusion qui est provoquée par un gradient de concentration. Le transport d’ions H+ à travers l’électrolyte, le Nafion, représente une importante perte résistive, ce qui réduit la performance de la PAC, de même qu’un faible transport de masse (pauvreté en réactifs ou accumulation des produits) dans la couche du catalyseur. De plus, l’énergie d’activation nécessaire à convertir les réactifs en produits fait diminuer la tension de sortie de la pile. 14 Traditionnellement, la performance d’une PAC est représentée par sa caractéristique courant-tension. En pratique, la tension de sortie de la PAC, notée V, est moindre que la tension idéale prédite par la thermodynamique, notée E thermo (voir figure 1.2). Fig. 1.2 - Caractéristique d’une PAC. La courbe de la figure 1.2 est influencée par trois types de perte : • Pertes dues à la barrière d’activation (réactions électrochimiques) ηact • Pertes ohmiques (conduction ionique et électronique) ηohmic • Pertes dues aux concentrations des réactifs ηconc Ainsi, la tension de sortie d’une pile est : V = E thermo - ηact - ηohmic - ηconc 1.2 (1.3) Equations générales de la PEMFC 1.2.1 Eléments de thermodynamique La thermodynamique est l’étude de la transformation d’une forme d’énergie en une autre. Quatre potentiels thermodynamiques sont reliés entre eux par les équations suivantes : 15 F = U – TS (1.4) H = U + pV (1.5) G = U – TS + pV (1.6) avec • U l’énergie interne du système1 en J/mol • H l’enthalpie2 en en J/mol • F l’énergie libre de Helmholtz3 en J/mol • G l’énergie libre de Gibbs4 en J/mol • T la température en K • p la pression en Pa • S l’entropie5 en J/mol.K • V le volume en m3 Pour une PAC, la tension réversible est la tension produite à l’équilibre thermodynamique. Dans les conditions standards, cette tension est notée E thermo . Dans des conditions nonstandards, la tension réversible varie avec : • La température • La pression 1 L’énergie interne U d'un système thermodynamique correspond à l’énergie associée au mouvement microscopique et d’interaction entre les particules aux échelles atomique et microscopique [O’Hayre 2006]. 2 La quantité de chaleur maximale mise en jeu pendant la transformation à pression constante (isobare) d'un système thermodynamique est donnée par l’enthalpie H [O’Hayre 2006]. 3 L’énergie libre F (appelée aussi ‘énergie libre de Helmholtz’) est, en thermodynamique, une fonction d'état dont la variation permet d'obtenir le travail utile susceptible d'être fourni par un système thermodynamique fermé, à température constante [Wikipédia]. 4 Le travail électrique maximal pouvant être obtenu d’une PAC opérant à température et pression constantes est donné par la variation de l’énergie libre de Gibbs de la réaction électrochimique [Fuel Cell Handbook - 1]. 5 L’entropie peut être interprétée comme la mesure du degré de désordre d'un système au niveau microscopique [Fuel Cell Technology Handbook]. 16 • La concentration. Pour la PAC H2 - O2, si la température augmente de 100 degrés, la tension de la pile diminue d’environ 23 mV ( Δsˆ est généralement négative pour les PAC) [O’Hayre 2006 – 1] : Δsˆ ⎛ dE ⎞ ⎜ ⎟ = ⎝ dT ⎠p nF (1.7) • s l’entropie du système • n le nombre d’électrons • F la constante de Faraday. La variation de la tension réversible E avec la pression est liée au changement de volume de la réaction. Si ce changement est négatif (moins de moles de gaz sont générées par la réaction que consommées), alors la tension de la pile augmentera avec l’augmentation de pression : ⎛ dE ⎞ Δvˆ ⎜ ⎟ =nF ⎝ dp ⎠T (1.8) L’équation de Nernst suivante donne la tension de sortie de la pile en fonction de l’activité chimique des espèces en présence: RT ∏ a produits V = E thermo ln vi nF ∏ a réactifs vi (1.9) • E thermo la tension électrochimique réversible à l’état standard (25°C) • R la constante des gaz parfaits • vi le coefficient stœchiométrique qui représente le coefficient affecté à une espèce dans l'équation chimique • a l’activité chimique de chaque espèce et qui représente sa concentration active dans le mélange. Pour une PEMFC, l’équation de Nernst peut aussi être écrite, de façon plus mécanique, en fonction de la pression des gaz en présence : 17 V(p) = E thermo - p H 2O RT ln 1 nF p H2 .pO2 (1.10) 2 Bien que l’équation de Nernst contienne la variable température, l’équation de Nernst ne tient pas compte complètement de l’effet de la variation de la température sur la tension réversible. C’est pourquoi, pour une température arbitraire T, l’équation de Nernst devient : V(T,p) = E T - p H 2O RT ln 1 nF p H2 .pO2 (1.11) Δsˆ (T - T0 ) nF (1.12) 2 avec E T = E thermo + Ainsi, la tension électrochimique réversible dépend de la température et l’activité chimique dépend des pressions partielles des réactifs et des produits de la réaction. Les dépendances des paramètres physiques que sont la température et la pression sont maintenant évidentes. Aussi, il existe une dépendance entre non seulement les pressions partielles des réactifs et produits mais aussi entre la pression de serrage et la température dans l’AME (Assemblage Membrane Electrode) de la pile. L’augmentation de la pression de serrage change la distribution de température interne et influe sur la tension de sortie de la pile [Fekrazad 2007]. D’autre part, l’efficacité réelle d’une PAC est donnée par la formule : ε réelle = ε thermo × ε tension × ε carburant • (1.13) ε thermo l’efficacité thermodynamique réversible ; elle reflète que l’enthalpie contenue dans une pile ne peut pas être exploitée entièrement pour fournir du travail utile. • ε tension l’efficacité de tension de la pile due aux pertes cinétiques ; elle est égale au rapport de la tension d’opération de la pile et de la tension thermodynamique réversible. • ε carburant l’efficacité du carburant due aux pertes de l’utilisation du carburant. D’où : 18 ⎛ Δgˆ ε réelle = ⎜ ⎜ Δhˆ ⎝ HHV ⎞ ⎛ V ⎞⎛ 1 ⎞ ⎟⎟ ⎜ ⎟⎜ ⎟ ⎠ ⎝ E ⎠⎝ λ ⎠ • G l’énergie libre de Gibbs • H l’enthalpie • λ le facteur stœchiométrique • HHV pour “higher heating value”6. (1.14) Cette efficacité thermodynamique diminue quand la température augmente. Remarque : Si on fournit à une PAC 1,5 fois plus de carburant que nécessaire pour une utilisation de 100% de carburant, on dit alors que la PAC opère à 1,5 fois stœchiométrique. 1.2.2 Cinétique des réactions Toutes les réactions électrochimiques impliquent le transfert de charges (électrons) entre une électrode et une espèce chimique adjacente à la surface de l’électrode. Le courant produit par une PAC (nombre d’électrons par unité de temps) dépend du taux de la réaction électrochimique (nombre de réactions par unité de temps). La densité de courant j, qui est le courant par unité de surface, est définie par : i = ∫ jds (1.15) A • A la surface en cm2 • i le courant en ampères. Pour que les réactifs soient convertis en produits, ils devraient vaincre une barrière d’énergie nommée énergie d’activation. Cette énergie d’activation détermine le taux de réaction. Le taux d’une réaction υ est la différence entre le taux de la réaction directe et du taux de la réaction inverse. 6 La valeur calorifique de l’hydrogène est la quantité générée par la combustion complète d’une mole d’hydrogène [Barbir 2005 - 1]. 19 υ=c fe * R 1 • -ΔG1++ RT * p 2 -(ΔG1++ - ΔG rxn ) RT -c f e (1.16) * c*R et c p sont les concentrations de surface des réactifs et des produits respectivement ΔG1++ et ΔG 2++ sont les barrières d’activation de la réaction directe et de la • réaction inverse, avec ΔG rxn = ΔG1++ − ΔG 2++ • f1 et f 2 sont les taux de dépérissement en produits et en réactifs respectivement. La densité de courant direct j1 et le taux de réaction υ sont reliés par : * R 1 -ΔG1++ RT j1 = nFc f e (1.17) La densité de courant inverse j2 et le taux de réaction υ sont reliés par : * p 2 -(ΔG1++ - ΔG rxn ) RT j2 = nFc f e (1.18) Une partie de la tension de la PAC est aussi sacrifiée pour diminuer l’énergie d’activation qui constitue une barrière pour la conversion des réactifs en produits. Cette tension perdue est nommée surtension d’activation et elle est notée ηact . La relation entre la densité de courant et la surtension d’activation est exponentielle, elle est décrite par l’équation de Butler-Volmer suivante : αnFηact RT - j = j0 (e -(1-α)nFηact e RT ) • j la densité de courant • j0 la densité de courant d’échange • α le coefficient de transfert de charges • ηact la surtension d’activation. (1.19) A l’équilibre, on a j1 = j2 = j0 (1.20) 20 A noter que l’équation de Butler-Volmer n’est autre que la différence entre les densités de courant direct j1 et inverse j2 , où α exprime comment le changement du potentiel électrique dans l’interface d’une réaction influe sur la barrière d’activation directe et inverse. L’augmentation de j0 améliore la « performance cinétique » d’une PAC. Il existe quatre moyens pour augmenter j0 : • Augmenter la concentration des réactifs c*R • Diminuer la barrière d’activation ΔG1++ • Augmenter la température T • Augmenter le nombre des réactions de site (augmenter la surface de l’interface par exemple). A noter que la densité de courant d’échange j0 dépend de la construction, du matériau et de la température. Les coefficients de transfert de charges (à l’anode et à la cathode) sont gouvernés par les processus de transfert d’électrons à l’interface électrode-électrolyte [Noren 2005]. De plus, les coefficients de transfert dépendent du choix du catalyseur [O’Hayre 2006 - 2]. Ainsi donc, dans l’équation de Butler-Volmer (1.19), il faudra définir une densité de courant d’échange qui dépend du champ de température et des coefficients de transfert qui sont fonctions du catalyseur. Dans le cas où ηact est très important, l’équation de Butler-Volmer peut être simplifiée et donne l’équation de Tafel. La surtension d’activation vaut alors : ηact = a + b log j • a=- • b= (1.21) RT ln j0 αnF RT αnF la pente de Tafel. Dans le cas où ηact est très petit, nous obtenons après développement par série de Taylor des termes exponentiels ( e x ≈ 1+x ): j = j0 nFηact RT (1.22) 21 Pour les PAC à basse température, le platine est utilisé comme catalyseur alors que les PAC à haute température utilisent des catalyseurs à base de Nickel ou de céramique. Un catalyseur est efficace s’il a une bonne activité ( j0 important), une bonne conductivité et s’il est stable. A noter que la cinétique de l’hydrogène est facile alors que celle de l’oxygène donne une importante surtension d’activation à faibles températures. 1.2.3 Transport des charges Pour une PAC, les ions et les électrons se déplacent sous l’effet de trois forces : • Les forces électriques de conduction, • Les forces chimiques de diffusion, • Les forces mécaniques de convection. Dans une PEMFC, les protons et les électrons s’accumulent à l’anode alors qu’ils sont consommés à la cathode. Cette accumulation/déplétion d’électrons aux deux électrodes crée un gradient de tension qui transporte les électrons de la réaction d’oxydation de l’anode vers la cathode à travers la charge électrique extérieure et y crée donc un courant électrique. Dans l’électrolyte, l’accumulation/déplétion de protons crée un gradient de tension et un gradient de concentration. Ces deux gradients transportent les protons de l’anode à la cathode. Notons que l’effet du gradient de tension est plus important que celui du gradient de concentration. Le transport de charges provoque une perte de tension pour la PAC, notée ηohmic . ηohmic = iR ohmic = i (R elec + R ionic ) • Rohmic la résistance ohmique (1.23) de la PAC (électrodes, électrolyte, interconnexions) • Relec la résistance électronique • Rionic la résistance ionique. 22 La résistance de l’électrolyte domine les autres résistances c’est pourquoi l’électrolyte doit avoir une grande conductivité ionique, une faible conductivité électronique et une grande stabilité (dans des environnements d’oxydation et de réduction). On a aussi : ηohmic = j (ASR ohmic ) (1.24) ASR ohmic est la résistance spécifique de surface. Elle vaut : ASR ohmic = A PAC R ohmic (1.25) Sachant que : R= L Aσ (1.26) • L l’épaisseur du conducteur • A la surface • σ la conductivité. Nous pouvons écrire aussi: ASR = L σ (1.27) Ainsi pour diminuer ASR, on utilise un électrolyte de largeur très mince. Le phénomène mécanique de compression de la pile fait diminuer l’épaisseur de la membrane et celle de la GDL. Les résistances ioniques et électroniques sont réduites respectivement, d’où une diminution des pertes ohmiques totales. De plus, la compression et le serrage de la pile réduit la résistance de contact GDL/PB et par suite les pertes des différentes interconnexions. D’autre part, la conductivité d’un matériau est donnée par la formule : σ i = ( z i F)ci u i (1.28) • ci la concentration molaire des porteurs de charge • ui la mobilité des porteurs de charge dans le matériau 23 • zi le nombre de charges des porteurs de charge. Dans le cas d’une PEMFC, on utilise le Nafion comme électrolyte. La conductivité du Nafion aux protons est dominée par la quantité d’eau présente. Une grande quantité d’eau donne une grande conductivité. Ainsi, la conductivité du Nafion peut être modélisée par la quantité d’eau dans la membrane. La quantité d’eau λ dans le Nafion est le rapport du nombre de molécules d’eau sur le nombre d’ions d’acide sulfonique SO3- H + du Nafion. Pour des PEMFC opérant aux alentours de 80°C , le contenu d’eau ou quantité d’eau λ est donnée par : λ = 0.0043 + 17,81aw − 39,85aw2 + 36, 0aw3 pour 0 < aw ≤ 1 (1.29) λ = 14 + 1, 4(aw − 1) pour 1 < aw ≤ 3 (1.30) aw est l’activité de la vapeur d’eau, elle est égale à : aw = pw pSAT • pw (1.31) l’actuelle pression partielle de la vapeur d’eau dans le système à la température d’opération • pSAT la pression de saturation de la vapeur d’eau du système à la température d’opération. Les pressions partielles dépendent de la température d'opération. Il est évident que cette équation locale est utilisée comme équation globale en supposant que la température d'opération est la même partout. En réalité une étude thermique montrerait que cette température est variable spatialement et que cette variation spatiale dépend de la géométrie de la pile et de la température extérieure. Il faudrait donc définir une activité moyenne de la vapeur d'eau qui pourrait dépendre de facteurs de formes géométriques et du champ de température extérieur. La conductivité du Nafion dépend considérablement de la quantité d’eau λ mais aussi de la température. 24 Avec ⎡ 1 ⎞⎤ ⎛ 1 σ(T,λ) = σ 303K (λ) exp ⎢1268 ⎜ - ⎟⎥ ⎝ 303 T ⎠ ⎦ ⎣ (1.32) σ 303K (λ) = 0,005193λ - 0,00326 (1.33) Cette température d’opération varie dans l’espace et le temps. Ainsi, il faudrait définir ici aussi une valeur moyenne de la quantité d’eau λ qui dépend de la température. Il reste à noter qu’il existe deux phénomènes de mouvement des molécules d’eau dans le Nafion: • Les protons emportent dans leur mouvement des molécules d’eau de l’anode à la cathode, c’est la traînée électro-osmotique (electro-osmotic drag). • Un mouvement de diffusion inverse (back diffusion) transporte les molécules d’eau de la cathode à l’anode. L’équation du flux molaire de l’eau est donc : J H2O = J H2O,drag + J H2O,back diffusion J H2O = 2n SAT drag j λ ρdry dλ D λ (λ) 2F 22 M m dz (1.34) (1.35) • J H2O le flux molaire de l’eau (mol/s.m2) • n drag le nombre de molécules d’eau accompagnant chaque proton • n SAT drag ≈ 2,5 • λ la quantité d’eau • ρ dry la densité sèche du Nafion (kg/m3) • Mm la masse équivalente du Nafion (kg/mol). • Dλ la diffusivité de l’eau dans la membrane (m2/s) D λ n’est pas constant, il dépend de la quantité d’eau λ qui, comme nous venons de l’écrire, pourrait dépendre du champ de température. Ce paramètre physique qu’est la température est influencé par des sollicitations mécaniques d’où un couplage multiphysique évident. 25 Dans la suite de ce chapitre, nous soulignerons des coefficients qui sont susceptibles d’être influencés par les couplages avec les phénomènes mécaniques. Ces coefficients primaires influenceront alors d’autres variables à travers les équations des modèles étudiés. Dans l’article [Seong 2006], une nouvelle équation d’état est établie afin de décrire l’activité de l’eau dans l’interface membrane polymère/eau dans une PAC à polymère solide. Le modèle développé est basé sur l’équation d’état PHSC (Perturbed Hard-Sphere-Chain)7 en y introduisant une perturbation obtenue du potentiel généralisé de Lennard-Jones (GLJ). Cette perturbation tient en fait pour la force d’interaction entre les particules. L’équation d’état de van der Waals est : P ⎛ P ⎞ =⎜ + ⎟ ρkT ⎝ ρkT ⎠reference ⎛ P ⎞ ⎜ ⎟ ⎝ ρkT ⎠ perturbation • P • ρ= • k la constante de Boltzmann • T la température. (1.36) la pression N V N étant le nombre de molécules et V le volume La perturbation est de la forme : ⎛ P ⎞ ρU = ⎜ ⎟ ⎝ ρkT ⎠ perturbation 2kT (1.37) • U = 4π ∫ W(r)r 2 dr l’énergie de perturbation • W(r) le potentiel généralisé de Lennard-Jones • r la distance entre les centres de deux molécules adjacentes. La pression P dépend du volume V qui est directement influencé par la compression de la pile. De plus, l’énergie de perturbation U dépend de P et de V. Donc dans le modèle étudié, 7 Dans la théorie des PHSC pour la représentation d’un modèle de polymères, la molécule est constituée de chaînes formées de corps durs, tangents et joints les uns aux autres ou segments. Le modèle a 3 paramètres caractéristiques : Le nombre de « hard-sphere » par chaînes, la taille des segments et l’énergie d’interaction. 26 une compression de la pile fera varier le volume de gaz et en conséquent la pression et l’énergie de perturbation. D’où un couplage avec le phénomène mécanique de compression. Toute variation de la perturbation modifie l’activité de l’eau dans la membrane et donc la performance de la pile Les résultats des simulations ont montré que l’activité de l’eau augmente quand le poids équivalent du polymère diminue (le poids équivalent du Nafion est le poids du Nafion par mole du groupe acide sulfonique). Dans l’article [Trabold 2006], l’auteur a montré que l’eau tend à s’accumuler dans les courbures des serpentins aux deux électrodes. Ceci a été démontré expérimentalement par la méthode de neutronographie qui est basée sur l’atténuation de la radiation suivant une loi exponentielle durant le passage à travers la matière. Les mesures ont été faites in situ (voir § 1.4.1). L’article [Chen 2005] étudie un humidificateur de membrane qui utilise les gaz d’échappement de la pile pour humidifier l’air entrant dans la pile. Les auteurs ont établi un modèle thermodynamique et ont ensuite conduit des simulations en régime stationnaire et en régime dynamique pour optimiser la conception de l’humidificateur en question. Le coefficient de diffusion de la membrane est : D w = Dλ e2416(1/303-1/Tm ) (1.38) où Tm est la température de la membrane. Le coefficient D λ est déterminé de façon empirique et dépend du contenu d’eau λ m de la membrane. Le coefficient de diffusion D w dépend de la température qui est asservie par un système de refroidissement et reste aux alentours de 80˚C. Mais cette température varie quand même spatialement. De plus, D w est aussi influencé par la compression de la PAC. L’activité de l’eau est définie par : 27 ak = p k,v (1.39) p k,sat • p k,v • p k,sat la pression de saturation de vapeur contenu dans le volume de la pression partielle de vapeur contrôle. Le transfert de chaleur est défini par le coefficient U: U= kmN Dh • km la conductivité thermique de la membrane • N le nombre de Nusselt8 • Dh le diamètre hydraulique du canal. (1.40) Le coefficient de diffusion et la conductivité thermique de la membrane dépendent de la température. L’activité de l’eau dépend de la pression. Or, comme nous l’avons vu avant, ces variables sont influencées par les couplages avec des phénomènes mécaniques comme les vibrations et le serrage des plaques bipolaires et ne doivent donc pas être gardées constantes lors de la modélisation. De plus, afin d’éviter la déshydratation ou l’inondation de la pile, il est crucial de contrôler l’humidité de la membrane. 1.2.4 Transport des espèces non chargées Le transport des espèces non chargées, c’est-à-dire le transport de masse, concerne la fourniture des réactifs et l’évacuation des produits dans une PAC. Ce type de transport se fait grâce à des forces de diffusion et de convection et non pas à des gradients de tension. Dans les électrodes, le transport de masse se fait par diffusion qui est provoquée par un gradient de concentration. La densité de courant limite jL correspond au point où la 8 Le nombre de Nusselt représente le rapport entre le transfert thermique par convection et le transfert par conduction. [Barbir 2005 - 2]. 28 concentration des réactifs tombe à zéro dans la couche du catalyseur. Une PAC ne peut jamais délivrer une densité de courant supérieure à jL : jL = nFDeff c0R δ (1.41) • Deff la diffusivité effective • c0R la concentration des réactifs dans la structure de flux • δ la largeur de l’électrode. La densité de courant limite jL dépend de la diffusivité effective et la largeur de l’électrode qui sont des coefficients primaires, influençables par les phénomènes mécaniques de compression ; d’où le couplage de jL avec les phénomènes en question. De plus, la concentration des réactifs dans la structure de flux varie dans le temps et dans l’espace et ne peut être gardée constante. Un faible transport de masse (pauvreté en réactifs ou accumulation des produits) dans la couche du catalyseur provoque une importante perte dans la performance de la PAC. Ici encore, une nouvelle perte de tension est due à la déplétion des réactifs dans la couche du catalyseur, elle est notée ηconc . La concentration affecte la performance d’une PAC à travers l’équation de Nernst : ⎛ RT 1 ⎞ ⎛ RT 1 ⎞ RT c0R ηconc = E 0Nernst - E*Nernst = ⎜ E thermo ln 0 ⎟ - ⎜ E thermo ln ⎟ = ln nF c R ⎠ ⎝ nF c*R ⎠ nF c*R ⎝ (1.42) La concentration des réactifs dans la couche du catalyseur, notée c*R , n’est pas constante tout le long du catalyseur et sa valeur varie en chaque point. La concentration c0R est la concentration « en gros » des réactifs. Pour j < jL , nous obtenons: ηconc = RT j ln L nF jL - j (1.43) 29 La concentration affecte aussi la performance d’une PAC à travers la cinétique de la réaction (équation de Butler-Volmer). Le calcul, déduit de l’équation de Butler-Volmer pour des densités de courant élevées, donne : ηconc = RT j ln L αnF jL - j (1.44) La perte totale due aux concentrations est donc la somme des « pertes de Nernst » et des « pertes de Butler-Volmer » : ηconc = RT j RT j ln L + ln L nF jL - j αnF jL - j ηconc = c ln jL jL - j (1.45) (1.46) c est une constante qui vaut : c= RT ⎛ 1 ⎞ ⎜1+ ⎟ nF ⎝ α ⎠ (1.47) Il est bon de noter que la valeur effective de la constante c est plus grande que celle prédite par l’équation précédente. Dans plusieurs cas, c est obtenu de façon empirique. Dans les structures de flux, le transport de masse se fait par convection, c.à.d. sous l’action d’une force mécanique. Dans une PAC, la convection est donnée par le nombre de Reynolds Re qui caractérise la viscosité du flux. Re = ρVL VL = μ υ (1.48) μ la viscosité cinématique ρ • υ= • ρ la densité du fluide • μ la viscosité dynamique du fluide • V la vitesse caractéristique du flux • L la longueur caractéristique du flux. 30 La viscosité μ dépend très faiblement de la pression mais dépend considérablement de la température. Deux lois sont données : μ ⎛T⎞ ≈⎜ ⎟ μ 0 ⎝ T0 ⎠ n loi simple (1.49) loi de Sutherland (1.50) 1,5 μ ⎛ T ⎞ T0 +S ≈⎜ ⎟ μ 0 ⎝ T0 ⎠ T+S n, μ 0 , T0 et S sont obtenus expérimentalement. Comme la température est un champ variable spatialement et temporellement dans la PAC, il serait juste de penser à remplacer la viscosité μ par un coefficient global qui dépendrait du champ de température. Une différence de pression est requise pour permettre l’écoulement du gaz à travers le canal de la structure. Cette chute de pression est principalement due à la friction entre le fluide et le canal. Ces frictions dépendent de la rugosité des parois des canaux et cet état de surface peut changer avec le temps suite à des corrosions. Cette variation de pression est donnée par la formule suivante : dp 32uμ = dx D 2h (1.51) Cette formule est facilement démontrée pour les canaux de section circulaire : dp 4 = τw dx D et f= (1.52) • D le diamètre du canal • τw l’effort de cisaillement. τw (1.53) 1 ρu 2 2 • f le facteur de friction • u la vitesse moyenne du flux. 31 Re = • et ρuD μ µ (1.54) la viscosité dynamique. f.Re = 16 pour les canaux circulaires (1.55) D’où finalement: dp 32uμ = dx D2 (1.56) Le profil de vitesse et donc la vitesse moyenne varie avec la rugosité des parois des canaux. Pour les autres sections, on définira un diamètre hydraulique. Par exemple, pour une section rectangulaire, on aura : Dh = 4A P (1.57) Où A est la section et P le périmètre. Dans les PAC, les gaz se déplacent le long du canal, ils peuvent être aussi transportés entre le canal et l’électrode. C’est le transport de masse convectif. Ce transport est caractérisé par un coefficient de transfert hm calculé à partir du nombre de Sherwood Sh9 : h m = Sh Dij (1.58) Dh • Dh le diamètre hydraulique • Dij le coefficient binaire de diffusion pour l’espèce i en présence de l’espèce j. Le choix du type d’arrangement des canaux affecte considérablement la valeur des pertes dues au transport de masse. De plus, à cause de la formation d’eau à la cathode, les PEMFC requièrent des flux à grande capacité d’élimination de l’eau. 9 Le nombre de Sherwood représente le rapport entre le transfert total de masse et le transfert par diffusion [Wikipédia]. 32 Il existe différents types d’arrangement de canaux : • canaux parallèles • canaux serpentins • canaux interdigitaux Les structures serpentines ou parallèles-serpentines sont les plus communes. Elles présentent un compromis entre la chute de pression et la capacité d’élimination de l’eau. Les matériaux utilisés dans les plaques bipolaires doivent répondre à plusieurs critères : • Grande conductivité électrique • Grande résistance à la corrosion • Grande conductivité thermique • Faible volume et masse • Facilité de fabrication • Prix adéquat. 1.3 Modèles de PEMFC 1.3.1 Modèle à une dimension Un modèle à une dimension de la PAC est proposé à partir d’une approche basée sur la balance des flux. Nous avons déjà vu que la tension de sortie réelle d’une PAC est donnée par : V = E thermo - ηact - ηohmic - ηconc (1.59) • Ethermo la tension thermodynamique de la PAC • ηact les pertes d’activation dues à la cinétique des réactions • ηohmic les pertes ohmiques des résistances ioniques et électroniques • ηconc les pertes de concentration dues au transport de masse. En utilisant l’approximation de Tafel qui n’est valable que pour j >> j0 , nous obtenons : 33 ⎛ j ⎞ V = E thermo - (a A + b A ln j) - (a C + bC ln j) - (j ASR ohmic ) - ⎜ c ln L ⎟ jL - j ⎠ ⎝ • (a A + b A ln j) - (a C + bC ln j) (1.60) les pertes d’activation de l’anode et de la cathode respectivement • j ASR ohmic • c ln jL jL - j les pertes ohmiques les pertes de concentration. Rappelons que j0 est la densité du courant d’échange et que la constante c dépend de α (Eq. 1.47). Il reste à noter que pour modéliser le comportement réel de la PAC, il faudra ajouter une perte parasite due à un courant de fuite jleak , à des réactions non voulues… tel que : (1.61) jgross = j + jleak • jgross le courant produit aux électrodes de la PAC • j le courant d’opération actuel de la PAC. La tension de sortie réelle devient : ⎛ ⎞ jL V = E thermo - [ a A + b A ln (j + jleak )] - [ a C + bC ln (j + jleak )] - (j ASR ohmic ) - ⎜ c ln ⎟ jL - (j + jleak ) ⎠ ⎝ (1.62) La perte de tension ohmique doit être basée sur la densité de courant j puisque uniquement le courant d’opération de la PAC est conduit à travers la pile. Le courant de fuite jleak est perdu par des réactions secondaires ou par des réactions non électrochimiques aux électrodes. Ainsi, le modèle de base d’une PAC requière quatre paramètres : • α et j0 décrivent les pertes cinétiques • ASRohmic décrit les pertes ohmiques 34 • jL décrit les pertes dues à la concentration. Voici maintenant un modèle plus sophistiqué de la PEMFC, basé sur la balance des flux. Dans ce modèle, la distribution de l’eau dans la pile est critique. La méthode de balance des flux repose sur l’idée « Ce qui rentre doit sortir » : J AH2O J HM2O J CH2O JM j A C H+ = = J H2 = 2J O2 = = = 2F 2 α α 1+α • J AH2 le flux net de H 2 à l’anode • α le coefficient de transfert de charge • JM H2O le flux net de H 2O à travers l’électrolyte... (1.63) Une série de suppositions est nécessaire afin de faciliter le travail de modélisation : • Le transport par convection est ignoré. La convection est typiquement le transport de masse dominant dans une PAC. Cependant, le modèle est un modèle 1D et la convection se fait suivant l’axe qui n’est pas l’axe d’étude. • Le transport par diffusion est ignoré (dans les structures de flux uniquement) • Les pertes ohmiques proviennent uniquement de l’électrolyte • Les pertes par activation proviennent de la concentration d’oxygène à la couche cathode-catalyseur • Les couches catalytiques sont si minces qu’elles se comportent comme des interfaces • L’eau existe à l’état gazeux uniquement (Cette supposition affecte cependant les résultats à la cathode). Equations qui gouvernent chaque domaine : Dans les électrodes, les processus de diffusion de H 2 , O 2 , H 2 O et N 2 doivent être modélisés : 10 Ji = • 10 -pDijeff dx i RT dz Ji (1.64) le flux molaire de diffusion de l’espèce i Le flux molaire de diffusion est Ji = -D.dci/dz or pv=nRT et ci = c.xi d’où l’équation 1.64. 35 • Dijeff la diffusivité effective • xi la fraction molaire de l’espèce i • p la pression totale du gaz à l’électrode tel que : p i = px i Noter que cette équation n’est valable que pour des processus impliquant deux espèces ; sinon il faudra appliquer l’équation de Maxwell-Stefan qui prend en compte des gaz à plusieurs composants. Dans l’électrolyte, on considère le flux de protons et le flux de l’eau. En se rappelant que dans le Nafion, il existe deux types de flux d’eau : la diffusion de retour et la traînée électroosmotique, on a alors: J H 2O = 2n drag j λ ρ dry dλ D λ (λ) 2F 22 M m dz • ρ dry • D λ (λ) la diffusivité de l’eau dans le Nafion. (1.65) la densité sèche du Nafion Le contenu en eau λ n’est pas constant, il est fonction de z. Dans le catalyseur, seule la cinétique de la réaction à la cathode est considérée. D’après Butler-Volmer simplifié, on a : ηcathode = jc0O2 RT ln 4αF j0 cO2 (1.66) 4 représente le nombre de valence pour une molécule d’oxygène. Dans toute la présentation qui précède, l’approvisionnement en hydrogène et oxygène est supposé infini. Cependant, en fonction du taux d’approvisionnement et de consommation, il peut y avoir épuisement des réactifs. Comme le modèle ignore les pertes de concentration à l’anode, seulement la déplétion en oxygène à la cathode est prise en compte (car son effet est plus important). Nous avons rencontré aussi différents types de modèles de la PEMFC dans la littérature. Nous allons les citer dans ce qui suit. 36 1.3.2 Différents modèles de la littérature Dans l’article [Busquet 2004], le modèle est empirique et permet de simuler les courbes V-J des PAC et des électrolyseurs dans le régime stationnaire. Le modèle est particulièrement adapté pour les PAC réversibles qui peuvent être utilisées en pile ou en électrolyseur. Le modèle a 4 degrés de liberté et les régions où les transferts de masse existent ne sont pas modélisées afin de ne pas détériorer la pile. La tension de la pile est : V(J) = E + b ln( J Jd )-2 +( b - Δ) x J 4J d (1.67) Les 4 paramètres (E, Jd, Δ and b) à déterminer dépendent de la température de la PAC et de la pression partielle d’oxygène. Cette influence est de la forme : K1 + K 2T + K 3 ln(p O2 ) (1.68) La méthode des moindres carrés est utilisée pour déterminer E, Jd, Δ and b. Il faut noter ici les couplages possibles avec les phénomènes mécaniques. Par exemple, la pression partielle d’oxygène pO2 est influencée par des connexions desserrées ce qui provoque une diminution de la pression et une augmentation des pertes de transport de masse d’où une dégradation de la performance de la pile. L’article [Xue 2006] établit une révision des différents modèles de la caractéristique V-J d’une PEMFC établis dans la littérature puis classe ces modèles en trois types : • Les modèles physiques formulés à partir de méthodes électrochimiques ou empiriques • Les modèles électriques formulés à partir du circuit électrique équivalent de la pile • Les modèles mathématiques qui sont établis pour convenir avec la caractéristique V- J en utilisant des approches mathématiques. L’auteur de l’article a établi un modèle mathématique unifié pour prédire les caractéristiques V- J des PAC en régime stationnaire et en régime dynamique. La méthode proposée est 37 basée sur la technique des moindres carrés et sur un ensemble d’équations électrochimiques qui représentent les PEMFC. La composante continue de la caractéristique V-J comprend la tension réversible et les pertes ohmiques. La composante transitoire comprend les pertes d’activation et de concentration. La réponse dynamique d’une PAC existe du fait que tout changement du courant n’entraîne pas un changement immédiat de la tension à cause de l’existence d’une double couche chargée dans la pile. L’introduction d’un temps de réponse du premier ordre décrit la réponse dynamique de la pile. La caractéristique V-J d’une PEMFC est : Vcell = E N - Va - Vc - Vohm (1.69) • Vcell la tension de sortie de la pile • EN la tension réversible • Va les pertes dues à l’activation de l’anode et de la cathode (surtensions d’activation) • Vc les pertes dues à la réduction des concentrations des réactifs ou du transport de masse de l’oxygène et de l’hydrogène (surtensions de concentration) • Vohm les pertes ohmiques dues à la résistance de l’électrolyte aux protons et à la résistance du parcours aux électrons (surtensions ohmiques). Les composantes transitoire et continue sont respectivement : Et Vtr = Va + Vc (1.70) Vst = E N - Vohm (1.71) Vcell = Vst − Vtr (1.72) Après un long calcul basé sur la technique des moindres carrés et après introduction d’un temps de réponse du premier ordre, l’équation unifiée est obtenue. 38 Le serrage des vis des plaques bipolaires fait varier l’épaisseur de la membrane et l’épaisseur et la porosité de la GDL. Avec le temps, des fissures peuvent apparaître dans le catalyseur ainsi que des trous dans la membrane. Ces phénomènes influent sur les surtensions ohmiques, d’activation et de concentration. Des couplages existent donc entre, d’un côté, les phénomènes mécaniques et les phénomènes naturels et de l’autre côté, les variables du modèle de pile. Les conséquences sur la performance de la pile sont multiples ; la diminution de l’épaisseur de la membrane réduit sa résistance et diminue les pertes ohmiques [Al-Baghdadi 2007]. L’augmentation de la porosité de la GDL améliore le transport de masse et fait diminuer ses pertes-là [Al-Baghdadi 2007], [Sun 2005], [P. Zhou 2006] et [P. Zhou 2007]. Les fissures dans le catalyseur augmentent sa résistance et baissent sa performance tout comme le délaminage entre la membrane et le catalyseur [Al Baghdadi 2009]. L’article [Tao 2005] introduit un modèle non linéaire basé sur les réseaux de neurones artificiels. Le modèle utilise des données expérimentales comme échantillons d’apprentissage. Les entrées étant la vitesse de l’air et la température d’opération de la PEMFC et les sorties étant la tension et la densité de courant de la pile. L’apprentissage est fait par l’algorithme LMBP (Levenberg-Marquardt Back Propagation). Ici aussi, toute variation de la vitesse de l’air et de la température d’opération de la pile, due à des phénomènes physiques résultant de vibrations mécaniques, modifie la tension de sortie et la densité de courant de la pile. Dans l’article [Fowler 2002], l’auteur ajoute au modèle électrochimique généralisé en régime statique (GSSEM : Generalised steady state electrochemical model for a PEM fuel cell) de nouveaux termes afin de pouvoir prendre en considération le vieillissement de l’assemblage membrane/électrode et d’incorporer ainsi la dégradation de la tension avec le temps dans le modèle. Ces termes sont la perte apparente de l’activité du catalyseur, l’augmentation de la résistance de la pile et la diminution du taux de transport de masse des réactifs. D’après les résultats des expériences, l’auteur a proposé une simple relation linéaire pour définir le taux de dégradation de la conductivité de la membrane et l’activité du catalyseur est considérée à travers un seul paramètre semi-empirique. La tension de sortie de la pile est donnée par : 39 Vcell = E thermo - ηact - ηconc - ηohm (1.73) Des fissures dans le catalyseur baissent la performance de l’activité catalytique. Le délaminage entre la membrane et le catalyseur augmente la résistance entre les différentes couches. L’apparition de trous dans la membrane crée des poches qui inondent la pile. Les surtensions d’activation, de concentration et ohmiques sont les variables impliquées par ces phénomènes physiques. Ici encore, les couplages ne peuvent être négligés dans l’étude de la performance de la pile. [Kundu 2006] présente une catégorisation des causes, modes et effets de la dégradation de la PEMFC. La référence examine aussi les défauts créés dans les membranes recouvertes de catalyseur (CCM catalyst-coated membranes) lors du processus de fabrication. Les facteurs qui influencent la fiabilité opérationnelle de la pile ont pour causes les propriétés des matériaux, la méthode d’assemblage de la pile et du stack ainsi que son intégration dans le système complet. De plus, l’environnement opérationnel (température, humidité, stœchiométrie) et la maintenance exercent une grande influence sur la dégradation de la pile. D’après leurs expériences effectuées sur la pile, les auteurs ont identifié six défauts majeurs sur des CCM : 1. Les fissures sur la surface du catalyseur lors de l’assemblage membraneélectrode (MEA Membrane Electrode Assembly) ou à travers une mauvaise manipulation de la membrane. Ces fissures augmentent la résistance du catalyseur et forment des poches pour la rétention d’eau. Cette inondation réduit le taux des réactions. De plus, les fissures exposent le catalyseur à l’érosion et donc à une perte graduelle de l’activité catalytique. 2. L’orientation des particules est due aux équipements d’élaboration du catalyseur. 3. Le délaminage entre le catalyseur et la membrane entraîne une augmentation de la résistance et crée des poches qui se remplissent d’eau et causent 40 l’inondation de la pile. Aussi, plus de courant ionique sera dirigé vers les zones voisines aux zones délaminées et ceci augmentera la chaleur qui dégradera la pile. Comme avec les fissures, les zones de délaminage sont susceptibles d’être écaillées. 4. Les groupements d’électrolyte, dus au mélange du catalyseur, créent des régions à grande résistance électrique. 5. Les groupements de catalyseur peuvent être formés suite à un mélange insuffisant du catalyseur, ceci réduit l’activité catalytique. 6. La variation de l’épaisseur de la couche du catalyseur et/ou de la membrane. Deux possibilités existent : Un électrolyte mince couplé à un catalyseur de grande épaisseur présentera une dégradation due à la chaleur des réactions qui sont favorisées par l’épaisseur du catalyseur. Les régions où l’électrolyte épais est couplé à un catalyseur mince présentent une résistance ionique importante. Ainsi, le contrôle du processus de fabrication des membranes recouvertes de catalyseur permet d’améliorer la qualité d’une pile. Dans l’article [Saisset 2006], l’auteur propose un modèle Bond Graph de la PEM, basé sur une approche énergétique car la pile est un système multidisciplinaire. Le formalisme Bond Graph est un outil graphique pour la description des échanges d’énergie à l’intérieur d’un système. De plus, l’approche énergétique consiste à définir et à modéliser les échanges d’énergie à l’intérieur du système. La résistance de la membrane et la capacité de la double couche sont obtenues par la spectroscopie d’impédance. Plusieurs hypothèses ont été faites (voir §2.1). Les potentiels à l’anode et à la cathode sont fonction de la variation de l’énergie libre de Gibbs qui dépend de la pression d’oxygène et d’hydrogène. A l’anode par exemple : E anode = ΔG anode nF (1.74) 41 et ΔG anode = ΔG 0 anode - RT ln(p O2 ) 2 (1.75) où n est le nombre de moles d’électrons échangés. Ainsi ces potentiels, fonctions de la pression des gaz, sont modifiés par une variation de pression due par exemple à une perte d’intégrité des joints du stack. Le couplage entre la pression des gaz et ce phénomène physique existe donc et les conséquences qui en découlent influent sur la performance de la pile. Dans l’article [Haddad 2006], l’auteur propose un modèle dynamique non linéaire pour une pile de type PEMFC, ce qui donne une approche plus réaliste du comportement dynamique de la pile. La pression, la température et le taux d’humidification sont pris en considération. Le modèle proposé est basé sur des aspects électriques et thermodynamiques. Une représentation mathématique d’espace-état est donnée. La simulation est faite sous MATLAB - Simulink. Une approche électrique donne la tension de la pile qui dépend des pressions d’oxygène et d’hydrogène : Vcell = 1.273 - 2.7645x10-4T + p H p O2 RT ln( 2 ) - ηact - ηconc - ηohm nF p H2O (1.76) La tension de sortie de la pile est directement liée à la pression des gaz à travers l’équation (1.76). Ici aussi, le couplage entre la pression et les phénomènes physiques existe : le serrage des PB, la compression inhomogène de la pile, les connexions desserrées influent sur la pression qui devient variable. La performance de la pile à travers sa tension de sortie s’en trouve modifiée. 1.4 Techniques de caractérisation des piles à combustible Cette section décrit les deux techniques majeures de caractérisation des PAC, dans le but de juger qualitativement et quantitativement les modèles. Les techniques de caractérisation se divisent en deux parties : 42 • Les techniques de caractérisation électrochimiques In Situ : utilisation de la tension, du courant et du temps pour caractériser la performance d’une PAC dans les conditions d’opération. • Les techniques de caractérisation Ex Situ : les composants sont enlevés de l’environnement de la PAC et sont caractérisés. 1.4.1 Techniques de caractérisation électrochimiques In Situ Dans les techniques de caractérisation électrochimiques In Situ, on trouve : • Les mesures courant-tension (j-V) : La technique potentiostatique dans laquelle la tension est contrôlée par l’utilisateur et le courant résultant est mesuré. La technique galvanostatique dans laquelle le courant est contrôlé par l’utilisateur et la tension résultante est mesurée. Ces deux techniques sont utilisées en statique et les conditions de test (échauffement, pression, température, taux du flux, force de compression) doivent être notées. Les trois méthodes suivantes sont réalisées en dynamique. • La spectroscopie d’impédance électrochimique : C’est la technique la plus utilisée pour déterminer les différentes pertes. L’impédance du système est mesurée et la courbe de Nyquist résultante est adaptée à un circuit équivalent à la PAC. Ainsi les pertes ohmiques, d’activation et de transport de masse sont résolues séparément. Pour la conduction ohmique, le circuit équivalent est simple : Z=R (1.77) Pour les réactions électrochimiques et à l’interface, le circuit équivalent est une résistance en parallèle avec un condensateur : Z= 1 1 Rf + jwCdl (1.78) 43 Pour le transport de masse, le modèle équivalent est l’élément de Warburg dont l’impédance est (pour une couche de diffusion d’épaisseur infinie) : Z= • σi (1 - j) w σi (1.79) le coefficient de Warburg pour l’espèce i. Le circuit équivalent et la courbe de Nyquist de l’impédance infinie de Warburg est donné à la figure 1.3. Fig. 1.3 - Circuit équivalent et courbe de Nyquist de l’élément de Warburg [O’Hayre 2006 - 3]. Dans le cas des PAC, nous utilisons le modèle de Warburg poreux, appelé aussi élément de diffusion O, qui a la forme suivante : Z= ⎛ jw ⎞ σi (1 - j) tanh ⎜⎜ δ ⎟⎟ w ⎝ Di ⎠ • δ l’épaisseur de la couche de diffusion • Di le coefficient de diffusion de l’espèce i. (1.80) L’impédance Z est fonction de l’épaisseur de la GDL et du coefficient de diffusion qui sont des paramètres physiques dépendant de couplages multiphysiques. Ainsi, les 44 pertes dues au transport de masse, représentées par l’élément de Warburg poreux d’impédance Z sont couplées aux sollicitations mécaniques. Le circuit équivalent et la courbe de Nyquist de l’impédance du modèle de Warburg poreux est donné à la figure 1.4. Fig. 1.4 - Circuit équivalent et courbe de Nyquist de l’élément de Warburg poreux [O’Hayre 2006 - 4]. Un circuit équivalent simple du modèle d’une PAC qui « souffre » de pertes d’activation à l’anode et à la cathode, de pertes ohmiques et de pertes de transport de masse à la cathode est donné figure 1.5. 45 • Rf représente la résistance cinétique du processus électrochimique à chaque électrode • Cdl représente la séparation de charges entre ions et électrons à travers l’interface • RΩ représente la conduction ohmique • ZW représente le transport de masse. Fig. 1.5 - Circuit équivalent simple et courbe de Nyquist d’une PAC [O’Hayre 2006 - 5]. • Les mesures par interruption de courant : Cette technique, attrayante pour les systèmes de haute puissance, distingue les processus ohmiques des processus non ohmiques d’une PAC. La montée immédiate de la tension après une interruption abrupte du courant est associée aux processus ohmiques alors que la montée de la tension qui dépend du temps est associée aux processus d’activation et de transport de masse. Voir la figure 1.6. 46 Fig. 1. 6 – Technique de mesure par interruption de courant [O’Hayre 2006 – 6]. (a) Circuit équivalent simplifié d’une PAC (b) Profil de l’interruption de courant appliqué à cette pile (c) Temps de réponse de la pile • La voltamétrie cyclique (CV) : Dans cette technique, le potentiel du système est balayé entre deux tensions limites et le courant est alors mesuré. Voir figure 1.7. En général, cette technique est utilisée pour déterminer l’activité du catalyseur in situ. 47 Fig. 1.7 – Technique de voltamétrie cyclique [O’Hayre 2006 – 7]. (a) Balayage de la tension entre deux valeurs limites (b) Courant résultant Dans l’article [Liu 2006], la durabilité de la PEMFC a été étudiée sous un courant cyclique et sous un courant constant. La cause de dégradation la plus dominante dans le cas d’un courant cyclique est le taux de mélange d’hydrogène (hydrogen crossover) dû à la formation de trous dans la membrane alors que les limitations dues au transport de masse sont la cause majeure de dégradation de performance dans le cas d’un courant constant. Le diagnostic comprend la polarisation de la pile avec des courants ascendants et descendants, la mesure d’impédance, la courbe de Tafel (tension en fonction du logarithme de la densité de courant), le taux de mélange d’hydrogène et la surface électrochimique active. Les courants cycliques simulent les conditions réelles de conduite sur la route. La formation de trous dans la membrane augmente sa résistance et crée des poches emprisonnant l’eau d’où une perte de la conductivité et inondation de la pile. Ce type de couplage dégrade la performance de la pile. Dans l’article [Kim 2005] et afin d’étudier la réponse du courant à de rapides changements de tension dans une PEMFC, des expériences ont été effectuées à des flux constants et pour deux types de parcours de flux : mono et triple. Dans le parcours de flux triple, il existe une descente après la montée de la densité de courant alors que cette descente n’existe pas dans le cas d’un parcours simple. Ceci provient du fait que le parcours de flux triple permet l’exposition d’une plus grande surface à la sortie. 48 Dans l’article [Friede 2004], l’auteur présente une méthode pour la mesure de la résistance de la membrane In-Situ. Cette résistance est ensuite tracée en mode d’opération transitoire durant des changements de paramètres et des variations de courant. Une approche 1D est faite. Dans les canaux, l’équation décrivant le mélange gazeux est : ⎞ dpi RT ⎛ p = . ⎜ J i,inlet - J i,consomme - i .(J inlet - J consomme ) ⎟ dt lch ⎝ p ⎠ • pi la pression partielle de chaque gaz • Ji le flux molaire de chaque gaz • lch la profondeur du canal. (1.81) Dans les GDL, la collision entre les molécules de différents gaz est décrite par l’équation de Maxwell-Stefan. : ⎛ ⎞ ⎜ pi J j - p jJ i ⎟ dpi = -RT. ⎜ ∑ ⎟ ε dt ij ⎜⎜ p Dij ⎟⎟ τ ⎝ ⎠ • τ la tortuosité du matériau • ε la porosité de la GDL • Dij le coefficient de diffusion binaire. (1.82) Dans les électrodes, l’équation de Tafel est utilisée pour lier les chutes de tension à la densité de courant i : p Oref2 RT ⎛ i ηi = ln ⎜ . 2αF ⎜⎝ f.i 0 (T) p O2 • f ⎞ ⎟⎟ ⎠ (1.83) est un facteur donnant le rapport entre la surface active et géométrique de la cellule 49 • i 0 (T) est fonction de l’énergie d’activation. Dans la membrane, le transport de protons et d’eau est exprimé en utilisant les gradients des potentiels électrochimiques en fonction de la conductivité de la membrane qui est un coefficient primaire, du coefficient de drag et celui de diffusion et de la température. Ces potentiels électrochimiques sont reliés au contenu d’eau de la membrane. Les chutes de tension aux électrodes varient avec les pressions partielles or la pression partielle dépend de la porosité et de la diffusion qui sont couplées à des phénomènes physiques et mécaniques comme les fuites de gaz, la compression inhomogène de la pile…. Par suite, toute modification des coefficients primaires de porosité et de diffusion par des phénomènes mécaniques va faire varier la tension de sortie de la pile. 1.4.2 Techniques de caractérisation électrochimiques Ex Situ Dans les techniques de caractérisation Ex Situ, on trouve : • L’analyse de la porosité : Pour être efficaces, les électrodes et le catalyseur doivent présenter une porosité. • La détermination de la surface : Les catalyseurs les plus efficaces présentent une grande surface réelle. • La mesure de la perméabilité, c.à.d. la facilité avec laquelle un gaz se déplace à travers un matériau. • L’inspection microscopique • Les analyses chimiques. 1.5 Phénomènes mécaniques influant la production de puissance d'une PEM 1.5.1 Composants de la pile dans les publications Les études présentées ci-dessous sont celles de P. Zhou et al., Y. Zhou et al., Zhang et al., Nitta et al., Fekrazad et Bergman, Hottinen et al., Sun et al. et Al-Baghdadi et al. • P. Zhou et al. examinent [P. Zhou 2006 et P. Zhou 2007] l’effet de la force de serrage sur la résistance de contact et la porosité globale de la GDL dans une 50 PEMFC. Ils proposent une forme optimale de dent pour la plaque bipolaire (section rectangulaire avec un rayon de courbure adéquat). La pression du gaz à l’entrée et à l’interface entre la GDL et la PB est gardée constante. Par la suite, ils étudient l’effet de la force de serrage sur la performance d’une pile PEMFC, la déformation de la membrane et du catalyseur n’étant pas prise en compte. La résistivité électrique de contact entre la GDL et la PB est donnée par : A⎛ B⎞ R= ⎜ ⎟ L⎝T⎠ C (1.84) A, B et C sont des paramètres constants pour chaque type de PB et de GDL. L est la longueur de la région de contact et T est la contrainte de traction qui n’est autre qu’une force divisée par une surface. R est un coefficient global de résistivité qui change avec la longueur du contact et la contrainte de traction. La porosité de la GDL après compression est ε , elle est donnée par : ε0 - ev 1- e v ε= (1.85) La perméabilité K dépend de la porosité : K = K0 • ε3 s1.5 2 (1- ε) s (1.86) représente la saturation de l’eau liquide. La compression de la pile modifie la porosité et la perméabilité de la GDL ainsi que la résistance de contact entre la GDL et la PB. Un couplage existe donc et les résultats numériques montrent qu’une plus grande force de serrage fait diminuer la porosité de la GDL sous la dent ce qui fait augmenter la résistance au transport des deux phases liquide et vapeur. Cependant, une plus grande force de serrage fait diminuer la résistance de contact et minimise les pertes dues à la résistance électrique dans la pile. 51 Il existe donc un optimum pour la force de serrage et qui fournit un maximum de densité de puissance. • Dans une autre étude, Y. Zhou et al. [Y. Zhou 2009] étudient l’effet de la combinaison des propriétés mécaniques et électrochimiques d’une PEMFC sur sa performance. La déformation de la GDL, le transfert de masse, la résistance électrique de contact ainsi que la porosité de la GDL sont modélisés suite à différentes pressions d’assemblage. Toutefois, une porosité moyenne globale est prise en compte, et non pas un champ de porosité variable localement. De plus, la perméabilité est maintenue constante. Un modèle 3D est développé. Le flux de gaz est obtenu par la résolution des équations de Navier-Stokes en régime continu : ∂u η − ∇.[η (∇u + (∇u)T )] + u + ∇p = 0 ∂t K ∇.u = 0 ρ • ρ la masse volumique • u la vitesse du flux • η la viscosité dynamique • K la perméabilité • p la pression. (1.87) Les flux de masse en phase gazeuse sont basés sur l’équation de diffusion/convection de Maxwell-Stefan : N ∇.[- ρωi ∑ Dij { j =1 ∇M M (∇ω j + ω j )} + ωi ρu ] = 0 Mj M • ωi la fraction massique de l’espèce i • D ij le coefficient de diffusion binaire (1.88) 52 • Mj la masse moléculaire de l’espèce j • M la masse molaire totale du mélange. Les diffusivités binaires sont telles que : 1.5 Dij = Dij ,0 p ⎛T ⎞ ⎜ ⎟ p0 ⎝ T0 ⎠ (1.89) Dans la GDL, la loi de Darcy modélise le flux : u=- K η ∇p (1.90) A la surface du catalyseur, la densité de courant est : ic = i0,refc ωO Fηact ,c exp(-α ) RT ωO ,0 2 (1.91) 2 et le flux massique d’oxygène consommé est : N O2 = - ic M O2 4F (1.92) La tension de sortie de la pile est : Vcell = 0.2329 + 0.0025T − η (1.93) Ici, le couplage avec les phénomènes mécaniques provient de la pression d’assemblage de la pile. En effet, la vitesse du flux dépend de la perméabilité K qui dépend de la porosité. Ainsi, toute pression sur la pile va modifier la porosité de la GDL, donc sa perméabilité et par suite la vitesse du flux. De plus, la compression de la pile influe sur le volume (donc sur la masse volumique) et sur les coefficients de diffusion. Les résultats ont montré qu’une plus grande pression d’assemblage augmente la résistance de la cellule de pile au transfert de masse et cause une distribution non homogène de la densité de courant. Cette pression réduit aussi la résistance de contact. Les surtensions sont alors modifiées ainsi que la tension de sortie de la pile. 53 • D’autre part, Zhang et al. [Zhang 2006] proposent deux méthodes (une expérimentale et une numérique) pour l’estimation de la résistance de contact entre les plaques bipolaires et les couches de diffusion de gaz. La relation constitutive de la résistance de contact est : Rcontact = 2.2163 + 3.5306 pcontact (1.94) Il est démontré dans cette étude que la résistance de contact est couplée avec la valeur moyenne de la pression de serrage durant l’assemblage de la pile et varie en conséquent. • Nitta et al. [Nitta 2007] étudient la compression non homogène de la couche de diffusion de gaz d’une PEMFC due à la structure canal/dent de la plaque bipolaire. Dans la GDL, la résistance du flux est caractérisée par la perméabilité du gaz : u=- K η ∇p (1.95) La compression affecte la résistance qui peut être exprimée par : RPB ,GDL ( z ) = z (1.96) σ GDL . A • σ GDL est la conductivité de la GDL • A la surface de la GDL. Dans cette étude, la pression est gardée constante et les résultats montrent qu’une plus grande compression de la GDL améliore sa conductivité, réduit la résistance de contact entre la GDL et l’électrode et influence le transport de masse (réactifs et eau) par modification de la porosité qui est considérée comme homogène au départ. Il existe donc un couplage de la compression non homogène de la GDL par la structure canal/dent de la PB avec la perméabilité (qui est fonction de la porosité) et la conductivité de la GDL. La vitesse du flux et la résistance de contact dépendent de 54 ces paramètres primaires et varieront en conséquent. Ceci devra être pris en considération dans la modélisation des PEMFC. • Fekrazad et Bergman [Fekrazad 2007] développent un modèle bidimensionnel de la PEMFC afin d’étudier l’effet des variations de la force de compression durant l’assemblage de la pile. Ces variations créent aussi des modifications dans les propriétés de la pile. La pression d’interface entre la GDL et la PB et la porosité de la GDL sont gardées constantes. La densité de puissance de la pile et les différences de température dans la membrane sont obtenues en fonction de la pression de serrage. Les densités de courants de transfert j aux électrodes sont calculées par l’équation de Butler-Volmer : ja = aj • ref ⎛ xH 2 ⎜⎜ ref ⎝ xH 2 aj ref 0.5 ⎞ ⎧ F F ⎛ ⎞ ⎛ ⎞⎫ ηa ⎟ − exp ⎜ −α c ηa ⎟ ⎬ ⎟⎟ ⎨exp ⎜ α a RT ⎠ ⎭ ⎝ RT ⎠ ⎝ ⎠ ⎩ (1.97) la densité de courant d’échange multipliée par la surface de l’électrode • xi la fraction molaire de l’élément i • αi le coefficient de transfert à chaque électrode. L’équation de Nernst donne la tension du circuit : Vcell ⎛ pH ⎜ 02 RT ⎜ p ln = 1.482 − 0.000845T + 2F ⎜ ⎜ ⎜ ⎝ • p0 ⎛ pO ⎞ . ⎜ 02 ⎟ ⎝ p ⎠ pH 2O p0 0.5 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (1.98) la pression de référence, égale à 101.3 kPa. Les densités de courant dépendent de la température et la tension de sortie de la PAC dépend de la température et de la pression des gaz. 55 Il existe aussi un couplage entre la variation de la pression de serrage p de la pile et les résistances de contact et la perméabilité de la couche de diffusion de gaz. R electrical ,cathode = ( −3.65*103 p 2 + 1.3*105 p + 1.3*105 ) R thermal ,cathode = ( −4.44 p + 144.2 p + 387 ) 2 −1 −1 K = (1.3*10−14 p 2 − 3.1*10−12 p + 2.15*10−11 ) (1.99) (1.100) Les résultats obtenus ont montré que l’augmentation de la pression de serrage réduit la résistance de contact et la perméabilité de la couche de diffusion de gaz. Ce type de couplage augmente la performance de la pile et change la distribution de température interne. • Hottinen et al. [Hottinen 2007] étudient les effets de la compression de la GDL sur les transferts de charge et de masse dans la pile. Le modèle utilise les propriétés physiques (la porosité, la perméabilité, la conductivité et la résistance de contact entre la GDL et l’électrode) obtenues expérimentalement en fonction de l’épaisseur de la GDL. Dans la GDL côté cathode, les équations de conservation de masse, de mouvement, d’espèces et de charges sont résolues : r ∇.( ρ u ) = 0 η r ∇p = − u K r ∇.N i = 0 ∇.(−σ GDL , x (1.101) ∂φGDL ,c r ∂φ r ex − σ GDL , y GDL ,c ey ) = 0 ∂x ∂y • φ le potentiel électronique • ρ la masse volumique • K la perméabilité • σ GDL , x la conductivité dans le plan de la GDL 56 • σ GDL , y la conductivité à travers le plan de la GDL • Ni flux molaire de l’espèce i. Dans la cathode/électrode, il ya résolution des équations de conservation de masse, d’espèces et de charges mais avec l’ajout de termes source car de l’oxygène est consommé et de l’eau est produite ainsi que le courant protonique est changé en courant électronique : jc M O2 jc M H 2O r ∇.( ρu ) = − + 4F 2F r j ∇.N O2 = − c 4F r j ∇.N H 2O = c 2F ∇.(−σ e∇φe,c ) = jc (1.102) ∇.(−σ m∇φm ) = − jc L’équation de Butler Volmer donne la production de courant jc : jc = j0,c cO2 ref cO2 exp( −α F η) RT • η est le potentiel de surtension • α est le facteur symétrique de réaction. (1.103) Dans la membrane, l’équation de conservation de charges est résolue : ∇ 2 .φm = 0 (1.104) • φm est le potentiel protonique. Dans le côté anode / GDL, l’équation de conservation de charges est résolue. ∇.(−σ GDL , x ∂φGDL , a r ∂φ r ex − σ GDL , y GDL ,a ey ) = 0 ∂x ∂y (1.105) A l’anode, le courant électronique est consommé et le courant protonique est produit : 57 ∇.(−σ e∇φanode ) = − ja (1.106) ∇.(−σ m∇φm ) = ja tel que le courant ja produit à l’anode est : ja = j0, a zF RT (φe , a − φm ) • z (1.107) est le nombre d’électrons transférés. Les équations de la variation de l’épaisseur h de la GDL, de la porosité ε, de la perméabilité K, de la conductivité σ et de la résistance de contact rcont entre la GDL et la dent des PB sont obtenues à partir des données expérimentales par interpolation. h( x) = 19.30314 log(( x − 0.0005) x106 + 1) x10−6 + hcomprimee ε ( x) = ε 0 h( x) − hmin h0 − hmin K ( x) = −1.700 x10−11 + 2.760 x10−7 h( x) − 1.484 x10−3 h( x) + 2.754h( x) 2 3 σ GDL , x ( x) = 6896 − 1.159 x107 h( x) σ GDL , y ( x) = 3285 − 8.385 x106 h( x) rcont = 5.83 x10−10 exp(2.06 x104 h( x)) (1.108) Ces grandeurs sont fonction de l’épaisseur de la GDL qui varie en fonction de la force appliquée sur la pile. Ainsi le phénomène de compression inhomogène de la GDL est étroitement couplé à la distribution de courant à cause de la variation des propriétés de la GDL (épaisseur, porosité, perméabilité, conductivité) et de la résistance de contact entre la GDL et l’électrode. • Sun et al. [Sun 2005] proposent un modèle bidimensionnel, en régime stationnaire et isotherme pour étudier l’influence des propriétés de la couche de diffusion de gaz et la géométrie des canaux de flux sur le taux de réaction local dans la couche cathode/catalyseur d’une PEMFC mais la pression à l’interface GDL/PB est gardée constante. 58 Le transport des espèces chimiques dans la GDL et le catalyseur est décrit par les équations de Maxwell-Stefan. N ⎡ M ∇M ⎤ (∇ω j + ω j )⎥ ∇.ni =∇. ⎢- ρωi ∑ Dijeff Mj M ⎥⎦ j =1 ⎢⎣ Dijeff = Dij ε GDL 1.5 (1.109) Dijeff,cathode = Dij ε cathode 1.5 Où ni est le flux massique des espèces A la cathode, la loi d’ohm décrit la migration des protons et le transport des électrons : i p = − k1eff ∇φ1 (1.110) ie ,C = − kSeff ∇φS • ip flux protonique • ie,C flux électronique dans le catalyseur • k1 conductivité protonique • kS conductivité électronique • φ1 potentiel de l’électrolyte • φS potentiel électrique. Ici aussi, la distribution des courants est influencée par les différentes conductivités et potentiels. Ainsi, les variables soulignées dans les équations du modèle proposé sont influencées par le couplage avec les phénomènes mécaniques comme le serrage des vis durant l’assemblage de la pile. La porosité, le volume (masse volumique) et les conductivités sont influencés par le phénomène de compression et par suite, les coefficients de diffusivité et les différents flux. 59 D’autre part, les résultats ont montré que la compression de la pile fait diminuer la porosité et augmenter la conductivité de la GDL sous la dent et une augmentation du rapport géométrique canal/dent améliore le transport d’eau et le taux de la réaction de réduction d’oxygène. • Al-Baghdadi et al. [Al-Baghdadi 2007 et Al-Baghdadi 2009] développent un modèle complet, tridimensionnel, non isotherme, en régime stationnaire et à une seule phase d’une pile PEMFC. Le modèle prend en compte les phénomènes majeurs de transport dans la pile : chaleur de convection et de diffusion, transport de masse, cinétique des électrodes et champs de potentiels. La porosité et les pressions à l’anode et à la cathode de la GDL sont considérées comme homogènes. Un calcul précis des surtensions d’activation est fait dans cet article. La température, la pression, le coefficient stœchiométrique, la largeur des canaux de flux des gaz, l’épaisseur, la porosité et la conductivité thermique de la couche de diffusion des gaz, l’épaisseur de la membrane et les potentiels sont étudiés afin de voir leurs effets sur la distribution de la densité de courant. L’ensemble des équations (1.111) suivantes est résolu dans le modèle. Les propriétés physiques de variation de l’épaisseur de la GDL, de variation des conductivités de la membrane et des électrodes, dus à des sollicitations mécaniques, sont couplées avec la tension de sortie de la pile à travers leur influence sur les densités de courant et les différentes surtensions. Vcell = E thermo - ηact - ηohm' - ηmem - ηconc 1 E thermo = 1.229 - 0.83x10-3 (T-298.15) + 4.3085x10-5T[ln p H 2 + ln p O2 ] 2 ⎛ C O2 ref i c = i 0,c ⎜⎜ ref ⎝ C O2 ⎞⎡ ⎛α F ⎞ ⎛ αF ⎞⎤ ⎟⎟ ⎢ exp ⎜ a ηact,c ⎟ + exp ⎜ - c ηact,c ⎟ ⎥ ⎝ RT ⎠ ⎝ RT ⎠⎦ ⎠⎣ ⎛ CH2 ref i a = i 0,a ⎜⎜ C ref ⎝ H2 ⎞ 2⎡ ⎛α F ⎞ ⎛ αF ⎞⎤ ⎟⎟ ⎢ exp ⎜ a ηact,a ⎟ + exp ⎜ - c ηact,a ⎟ ⎥ ⎝ RT ⎠ ⎝ RT ⎠⎦ ⎠ ⎣ 1 ∇.(k e∇ηohm' ) = 0 ∇.(k m∇ηmem ) = 0 60 ηconc,c = RT ⎛ i c ⎞ ln ⎜1⎟ 2F ⎜⎝ i L,c ⎟⎠ ηconc,a = RT ⎛ i a ln ⎜12F ⎜⎝ i L,a i L,c = i L,a = ⎞ ⎟⎟ ⎠ 2FD O2 cO2 δ GDL 2FD H 2 c H 2 δ GDL • ke la conductivité électronique effective de l’électrode • km la conductivité ionique de la membrane • δ GDL l’épaisseur de la GDL • D le coefficient de diffusion • c la concentration. En résumé : Des températures d’opération élevées augmentent la puissance et l’efficacité de la pile, ainsi le coût serait moindre. De plus, les gradients de température à l’intérieur de la pile seraient réduits. De même pour des pressions d’opération élevées et des flux à coefficients stœchiométriques élevés. Cependant dans ce cas, il y aura une puissance parasite pour comprimer les réactifs et le stack devra supporter une plus grande pression ce qui augmenterait le coût de la pile. Dans le cas de pressions élevées et dans le cas de coefficients stœchiométriques élevés, il doit y avoir un optimum entre la performance de la pile et le prix du souffleur. Un canal de flux des gaz plus étroit augmente la puissance et l’efficacité de la pile; le coût serait moindre et le gradient maximal à l’intérieur de la pile diminue mais la pression augmenterait avec un canal étroit. Pour ce qui est de la couche de diffusion des gaz, une couche mince augmenterait faiblement la puissance et l’efficacité de la pile et diminuerait le prix mais le gradient maximal à l’intérieur de la pile augmenterait. Pour des porosités et des conductivités 61 thermiques élevées de la couche de diffusion, le gradient maximal de température à l’intérieur de la pile serait moindre. Une membrane mince augmenterait la puissance et l’efficacité de la pile et diminuerait le prix et le gradient maximal à l’intérieur de la pile diminuerait. Cependant la porosité est considérée comme homogène. Le même modèle est ensuite développé en un modèle multiphasique qui étudie la déformation et les déplacements des éléments à l’intérieur d’une cellule suite à des changements de température et d’humidité relative. 1.5.2 Tableau comparatif de modèles relatifs au sujet de thèse Afin de mieux comparer les différents modèles rencontrés dans l’étude bibliographique, le tableau 1 a été élaboré. Mais en tenant compte du fait que le sujet de la thèse est la modélisation de la dégradation de puissance d’une PEMFC suite aux sollicitations mécaniques appliquées, le tableau résume uniquement les modèles qui ont rapport à ce sujet. 62 Modèle 1 [P. Zhou 2006] Entrées E Propriétés de la GDL, géométrie de la plaque bipolaire et force de serrage Sorties S Résistance de contact GDL/PB et porosité Formules S(E) La résistance de contact entre la GDL et la plaque bipolaire est considérée comme plusieurs résistances en parallèle. Résultats typiques La géométrie de la plaque bipolaire, la forme de la dent, les propriétés physiques de la GDL et de la plaque bipolaire et la force de serrage ont tous une influence sur la résistance électrique et la porosité. Propriétés spécifiques Résolution des équations par un logiciel d’éléments finis. Modèle 2 [P. Zhou 2007] Entrées E Pression de compression Sorties S Potentiel de la pile, densité de courant et densité de puissance Formules S(E) 1. Résolution des équations de continuité, de conservation du mouvement et d’espèces. 2. Résolution des équations du taux de transfert de masse de l’eau, des flux massiques d’oxygène et de l’eau (liquide et vapeur). 3. Résolution de l’équation de Butler-Volmer simplifié qui décrit la densité de courant local le long du catalyseur. 4. Résolution de la pression capillaire entre la phase liquide et gazeuse par la théorie des flux non saturés. 5. Calcul du potentiel de la pile (la surtension d’activation est donnée par l’équation de Tafel). 63 Résultats typiques Une plus grande force de serrage diminue la porosité de la GDL ce qui augmente la résistance au transport des deux phases liquide et vapeur. Une plus grande force de serrage diminue la résistance de contact et minimise les pertes dues à la résistance électrique dans la pile. Propriétés spécifiques Modèle à 2 phases. Le flux est stationnaire, isotherme et laminaire dans la GDL. L’eau générée est liquide. Pile PEMFC avec des distributeurs de gaz interdigitaux. Modèle 3 [Y. Zhou 2009] Entrées E Déformation de la GDL, transfert de masse et résistance électrique de contact dus à différentes pressions de serrage Sorties S Tension et densité de courant de la pile. Formules S(E) 1. Dans les canaux de flux des gaz : Résolution des équations de Navier-Stokes Le flux de masse est calculé par l’équation de diffusion et de convection de Maxwell-Stephan 2. Dans la GDL poreuse: Résolution de l’équation de Darcy 3. Dans les couches de catalyseurs : Résolution de l’équation de la distribution de la densité de courant Résolution des équations électrochimiques de flux Résultats typiques Une plus grande pression d’assemblage augmente la résistance de la cellule de pile au transfert de masse et cause une distribution non homogène de la densité de courant. Cette pression réduit cependant la résistance de contact. Une pression d’assemblage optimale est observée. 64 Propriétés spécifiques Modèle multiphysique 3D qui prend en compte la combinaison des propriétés mécaniques et électrochimiques d’une PEMFC. Modèle 4 [Zhang 2006] Entrées E Relation constitutive entre la résistance de contact et la pression Sorties S Résistance de contact entre la plaque bipolaire et la GDL Formules S(E) Relation géométrique expérimentale entre la résistance de contact et la pression Résultats typiques La résistance de contact est influencée par la valeur moyenne de la pression de serrage durant l’assemblage de la pile. Propriétés spécifiques Deux méthodes semi-empiriques pour l’estimation de la résistance de contact entre plaque bipolaire et GDL. Modèle 5 [Nitta 2007] Entrées E Compression non homogène de la GDL due à la structure canal/dent de la plaque bipolaire Sorties S Perméabilité, conductivité, résistance de contact Formules S(E) Les équations de la perméabilité, la conductivité et la résistance de contact sont fonction de l’épaisseur de la GDL Résultats Une plus grande compression de la GDL améliore sa conductivité, réduit la résistance de contact entre la GDL et l’électrode et 65 typiques entrave le transport de masse (réactifs et eau). Propriétés spécifiques Propriétés de la GDL évaluées expérimentalement. Modèle 6 [Fekrazad 2007] Entrées E Pression de serrage Sorties S Densité de puissance de la pile, température dans l’AME de la pile Formules S(E) 1. Dans les plaques bipolaires : Résolution des équations de l’énergie thermique et du transport d’électrons. 2. Dans les couches de diffusion des gaz : Résolution des équations de conservation de masse, de mouvement, d’énergie thermique, d’espèces et de charges électroniques. 3. Dans la membrane : Résolution des équations de conservation de masse, d’énergie thermique, d’espèces et de charges ioniques. 4. Dans le catalyseur : Résolution des équations de conservation de masse, de mouvement, d’énergie thermique, d’espèces et de charges. Un terme source apparaît dans les équations de conservation de masse, d’énergie thermique, d’espèces et de charges. Résultats typiques L’augmentation de la pression de serrage réduit la résistance de contact et la perméabilité de la couche de diffusion de gaz ce qui augmente la performance de la pile et change la distribution de température interne. Propriétés spécifiques Modèle statique, bidimensionnel. Les équations des résistances de contact électriques et thermiques ainsi que la perméabilité sont obtenues à partir des données expérimentales par interpolation. 66 Résolution des équations par un logiciel d’éléments finis et résolveur FEMLAB. Modèle 7 [Hottinen 2007] Entrées E Compression inhomogène de la GDL Sorties S Courbe de polarisation de la pile, fractions molaires de l’oxygène, distribution de la densité de courant Formules S(E) 1. Côté cathode GDL : Résolution des équations de conservation de masse, de mouvement, d’espèces et de charges. 2. Côté cathode électrode : Résolution des équations de conservation de masse, de mouvement, d’espèces et de charges mais avec l’ajout de termes source car de l’oxygène est consommé et de l’eau est produite ainsi que le courant protonique est changé en courant électronique. 3. Dans la membrane : Résolution de l’équation de conservation de charges. 4. Côté anode GDL : Résolution de l’équation de conservation de charges. Résultats typiques L’effet de la compression non homogène de la GDL sur les transferts de charge et de masse est étudié. Ceci influe sur la distribution de la densité de courant. Aussi, des gradients de température apparaissent et diminuent la durée de vie et la performance de la pile. Le transfert de masse côté anode est négligé. Propriétés spécifiques Modèle bidimensionnel, isotherme avec propriétés homogènes de la GDL et modèle bidimensionnel avec prise en compte de la compression non homogène de la GDL. Les équations de la variation de l’épaisseur de la GDL, de la porosité, de la perméabilité, de la conductivité et de la résistance de contact sont obtenues à partir des données expérimentales par interpolation. Résolution des équations par un logiciel d’éléments finis COMSOL Multiphysics version 3.2b. 67 Modèle 8 [Sun 2005] Entrées E Propriétés de la GDL et géométrie des plaques Sorties S Densité de courant et taux de la réaction d’oxygène Formules S(E) 1. Dans la GDL : Résolution des équations de Maxwell-Stefan de transport des espèces chimiques et de conservation de charges. 2. Dans la cathode : Résolution des équations du taux de la réaction d’oxygène et de transport de protons et d’électrons. 3. Dans le catalyseur : Résolution des équations de Maxwell-Stefan de transport des espèces chimiques, de conservation d’espèces et de charges. Résultats typiques Une augmentation du rapport canal-dent améliore le transport d’eau et le taux de la réaction de réduction d’oxygène. La compression de la GDL suite à l’assemblage de la pile fait diminuer la porosité et augmenter la conductivité sous la dent. Une augmentation de la conductivité de la GDL augmente la densité de courant produite dans la région sous le canal alors qu’une modeste augmentation est produite sous la dent. Propriétés spécifiques Modèle bidimensionnel, en régime stationnaire et isotherme. Résolution des équations par un logiciel d’éléments finis et le solveur d’équations à dérivées partielles FEMLAB. Modèle 9 [Al-Baghdadi 2007] Entrées E Densité de courant nominal désirée 68 Sorties S Formules S(E) Surtensions d’activation et donc potentiel de la pile 1. Dans les canaux de flux des gaz : Résolution de l’équation de continuité de Navier-Stokes en régime stationnaire et de l’équation du mouvement. Résolution de l’équation de transport de masse en statique (diffusion et convection). Résolution de l’équation de l’énergie convective (champ de température). 2. Dans les couches de diffusion des gaz : Résolution de l’équation de continuité et de l’équation du mouvement. Résolution de l’équation de transport de masse en milieu poreux. Résolution de l’équation de transfert de chaleur. Résolution de l’équation de distribution du potentiel. 3. Dans le catalyseur : Résolution de l’entropie spécifique de l’oxygène, de l’hydrogène et de production d’eau. Résolution de l’équation de génération de chaleur. Résolution des équations de Butler-Volmer de distribution de la densité de courant local. 4. Dans la membrane : Résolution de l’équation du flux net d’eau à travers la membrane (équilibre entre le mouvement d’eau de l’anode vers la cathode et de la cathode vers l’anode). Résolution de l’équation de transfert de chaleur. Résolution de l’équation de perte de potentiel due à la résistance de la membrane au transport protonique. 5. Potentiel de la cellule : Résolution de l’équation du potentiel. Les surtensions d’activation à la cathode et à l’anode sont calculées. 69 Les surtensions de diffusion à la cathode et à l’anode sont calculées. 6. Conditions aux limites : Les valeurs d’entrée à la cathode et à l’anode sont prévues pour la vitesse, la température et les concentrations des espèces. Les vitesses d’entrée de l’air et du carburant sont calculées à partir de la valeur désirée de la densité de courant. A la sortie, la pression seulement est prévue comme étant la pression désirée de l’électrode. Résultats typiques Des températures d’opération élevées augmentent la puissance et l’efficacité de la pile, ainsi le coût est moindre. De plus, les gradients de température à l’intérieur de la pile sont réduits. De même pour des pressions d’opération élevées et des flux à coefficients stœchiométriques élevés. Cependant dans ce cas, il y a une puissance parasite pour comprimer les réactifs et le stack devra supporter une plus grande pression ce qui augmente le coût de la pile. Dans le cas de pressions élevées et dans le cas de coefficients stœchiométriques élevés, il doit y avoir un optimum entre la performance de la pile et le prix du souffleur. Un canal de flux des gaz plus étroit augmente la puissance et l’efficacité de la pile ; le coût est moindre et le gradient maximal à l’intérieur de la pile diminue mais la pression augmente avec un canal étroit. Pour ce qui est de la couche de diffusion des gaz, une couche mince augmente faiblement la puissance et l’efficacité de la pile et diminue le prix mais le gradient maximal à l’intérieur de la pile augmente. Pour des porosités et des conductivités thermiques élevées de la couche de diffusion, le gradient maximal de température à l’intérieur de la pile est moindre. Une membrane mince augmente la puissance et l’efficacité de la pile et diminue le prix et le gradient maximal de température à l’intérieur de la pile diminue. Propriétés spécifiques Modèle tridimensionnel, non isotherme, statique et à une seule phase. Discrétisation des équations par volumes finis et résolution par un code CFD (Computational Fluid Dynamics) Tableau 1 - Tableau récapitulatif des différents modèles étudiés 70 1.5.3 Le système pile Un ensemble de plusieurs cellules est nécessaire pour obtenir la tension désirée puisqu’une cellule unique ne fournit que 0,6 - 0,7 V durant son opération. L’assemblage de plusieurs cellules ainsi que les sous-systèmes peuvent provoquer des couplages mécaniques/multiphysiques influençant les performances de la pile entière. C’est pourquoi la description des sous-systèmes a été incluse dans l’étude bibliographique. ƒ Le sous-système PAC Ce sous-système convertit un gaz riche en hydrogène et un oxydant en électricité et chaleur. Connectées en série pour répondre aux exigences de tension des applications réelles, les cellules sont empilées. Cette empilement ou stack doit être simple et de fabrication non onéreuse, les connexions électriques entre les cellules à faibles pertes, le collecteur des réactifs gazeux efficace, le système de refroidissent efficace et le scellage entre les différentes cellules fiable. La forme la plus répandue de l’interconnexion entre les cellules est l’empilement vertical ou à plaques bipolaires. Ainsi les plaques bipolaires connectent électriquement les piles entre elles et fournissent un chemin pour le flux des gaz. Cependant cette configuration bipolaire est difficile à sceller ; les gaz fuiront des électrodes poreuses. Pour cela, l’électrolyte est conçu un peu plus large que les électrodes et un scellage adéquat est mis autour des électrodes selon la figure 1.8. Fig. 1.8 - Méthode de scellage [O’Hayre 2006 - 8]. 71 D’autres configurations existent. Dans tous les cas, le meilleur empilement est celui qui minimise le nombre de scellages. Notons que la performance moyenne d’une cellule dans l’empilement est moindre que celle d’une cellule unique opérant dans les mêmes conditions [Baschuk 2004]. Dans l’article [Barreras 2005], la distribution de l’écoulement dans la structure de flux d’une PEMFC commerciale est étudiée des points de vue expérimental et numérique. La structure est bipolaire et a une surface de 50 cm2, elle est formée de 16 canaux centraux parallèles et 4 canaux entourent la surface d’écoulement. Cette structure, communément faite de graphite ou de certains métaux tels que le titane, guide l’écoulement des réactifs et fournit un chemin de conduction pour les électrons d’une cellule à une autre. Le centre de cette étude est la visualisation de l’écoulement de l’hydrogène. Cependant, des vapeurs telles que l’acétone nécessaires pour visualiser l’écoulement ne peuvent être utilisées avec des matériaux plastiques rapidement dégradables. La visualisation de l’écoulement dans la structure est réalisée par PLIF « Planar Laser-Induced Fluorescence » et la mesure de la vitesse est faite par la trace d’un colorant injecté dans la structure. Pour cela, une étude a été réalisée afin d’utiliser un liquide à la place du gaz et qui conserve la même valeur du nombre de Reynolds ; le liquide utilisé est un mélange de glycérine et d’eau (49.6 % de glycérine et 51.4 % d’eau). Le système optique pour la visualisation est formé d’un rayon laser, d’une lentille, de la structure de flux et d’une caméra. Pour la simulation numérique bidimensionnelle, les équations de conservation de masse et de mouvement de Navier-Stokes en régime stationnaire sont utilisées : r ∇( ρ v) = 0 rr r ur ∇( ρ v.v) = −∇p + ∇( μ∇v) + ρ g • v la vitesse de l’écoulement (m s-1) • ρ la densité (kg m-3) • μ la viscosité absolue (Pa s). (1.112) 72 La méthode des volumes finis est utilisée pour la discrétisation des équations cidessus. Les résultats ont montré une très faible efficacité de la structure de flux vis-à-vis de la distribution de l’écoulement et une utilisation non équilibrée du catalyseur. Le chemin de préférence est formé des canaux latéraux alors que l’écoulement est très lent dans les canaux centraux. Ceci peut être expliqué par la formation de bulles aux canaux d’entrée dans la zone centrale de la structure principalement. Ainsi, les résultats expérimentaux et les prédictions numériques sont conformes : La distribution de l’écoulement dans la structure est non homogène et confère une performance limitée à la pile. Aussi l’écoulement peut être influencé par les vibrations que peut subir une PAC. Ce type de couplage modifie la performance de la pile. Plusieurs travaux de recherche dans le domaine des plaques bipolaires dont le but est l’amélioration de la performance d’une PEMFC ont été publiés. Dans l’article [Yen 2006], l’auteur a développé des plaques bipolaires composites de faible masse et de performance élevée. Ces plaques sont composées d’une résine de vinyle ester, de graphite et d’argile organique (organoclays). L’argile organique a été préparée par échange ionique de montmorillonite (MMT) avec trois masses moléculaires différentes d’agents intercalaires. Cette manière avec laquelle a été préparée l’argile organique a amélioré les propriétés des plaques composites telles que la porosité, la protection anticorrosive, les forces d’impact… Les tests I-V et I-P des plaques bipolaires composites à composition optimale ont montré une performance équivalente à celle des plaques bipolaires à graphite seul. Dans un autre article [Blunk 2006 - 1], il a été démontré que des plaques composites à grande concentration de graphite pourraient satisfaire les attentes des PEMFCs du point de vue conductivité électrique, mais ont tendance à devenir fragiles et à générer des taux d’infiltration d’hydrogène quand l’épaisseur des plaques est très faible (cette faible épaisseur est requise pour atteindre les densités de puissance volumétriques requises dans le cas d’applications véhiculaires). Actuellement, les fournisseurs de plaques composites formulent leur matériel à des concentrations de graphite légèrement en-dessous du CPVC (critical pigment 73 volume concentration) nécessaire pour une grande conductivité électrique et une force adéquate des plaques. Dans l’article [Blunk 2006 - 2], l’auteur étudie la possibilité d’obtenir une grande conductivité des plaques avec une faible concentration de graphite expansé afin de satisfaire les exigences requises pour les PEMFC dans leur application transport (i.e. une résistance de surface spécifique des plaques de moins de 20 mΩ.cm2 et une largeur des plaques de moins de 2 mm pour atteindre la densité de puissance volumétrique relative aux véhicules qui doit être supérieure à 2 kW.l-1). Les plaques en métal comme le titane ou l’acier inoxydable ont de bonnes caractéristiques mécaniques et une faible infiltration des gaz. De plus, ils ont une excellente résistance à la corrosion car un film stable d’oxyde se forme à leurs surfaces ; ce film est enlevé et une couche conductrice est appliquée. Pour les plaques polymères, elles ont une bonne stabilité électrochimique dans l’environnement agressif des PEMFC mais satisfont avec difficulté aux exigences de résistance et d’épaisseur simultanément que de bonnes propriétés mécaniques. Dans [Wolf 2006], l’auteur a expérimenté de nouveaux matériaux composites pour les plaques bipolaires, de densité plus faible que celle des plaques commerciales à base de polymères et de carbone. Ces nouveaux matériaux, de conductivité électrique élevée, sont des polymères LCP (Liquid Crystal Polymer) et de fibres de carbone. La conductivité électrique et la perméabilité d’hydrogène ont été mesurées. Dans [Antoni 2007], une publication de PSA Peugeot Citroën et du Commissariat à l’Energie Atomique (CEA), les auteurs ont conçu une pile permettant de franchir une nouvelle étape en termes de puissance électrique, de rendement énergétique et de compacité. C’est la pile GENEPAC qui est constituée de quatre modules de 20 kW. Ils ont utilisé des plaques minces en acier inoxydable embouties, permettant ainsi de réduire le coût et le volume de la pile par rapport aux plaques traditionnelles en graphite ou les composites constituées de polymère chargé de graphite. Ainsi, le métal embouti diminue considérablement les épaisseurs par rapport au graphite. De plus, la fabrication des plaques nécessite des procédés facilement industrialisés. 74 ƒ Le sous-système de gestion thermique Dans une PAC, l’énergie qui n’est pas convertie en électricité est dissipée sous forme de chaleur. La pile nécessite donc un refroidissement actif. Les piles à faible température d’opération comme la pile PEMFC sont parfois refroidies de manière passive, par convection naturelle. Un refroidissement actif nécessite l’addition de canaux de refroidissement dans les plaques bipolaires. L’efficacité d’un système de refroidissent est évaluée par : efficacité = débit d'élimination de la chaleur puissance consommée par les ventilateurs... (1.113) Le refroidissement peut aussi se faire par un liquide mais dans ce cas, le liquide doit être recyclé et dans le cas de l’eau, il doit être de-ionisé. Ainsi, le refroidissement de la pile prévient son surchauffage et la formation de larges gradients de température à l’intérieur de l’empilement. De plus, la chaleur dégagée est utile pour le chauffage interne de la pile (préchauffage des gaz entrant, vaporisation de l’eau pour humidifier ces gaz entrant dans la pile) et même pour le chauffage externe (chauffage de l’habitacle de la voiture, chauffage et eau chaude de l’immeuble). Tout excès de chaleur est rejeté dans l’atmosphère. ƒ Le sous-système de traitement et de fourniture du combustible Le combustible est soit de l’hydrogène pur, soit produit à partir d’un porteur d’hydrogène. Ce sous-système fournit donc le combustible H2 mais le cas échéant, purifie aussi le gaz du monoxyde de carbone et des composés sulfureux. Les porteurs d’hydrogène typiques sont le méthane, le méthanol… Cependant, ces porteurs doivent être traités afin de produire de l’hydrogène. Le traitement peut se faire par électro-oxydation directe ou par reformage (reformage externe et reformage interne11). L’électro-oxydation directe diminue la performance de la pile à cause de complications cinétiques. Le reformage externe produit des espèces telles le CO et le CO2 qu’il faudra éliminer. 11 Le reformage interne, qui se fait à l’intérieur de la pile à la surface du catalyseur de l’anode, n’est valable que pour les piles à haute température d’opération. Le CO produit n’est pas alors un problème car il peut être utilisé directement comme carburant [O’Hayre 2006]. 75 L’électro-oxydation directe est la mieux adaptée pour les applications portables alors que le reformage est appliqué dans les applications stationnaires. Le reformage externe peut se faire de plusieurs façons : ¾ vapeur (Steam Reforming SR). C’est une réaction endothermique. 1 C x H y +xH 2O(g) ←⎯ → xCO+( y+x)H 2 ⎯⎯ → CO, CO 2 , H 2 , H 2 O 2 (1.114) Le gain en H 2 ( H 2 yield) est le plus élevé et suivant le principe de Le Chatelier12, nous pouvons augmenter ce gain et diminuer le gain en CO (CO yield) par la réaction suivante de décalage de l’eau en phase gazeuse (Water gas shift). CO+H 2 O(g) ←⎯ → CO 2 +H 2 (1.115) Le gain en H 2 qui représente le pourcentage de moles d’hydrogène dans le flot reformé à la sortie du reformeur de combustible est y H2 : y H2 = n H2 (1.116) n • n H2 le nombre de moles d’hydrogène produit par le reformeur • n le nombre total de moles de tous les gaz à la sortie. ¾ oxydation partielle. C’est une réaction exothermique. 1 1 C x H y + xO 2 ←⎯ → xCO+ yH 2 2 2 (1.117) De même que dans le reformage par vapeur, le gain en H 2 peut être augmenté et le gain en CO diminué par la réaction de décalage de l’eau en phase gazeuse (Water gas shift). ¾ Reformage auto-thermique. Ce type de reformage combine les réactions de reformage par vapeur, d’oxydation partielle et de décalage d’eau gaz dans un même processus. La réaction du reformage auto-thermique est neutre. 12 Le principe de Le Chatelier : «Si une contrainte est appliquée à un système en équilibre, le système alors se réajuste, si possible, pour réduire cette contrainte [Fuel Cell Handbook - 2]. 76 1 1 → xCO 2 +( y+z)H 2 ⎯⎯ → CO, CO 2 , H 2 , H 2 O C x H y +zH 2 O(l) + (x − z)O 2 ←⎯ 2 2 (1.118) A noter que la quantité de CO dans le flot reformé reste élevée pour les piles à faibles températures, même après la réaction de décalage d’eau en phase gazeuse. Ici, le réacteur de nettoyage de CO réduit le gain en CO à des niveaux très faibles soit par réaction chimique (méthanation sélective du CO ou oxydation sélective du CO : un catalyseur encourage une réaction qui enlève le CO sur une autre qui consommerait l’hydrogène), soit par séparation physique (par adsorption qui enlève toutes les espèces sauf l’hydrogène et ceci grâce à sa faible masse moléculaire ou par séparation par une membrane de palladium qui laisse diffuser l’hydrogène à des vitesses plus grandes que les autres espèces). Ainsi, la gestion de la chaleur dans une PAC est nécessaire pour maximiser sa restitution. Note sur le stockage de l’hydrogène : L’efficacité du stockage de l’hydrogène, qui est un point critique pour les PAC portables ou mobiles, se mesure par la densité d’énergie gravimétrique et par la densité d’énergie volumétrique. enthalpie emmagasinée du combustible masse totale du système enthalpie emmagasinée du combustible Densité d'énergie volumétrique = volume total du système Densité d'énergie gravimétrique = (1.119) Le stockage de l’hydrogène très peu dense peut se faire : ¾ sous forme de gaz comprimé à haute pression dans des cylindres spéciaux, cependant ceci pose des problèmes de sécurité, ¾ comme un liquide car l’hydrogène refroidi à 22 K devient liquide, ¾ dans un métal hybride car sous un faible chauffage, l’hydrogène est libéré. 77 ƒ Le sous-système de l’électronique de puissance Ce sous-système est formé de : ¾ la régulation. Sachant que la tension d’une pile dépend largement de la température, la pression, l’humidité et le débit des réactifs, des convertisseurs courant continu/courant continu assurent la régulation de la tension. ¾ la conversion de puissance. Afin de convertir la tension continue fournie par la pile en tension alternative monophasée ou triphasée, des convertisseurs courant continu/courant alternatif assurent la conversion. ¾ le contrôle et la surveillance. Plusieurs variables comme la température, le débit des gaz, la tension de sortie, le refroidissement et le reformage nécessitent un contrôle par des valves, pompes, capteurs… ¾ la gestion de l’alimentation. Il s’agit de faire correspondre la tension de sortie de la pile avec la demande. Pour assurer la demande électrique d’une certaine charge, un système de PAC est mis sur des unités de stockage d’énergie comme la batterie, les condensateurs ou même le grillage alternatif. La pile charge l’unité de stockage lorsque la demande d’électricité est faible ou « vend » l’excès d’électricité au réseau. ƒ Le sous-système de gestion de l’air Dans l’article [Blunier 2007], l’auteur insiste sur l’importance de la gestion de l’air dans une PEMFC. La température de sortie du système de fourniture d’air doit correspondre à celle du stack. Des modèles contrôlés de pile doivent être développés afin de prendre en compte les effets dynamiques de variation du courant. Ainsi donc, la compression de la pile et l’humidification des gaz ne doivent pas être découplées. 1.5.4 Synthèse et conclusion de l’étude bibliographique Depuis sa fabrication et son processus d’assemblage, la pile subit des sollicitations mécaniques non négligeables dues principalement au serrage des plaques bipolaires. Ces pressions appliquées compriment la pile de façon hétérogène. La GDL est comprimée sous la dent de la PB et gonfle sous le canal. 78 D’autres phénomènes physiques peuvent aussi apparaître dans la pile comme la variation de l’épaisseur de la membrane et de la GDL, l’apparition de trous dans la membrane, le délaminage entre la membrane et le catalyseur et les fissures dans le catalyseur. Aussi, les propriétés physiques de la pile comme la porosité et la perméabilité de la GDL ainsi que la résistance de contact entre la PB et la GDL sont modifiées. D’autre part, dans son environnement d’utilisation, la pile subit des vibrations mécaniques dues, par exemple, aux aspérités de la route pour les applications transport. Ces vibrations font varier la pression des gaz qui s’écoulent dans les canaux des PB. Dans les études présentées ci-haut, des valeurs globales ou moyennes sont prises pour la porosité et la perméabilité de la GDL, la pression des gaz et la résistance de contact entre la PB et la GDL. Notre apport personnel consiste à considérer des champs de porosité, perméabilité, pression des gaz, et résistance de contact variables localement afin de modéliser plus correctement l’effet de la compression inhomogène de la pile. La figure 1.9 montre la localisation spatiale des phénomènes qui seront étudiés dans les différents chapitres de ce mémoire. Il s’agit de la variation de la porosité et de la perméabilité de la GDL, de la modification des pressions des gaz à l’interface entre la GDL et la PB et de la variation de la résistance de contact entre la GDL et la PB du côté de la cathode. 79 Fig. 1.9 - Localisation spatiale des phénomènes étudiés dans ce mémoire. Le tableau 2 récapitule les différents phénomènes qui résulteraient de sollicitations mécaniques appliquées sur la pile et qui pourraient détériorer la membrane, le catalyseur, les électrodes, les couches de diffusion des gaz et affecter donc la puissance de sortie de la pile. Leurs effets sur la pile, les variables impliquées et les conséquences possibles sur la performance sont aussi explicitées. La définition des variables impliquées se trouve dans l’annexe A à la fin de ce mémoire. Variation de l’épaisseur de la membrane Phénomène physique dû à des sollicitations mécaniques Conséquences Déséquilibre dans le flux des ions H+ à travers la membrane Variables impliquées ηohm , k m Modèles intégrant [Al-Baghdadi 2007] ces variables Conséquences Une diminution de l’épaisseur de la membrane réduit sa résistance et 80 attendues sur les diminue les pertes ohmiques. performances et la durabilité. Variation de l’épaisseur de la couche de diffusion des gaz Phénomène physique dû à des sollicitations mécaniques Conséquences Influence sur le transport de masse de la GDL Variables impliquées ηconc Modèles intégrant [Al-Baghdadi 2007] et [Sun 2005] ces variables Une diminution de l’épaisseur de la couche de diffusion augmente le Conséquences attendues sur les transport de masse à travers elle et réduit ainsi les pertes de transport de performances et masse. la durabilité. Variation de la porosité de la couche de diffusion des gaz Phénomène physique dû à des sollicitations mécaniques Conséquences Influence sur le transport de masse Variables impliquées ηconc , ε Modèles intégrant [Al-Baghdadi 2007], [Sun 2005], [P. Zhou 2006] et [P. Zhou 2007] ces variables Une augmentation de la porosité de la couche de diffusion de gaz Conséquences attendues sur les améliore le transport de masse et diminue ses pertes-là. performances et la durabilité. Variation de la conductivité thermique de la couche de diffusion des gaz Phénomène physique dû à des sollicitations mécaniques Conséquences Influence sur la température de la pile Variables impliquées k mem Modèles intégrant [Al-Baghdadi 2007] 81 ces variables Une plus grande conductivité thermique améliore la puissance Conséquences attendues sur les d’opération de la pile. performances et la durabilité. Fissures dans le catalyseur Phénomène physique dû à des sollicitations mécaniques Conséquences Augmentation de la résistance du catalyseur Variables impliquées ηact , ηohm Modèles intégrant ces variables Baisse de la performance de l’activité catalytique, dégradation de la Conséquences attendues sur les performance et diminution de la tension de sortie de la pile. performances et la durabilité. Délaminage entre la membrane et le catalyseur Phénomène physique dû à des sollicitations mécaniques Conséquences Augmentation de la résistance entre les différentes couches, accumulation d’eau à la cathode et inondation de la pile Variables impliquées ηohm , ηconc Modèles intégrant [Al-Baghdadi 2009] ces variables Dégradation de la performance de la pile et dégradation du taux de Conséquences attendues sur les transfert de masse. performances et la durabilité. Apparition de trous dans la membrane Phénomène physique dû à des sollicitations mécaniques Conséquences Augmentation de la résistance de la membrane d’où une perte de la conductivité. 82 Création de poches emprisonnant l’eau d’où une inondation de la pile. Diminution de la perméabilité de la membrane au passage des ions H+ Variables impliquées ηohm , K Modèles intégrant ces variables Dégradation de la performance et diminution de la tension de sortie de Conséquences attendues sur les la pile. performances et la durabilité. Fuites de gaz Phénomène physique dû à des sollicitations mécaniques Conséquences Variation du flux des réactifs Variables impliquées ηconc , cO , c H 2 2 Modèles intégrant ces variables Dégradation de la performance. Conséquences attendues sur les performances et la durabilité. Perte d’intégrité des joints du stack Phénomène physique dû à des sollicitations mécaniques Conséquences Connexions desserrées Variables impliquées p H 2 , p O2 Modèles intégrant ces variables Une diminution de la pression augmente les pertes de transport de masse Conséquences attendues sur les d’où une dégradation de la performance de la pile performances et la durabilité. 83 Serrage des vis des plaques bipolaires Phénomène physique dû à des sollicitations mécaniques Conséquences Compression non homogène de la GDL Variables impliquées K , ηconc , rcont Modèles intégrant [Nitta 2007], [Hottinen 2007], [Zhang 2006], [P. Zhou 2006], [P. Zhou 2007] et [Fekrazad 2007] ces variables Une forte compression diminue la perméabilité des gaz, améliore la Conséquences attendues sur les conductivité et diminue les résistances de contact GDL/PB. Ceci fait performances et diminuer les pertes totales. la durabilité. Variation locale significative des transports de charge et de masse dans la GDL. Pression d’assemblage Phénomène physique dû à des sollicitations mécaniques Conséquences Résistance de contact entre la GDL et les plaques bipolaires. Différences de température dans la membrane Variables impliquées rcont Modèles intégrant [Zhang 2006], [Fekrazad 2007], [Sun 2005] et [Y. Zhou 2009] ces variables Une réduction de la résistance de contact entre la GDL et les plaques Conséquences attendues sur les bipolaires entraîne une augmentation de la performance du stack. performances et la durabilité. Tableau 2 - Tableau récapitulatif des différents phénomènes physiques pouvant influer sur la production de puissance d’une pile 84 Chapitre 2 Les outils pour la modélisation 2.1 Simulateurs systèmes 1D AMESIM DYMOLA 2.2 Simulation mécanique 2D ou 3D SAMCEF FLUENT 2.3 Simulation multiphysique 2D ou 3D COMSOL 2.4 Interface entre logiciels MATLAB 85 Dès le début de la thèse, il était clair que la simulation serait l’outil par excellence pour la découverte de l’influence de paramètres mécaniques multiphysiques sur les performances. Cependant, les logiciels à utiliser n’avaient pas encore été définis car les spécifications devaient être élaborées. Une étude sur les logiciels a été menée au début de la thèse pendant que parallèlement, des études étaient menées à FCLAB pour l’acquisition de logiciels en fonction des objectifs de l’institut. Finalement, il est apparu que les simulateurs 0D (sans dimension spatiale) comme AMESIM et SIMULINK étaient inadaptés pour le cœur de la thèse. Plusieurs outils logiciels sont cependant concernés par la thèse à différents niveaux. En effet tout d’abord, par delà l’objectif de connaître l’incidence sur les performances de paramètres mécaniques et multiphysiques, les résultats doivent pouvoir être intégrés dans les différents outils de simulation de l’institut FCLAB. En particulier les simulateurs 0D de type AMESIM ou DYMOLA étendus à une description 1D (en tenant compte d’une dimension, l’épaisseur de la GDL par exemple) doivent permettre d’optimiser les systèmes piles. Un travail de traduction des phénomènes étudiés en 3D et 2D au cours de la thèse doit être effectué sous forme d’élaboration et de détermination de paramètres pour les modèles 1D. D’autre part, les phénomènes étudiés par des simulations 2D et 3D nécessitent l’utilisation de logiciels. Ici le pragmatisme a prévalu. En ce sens que le logiciel COMSOL a été choisi car il correspond à nos attentes pour la modélisation multiphysique, dont c’est la raison d’être, et les calculs en mécanique des fluides et des structures. Cependant, pour d’autres objectifs, des collaborateurs des deux institutions, USEK et FCLAB possédaient une expérience avec d’un côté FLUENT pour les problèmes en mécanique des fluides et SAMCEF pour des problèmes liés à la mécanique structurale des piles à combustible. Ces compétences ont été exploitées pour avancer plus rapidement dans l’analyse des phénomènes concernés par la thèse, même si l’objectif final est de tout rassembler sous COMSOL par souci d’efficacité pour la suite. Notons aussi l’usage intensif de MATLAB, surtout à l’interface de tous ces logiciels. 86 Ce qui suit est une description résumée des logiciels déjà cités. 2.1 Simulateurs systèmes 1D Au début de la thèse, le choix de logiciel pour le simulateur 1D n’avait pas encore été effectué par FCLAB, mais AMESIM apparaissait comme le choix le plus probable, les développeurs basés en France ayant manifesté leur intérêt pour le développement de modules « Piles à Combustible ». Cette option ne s’est pas concrétisée pour des raisons stratégiques à FCLAB et finalement, c’est le logiciel DYMOLA qui a été choisi et mis en œuvre. AMESIM et DYMOLA ont des points communs : • Ce sont des concurrents sérieux de SIMULINK pour la simulation de systèmes dynamiques. • Ils ont des interfaces utilisateurs orientés objets, les objets devant être reliés entre eux conformément au système que l’on veut décrire. Mais contrairement à SIMULINK où les blocs sont plutôt des fonctions mathématiques, AMESIM et DYMOLA proposent des blocs physiques. AMESIM AMESIM est une plateforme de simulation de systèmes multi-domaines unidimensionnels. Elle utilise des symboles pour représenter les composants individuels du système. La modélisation est faite suivant les quatre étapes suivantes : • Sketch mode : les différents composants sont reliés entre eux. • Submodel mode : les sous-modèles physiques associés à chaque composant sont choisis. • Parameter mode : les paramètres de chaque sous-modèle sont fixés. • Run mode : la simulation est débutée. Un modèle AMESIM est formé de plusieurs sous-modèles reliés entre eux en respectant des règles de causalité. Chaque sous-modèle possède des ports qui peuvent recevoir plusieurs entrées et sorties. Ainsi, AMESIM est basé sur la théorie des Bond Graph. 87 L’outil Bond Graph (graphe de liaison ou graphe à liens), défini initialement par Paynter en 1959 aux Etats-Unis, est un outil pluridisciplinaire permettant de représenter graphiquement, à l’aide d’un langage unifié, la nature des échanges de puissance dans un système (apport, stockage et dissipation d’énergie). Dans la théorie Bond Graph, il s’agit de multiports qui sont les nœuds d’un graphe, de ports représentés par des segments et de liens qui relient une paire de multiports. Les liens sont des liens de puissance et représentent le flux d’énergie échangé entre les deux multiports connectés. Cette puissance s’exprime comme le produit de deux variables généralisées d’effort e et de flux f. Dans le formalisme des Bond Graph, l’énergie est un concept unifiant plusieurs domaines de la physique. En particulier, concernant les PAC qui sont des systèmes pluridisciplinaires par excellence, une approche énergétique via les Bond Graph englobera tous les domaines que couvre une PAC (électrochimique, électrique, thermique, hydraulique…). Par exemple [Saisset 2006], du point de vue hydraulique, la pression des gaz est l’effort et le débit des gaz le flux. La transition entre les domaines chimiques et électrochimiques est réalisée par un transformateur dont le rapport de transformation relie le potentiel de chaque électrode à l’énergie libre de Gibbs. Les différentes pertes d’activation et de diffusion sont aussi modélisées par un élément non linéaire dissipatif R. Le phénomène capacitif qui apparaît entre l’électrode et l’électrolyte est représenté par un élément C. Les pertes ohmiques sont représentées par un élément linéaire dissipatif. Cependant, cette modélisation n’est pas complète et beaucoup d’hypothèses ont été faites : Le modèle est unidimensionnel, à électrodes dissociées, les réactifs sont de l’hydrogène et du dioxygène purs, la concentration des gaz dans les canaux est uniforme et il n’existe pas de perte de pression et de réactions parasites, la diffusion des gaz est résolue en régime stationnaire. Ainsi, l’approche énergétique propose une description globale et non pas locale des phénomènes physiques. C’est pourquoi la modélisation des PAC par Bond Graph doit être incluse dans d’autres modélisations. 88 DYMOLA DYMOLA, Dynamic Modeling Laboratory, est un outil complet pour la modélisation et la simulation de systèmes complexes. Il permet aussi d’établir des interactions entre des systèmes de domaines différents comme l’électrique, mécanique, thermodynamique, pneumatique et thermique. Ainsi, les modèles décrivent la réalité physique de façon plus concrète. Comme DYMOLA utilise le langage MODELICA pour la description des objets, l’utilisateur peut créer ou modifier les modèles déjà existant selon le besoin. Des équations générales sont utilisées pour modéliser des phénomènes physiques. A l’institut FCLAB, des travaux de recherche utilisant DYMOLA / MODELICA ont permis de mettre à jour la bibliothèque relative aux PAC et ceci dans le but de modéliser une cellule de pile par une description 1D nécessaire à l’introduction des aspects mécaniques. Un objectif à terme des travaux menés dans cette thèse est de structurer les coefficients primaires en fonction des contraintes externes du genre mécanique. 2.2 Simulation mécanique 2D ou 3D SAMCEF Le logiciel SAMCEF est un logiciel de calcul des structures par éléments finis. Il est formé de plusieurs modules indépendants. SAMCEF Field est l’interface graphique dans laquelle l’utilisateur définit son modèle à partir d’une géométrie créée de façon assistée par ordinateur. Le module principalement utilisé ici est SAMCEF Linear dédié à l’analyse linéaire statique et dynamique des structures mécaniques. Le logiciel SAMCEF Field a été utilisé à FCLAB pour les modélisations mécaniques de cellules dans le cadre de divers projets. La thèse a pu profiter de résultats partiels pour le calcul des déformations de la GDL et de la PB quand la pile est comprimée par des charges uniformément réparties. 89 FLUENT FLUENT est un logiciel de mécanique des fluides basé sur la méthode des volumes finis pour simuler les écoulements de fluides et les transferts thermiques ainsi qu’une multitude de phénomènes associés incluant les écoulements turbulents, réactifs et multiphasiques. En analyse numérique, la méthode des volumes finis est utilisée pour résoudre numériquement des équations aux dérivées partielles en faisant des approximations d'intégrales. L’équation aux dérivées partielles est résolue de manière approchée sur un maillage constitué de volumes finis. Les volumes finis sont de petits volumes en 3D, des surfaces en 2D et des segments en 1D dont la réunion forme le domaine d'étude. Comme de nombreuses autres méthodes numériques, la qualité de la discrétisation est importante : unicité de la solution, stabilité, convergence et mesure d'erreur. Cette méthode permet donc de résoudre de manière discrète une équation aux dérivées partielles dont on cherche une solution approchée « suffisamment » fiable. Dans FLUENT, la modélisation consiste d’abord à dessiner ou importer la géométrie, définir ensuite les paramètres et les conditions du modèle et enfin traiter les résultats obtenus. Ce logiciel a été retenu pour l’étude de l’écoulement de gaz dans le canal des PB. 2.3 Simulation multiphysique 2D ou 3D COMSOL COMSOL Multiphysics (anciennement connu sous le nom de FEMLAB) est un logiciel de résolution basé sur une analyse par éléments finis pour différentes applications physiques et spécialement pour les phénomènes couplés ou multiphysiques. C’est un environnement interactif pour la modélisation et la résolution de problèmes scientifiques basés sur les équations aux dérivées partielles. Le logiciel résout des phénomènes physiques couplés et permet la construction de modèles en définissant les propriétés des matériaux, les charges, les contraintes, les sources et les flux plutôt que leurs équations grâce aux modes physiques implémentés. 90 COMSOL peut être utilisé dans plusieurs domaines d’application et comprend plusieurs modules (module chimique, module acoustique, module science de la terre, module transfert de chaleur, module mécanique…). Aussi, ce logiciel a été retenu pour la modélisation multiphysique d’un système de PEMFC. 2.4 Interface entre logiciels MATLAB MATLAB tient pour MATrix LABoratory. C’est un langage informatique de haut niveau et un environnement interactif pour le développement d’algorithmes, la visualisation et l’analyse des données ainsi que pour le calcul numérique. MATLAB permet la manipulation de matrices, d’équations différentielles, de polynômes, le traitement de signal ainsi que d’autres applications. Il permet aussi l’interfaçage avec des programmes écrits dans d’autres langages comme le langage C ou C++. MATLAB a été retenu pour tous les calculs numériques d’interpolation effectués dans notre travail de recherche. 91 Chapitre 3 Influence de la porosité et de la perméabilité locales de la couche de diffusion des gaz sur les performances des PEMFC 3.1 Approche générale 3.1.1 Le modèle partiel de représentation mécanique de la couche de diffusion 3.1.2 Modélisation multiphysique 3.2 Les modélisations proposées 3.2.1 Le modèle mécanique partiel pour le calcul des déformations 3.2.2 Le calcul des champs locaux de porosité et de perméabilité 3.2.3 Le modèle multiphysique initial 3.2.4 Le modèle multiphysique modifié 3.3 Les résultats de l'analyse numérique 3.3.1 Les déformations mécaniques de la PB et de la GDL 3.3.2 Les champs de porosité et de perméabilité de la GDL 3.3.3 Influence de la géométrie des plaques bipolaires en graphite 3.3.4 Influence de la compression dans le cas de plaques bipolaires en graphite 3.3.5 Influence de la géométrie des plaques bipolaires en acier 3.3.6 Comparaison des PB en graphite et en acier 3.3.7 Influence de la géométrie des PB en graphite sur la densité locale de courant 3.4 Conclusion 92 Pour bien répartir sur la membrane les gaz concentrés dans les canaux des plaques bipolaires, des couches de diffusion des gaz ont été intercalées. La porosité du matériau est nécessaire à la diffusion de gaz à travers la couche de diffusion. La qualité de cette diffusion de gaz est un élément duquel dépendent la densité de puissance et plus généralement la performance de la pile. Il est démontré ici que les charges mécaniques déforment la GDL. La déformation de cette dernière affecte ses propriétés physiques et conduit à des champs de porosité et de perméabilité hétérogènes. L’influence sur la performance d’une pile de ce manque d’uniformité pour la porosité et la perméabilité est l’objet d’étude dans ce chapitre. La variation du champ de porosité dépend de la géométrie de la plaque bipolaire qui est en contact avec la GDL. Dans cette étude, les trois types de section (rectangulaire, trapézoïdale et trapézoïdale avec congés de raccordement) et les deux types de matériau (graphite et acier) les plus courants ont été sélectionnés pour les canaux des plaques bipolaires (figures 3.1 et 3.2) qui sont considérées comme déformables. Ainsi, ce chapitre rend compte du couplage multiphysique entre les phénomènes physico-chimiques et le comportement mécanique des composants. Cette analyse s’effectue en tenant compte des conditions réelles d’usages d’une PEMFC. Dans les travaux suivants : [P. Zhou 2006 et P. Zhou 2007], [Y. Zhou 2009], [Zhang 2006], [Nitta 2007], [Fekrazad 2007], [Hottinen 2007], [Sun 2005], [AlBaghdadi 2007 et Al-Baghdadi 2009], la porosité de la GDL d’une PAC est gardée constante. Ici par contre, nous veillons à prendre en compte un champ de porosité local variable afin de simuler correctement l’effet de compression inhomogène dans la pile. Pour atteindre cet objectif, la méthode des éléments finis est utilisée pour modéliser et analyser les distributions locales de porosité et de perméabilité de la GDL. Différents niveaux de compression sur les plaques bipolaires sont étudiés. Les effets de ce champ de porosité et de cette perméabilité variables sur les courbes de polarisation et de densité de puissance sont rapportés. La densité de courant local est aussi étudiée. 93 3.1 Approche générale Les propriétés physiques et géométriques de la cellule étudiée sont dans les deux tableaux 3 et 4 suivants. Tableau 3 : Paramètres de transport électrochimiques du modèle Paramètre Symbole Valeur Unité Porosité initiale GDL ε 0.75 _ Perméabilité initiale de la GDL K 1.35e-10 m2 Conductivité ionique de la membrane km 17.12 S/m 0.2 _ Porosité de la membrane Coefficient de water drag drag H2O 3 _ Densité de courant d’échange à l’anode ref i 0,a 1000 A/m2 Densité de courant d’échange à la cathode ref i 0,c 1 A/m2 Coefficient de transfert symétrique α 0.5 _ Module d’Young de la GDL 60 MPa Coefficient de Poisson de la GDL 0.09 PB en graphite : module d’Young 10 PB en graphite : coefficient de Poisson 0.3 PB en acier : module d’Young 200 PB en acier : coefficient de Poisson 0.3 GPa GPa 94 Tableau 4 : Paramètres géométriques et d’opération Paramètre Symbole Valeur Unité Fraction massique initiale d’hydrogène ωH2 ,0 0.1 _ Fraction massique initiale d’oxygène ωO2 ,0 0.168 _ Fraction massique initiale d’eau ωH2O,c,0 0.2 _ Température d’entrée des carburants T 353 K Pression côté anode pa 1.1 atm Pression côté cathode pc 1.1 atm Viscosité dynamique de l’oxygène η 2.1e-5 Pa.s Epaisseur de la membrane δm 0. 1e-3 m Epaisseur de la GDL δGDL 0.2e-3 m Largeur du canal wc 1e-3 m Largeur de la dent wl 0.5e-3 m Hauteur du canal hc 0.5e-3 m Hauteur de la plaque bipolaire h BP 1.5e-3 m 3.1.1 Le modèle partiel de représentation mécanique de la couche de diffusion Le modèle partiel pour la représentation mécanique de la couche de diffusion est une structure 2D formée d’une moitié de dent et d’une moitié de canal côté cathode d’une PEMFC. La figure 3.1 représente des PB formées d’un solide plein en graphite. 95 Section rectangulaire Section trapézoïdale Section trapézoïdale avec un rayon de courbure de 0.2 mm Fig. 3.1 – Profils de dent de plaques bipolaires en graphite en contact avec la GDL. Dans la figure 3.2, la PB est formée de deux solides en acier (partie hachurée et partie pleine) dont un est en tôle mince (partie pleine). L’épaisseur des PB en tôle mince d’acier est de 0.07 mm et leurs dimensions externes sont les mêmes que celles des PB en graphite. Comme dans la figure précédente, la partie inférieure représente la membrane et celle du milieu la GDL. 96 Section rectangulaire Section trapézoïdale Section trapézoïdale avec un rayon de courbure de 0.2mm Fig. 3.2 – Profils de dent de plaques bipolaires en acier en contact avec la GDL. Ce modèle est développé pour étudier les effets de différentes pressions de compression sur la déformation de la GDL. Le domaine de calcul est ainsi limité car la pile présente une périodicité géométrique. La surface inférieure de la GDL est fixée et les côtés verticaux ont des conditions de symétrie. L’axe X est dans la direction transversale de la pile de façon à être transversal par rapport à l’écoulement dans les canaux des PB. L’axe Y est dans la direction longitudinale de la pile, c.-à-d. perpendiculaire aux PB. La porosité initiale de la GDL est prise égale à 0.75. C’est une valeur moyenne que nous avons adoptée dans notre étude. En effet, la porosité de la GDL de type TGP-H-120 est de 0.78 [Zhang 2006] ; de plus, dans [P. Zhou 2007], elle est de 0.7. Des simulations sont effectuées correspondant au plan d’expérience numérique suivant : Trois pressions uniformes (1 MPa, 5 MPa et 10 MPa) sont appliquées sur des plaques bipolaires en graphite et en acier et ceci pour différents types de dents : dent à section rectangulaire, dent à section trapézoïdale et dent à section trapézoïdale avec congés de raccordement (figures 3.1 et 3.2). Ainsi le plan d’expériences numériques se compose de 18 expériences. 97 Les pressions 5 MPa et 10 MPa bien que supérieures à celles dans la réalité facilitent la visualisation de la déformation de la GDL. Le paragraphe suivant présente les étapes de la modélisation multiphysique proposée. 3.1.2 Modélisation multiphysique La géométrie considérée est formée de la membrane et des GDL côté anode et côté cathode (figure 3.3). Par symétrie à notre modèle mécanique de déformation de la GDL, nous avons considéré une dent et deux moitiés de canaux de chaque côté de la dent. Les électrodes sont considérées comme des frontières. Partant du fait que le taux de la réaction d’oxydation d’hydrogène est très rapide, seule la déformation de la GDL côté cathode est modélisée. Cette GDL est comprimée sous la dent et gonfle sous les canaux. La cathode est alimentée par de l’oxygène, de l’azote et de l’eau. L’anode est alimentée par de l’hydrogène humidifié. Une différence de pression existe entre les canaux à gauche et à droite de la dent. Plusieurs autres hypothèses ont été prises en compte : • L’eau générée par la réaction de réduction de l’oxygène et par la migration d’eau de l’anode vers la cathode est liquide, • L’air et l’hydrogène sont des gaz idéaux, saturés en vapeur d’eau, • La GDL et le catalyseur sont homogènes et isotropes avant compression de la pile, • Les écoulements sont laminaires, à faible nombre de Reynolds, • La température est supposée constante. La figure 3.3 résume les conditions aux limites du modèle de pile PEMFC étudiée. 98 a- Entrée d’air humidifié; fractions massiques et pressions imposées. b- Sortie d’air humidifié; flux convectif à pression de référence. c- Entrée d’hydrogène humidifié; fraction massique d’hydrogène et pression imposées. d- Sortie d’hydrogène humidifié; flux convectif à pression de référence. Fig. 3.3 – Conditions aux limites du modèle étudié. Les courbes de polarisation et de densité de puissance ainsi que les densités de courant local sont obtenues grâce à un plan d’expériences numérique identique au précédent mais où les densités de courant sont obtenues par des simulations successives dans lesquelles la tension de la pile est fixée entre 0 et 1V avec des pas de 0.05V. 99 3.2 Les modélisations proposées 3.2.1 Le modèle mécanique partiel pour le calcul des déformations13 L’objectif de ce sous paragraphe est de définir la déformation de la GDL comprimée par une force répartie sur la PB. Le modèle en 3D comporte deux solides obtenus par l’extrusion des profils. L’un est la PB qui se présente comme une moitié du profil d’une dent et d’un canal pour des raisons de symétrie. L’autre solide représente la GDL et il est en contact avec le premier. Les trois géométries du profil et les deux matériaux de la PB introduits aux figures 3.1 et 3.2 sont étudiés. La figure 3.4 présente les modélisations éléments finis de deux types de sections pour les deux matériaux sélectionnés. La GDL est représentée dans la partie inférieure des figures et les PB (pleines en graphite et en tôle d’acier) dans la partie supérieure. 13 Ces modélisations ont été effectuées à l’aide du logiciel SAMCEF par Marie-Christine ILTCHEV en marge de son travail pour le projet ANR OPTIMECA. 100 PB en Graphite – Section rectangulaire PB en mince rectangulaire tôle d’acier – PB en Graphite – Rayon de courbure à 45 degrés Section PB en mince tôle d’acier – Rayon de courbure à 45 degrés Fig. 3.4 - Les modèles éléments finis relatifs aux différentes sections. Une série de pressions entre 1, 5 et 10 MPa est appliquée sur la face supérieure de la PB. La face inférieure de la GDL est bloquée. Les faces avant et arrière des solides sont aussi bloquées en profondeur. Des conditions de symétrie sont appliquées sur les faces verticales latérales de la PB et de la GDL. 101 La PB et la GDL sont assemblées par « contact ». Les faces du maillage en contact géométrique (avant application d’effort) sont prises en compte, mais également les faces susceptibles d’entrer en contact après déformation. Le modèle consiste en une seule rangée d’éléments 3D de type hexa-penta en forme de prismes ou de blocs. Des raffinements de maillage sont pratiqués dans les zones en contact et celles susceptibles d’entrer en contact après déformation car la GDL, comprimée sous la dent et gonflée sous le canal, peut coller sur le flanc de la dent des PB. La déformation est relevée sur les coordonnées des nœuds à la surface supérieure de la GDL. 3.2.2 Le calcul des champs de porosité et de perméabilité Un demi-modèle de GDL a été créé sous COMSOL Multiphysics. Il s’agit d’une structure 2D formée d’une moitié de dent et d’une moitié de canal du côté cathode d’une PEMFC. La déformation surfacique de la GDL est introduite dans COMSOL sous forme de déplacement pré-décrit. Le champ des déformations est ensuite calculé ev ( x, y ) = ex ( x, y ) + ey ( x, y ) + ez ( x, y ) (3.1) avec ev, la dilatation volumique, somme des déformations normales : e v = e x + e y + e z (3.2) Le champ de porosité est calculé selon la formule suivante (voir Annexe D): ε (x , y ) = ε 0 + ev (x , y ) 1 + e v (x , y ) (3.3) avec ε 0 la porosité initiale de la GDL avant compression. La perméabilité K est calculée d’après l’équation de Kozeny-Carman: [ε (x , y )]3 K (x , y ) = C [1- ε (x , y )]2 (3.4) 102 C est une constante donnée pour différents types de GDL. Elle vaut 1.276e-11 m 2 pour une GDL de type TGP-60-H et 3.952e-11 m 2 pour le type SGL31BA [Feser 2006]. Une valeur moyenne est C = 2e-11 m 2 . Par symétrie, les champs ont été étendus à une dent complète à l’aide de MATLAB. Une grille régulière de la dent complète est ensuite construite et les valeurs du champ de porosité aux nœuds de la grille sont calculées par interpolation linéaire. 3.2.3 Le modèle multiphysique initial Un modèle de cellule de PEMFC fondé sur les équations suivantes peut être trouvé dans COMSOL. Ces équations, bien connues des spécialistes de la PAC, gouvernent les différents phénomènes dans la pile [Comsol 3.5a]. La géométrie est celle de la figure 3.3 mais la GDL côté cathode n’est pas déformée. La porosité, la perméabilité et les coefficients de diffusion sont tous constants. Les symboles utilisés sont définis dans l’annexe A, à la fin de ce mémoire. Dans les GDL côté anode et côté cathode, la distribution du potentiel est donnée par : ∇.(- k e ∇V e ) = 0 (3.5) Dans la membrane, la distribution du potentiel est donnée par : ∇.(- k m ∇V m ) = 0 (3.6) La différence de potentiel entre la cathode et l’anode correspond au potentiel total de la cellule. Dans les GDL, le transport est modélisé par un transport en milieu poreux. L’équation de continuité dans les GDL est : ∇.( ρεu ) = 0 (3.7) La loi de Darcy modélise le flux dans les milieux poreux: u =− K η ∇p (3.8) La vitesse des gaz aux électrodes est donnée par les équations (3.9) et (3.10) suivantes : Du côté anode: 103 u =− ia M H2 ( ) ρF 2 (3.9) A l’anode, le mouvement de diffusion inverse (back diffusion) d’eau n’est pas pris en compte. Du côté cathode : u =− ⎤ i c ⎡ M O2 ⎛ 1 ⎞ − ⎜ + drag H 2O ⎟ M H 2O ⎥ ⎢ ρF ⎣ 4 ⎝2 ⎠ ⎦ (3.10) Au niveau du catalyseur, il ya consommation d’hydrogène et d’oxygène et production d’eau. Les flux massiques d’hydrogène à l’anode et d’oxygène et d’eau à la cathode sont alors: ia M H2 2F i N O2 = − c M O2 4F i 1 N H 2O = c ( + drag H 2O )M H 2O F 2 N H2 = − (3.11) A la cathode, la distribution de la densité de courant local ic dans la couche du catalyseur est donnée par l’équation de Butler-Volmer. ωO Fηact ,c exp(-α ) ωO ,0 RT i c = i 0,refc 2 (3.12) 2 avec α le coefficient de transfert symétrique. Il en est de même à l’anode mais il faut noter que la surtension d’activation à la cathode est généralement plus importante que celle à l’anode. L’équation de transport de masse dans les milieux poreux prend en compte les différentes espèces mises en jeu d’après l’équation de diffusion convection de MaxwellStephan: N ∇.[-ρωi ∑ D ij { j =1 ∇M M (∇ω j + ω j )} + ωi ρu ] = 0 Mj M (3.13) Les coefficients de diffusion effectifs, fonction de la porosité locale moyenne sont calculés d’après: 104 Deff O 2- N 2 = 0.22e - 4*(T / 293.2) ^1.5* ε ^1.5 Deff D O 2- H 2 O = 0.282e - 4 /(T / 308.1) ^1.5* ε ^1.5 eff H 2- H 2 O Deff = 0.915e - 4*(T / 307.1) ^1.5* ε ^1.5 (3.14) = 0.256e - 4*(T / 307.5) ^1.5* ε ^1.5 H 2O - N 2 3.2.4 Le modèle multiphysique modifié Dans notre étude, nous avons traité le même modèle précédent dans lequel nous avons modélisé la déformation de la GDL côté cathode et nous avons considéré des variables locales c.-à-d. qui dépendent de x et y, z étant l’axe de la pile. Les champs locaux de porosité (équation 3.3), de perméabilité (équation 3.4) et de coefficients de diffusion effectifs (équation 3.15) obtenus sont introduits dans le modèle multiphysique 2D. Rappelons que les coefficients de diffusion effectifs locaux, fonction de la porosité locale sont calculés d’après : Deff O 2-N 2 (x , y ) = 0.22e - 4*(T / 293.2) ^1.5* ε (x , y ) ^1.5 Deff O 2- H 2O (x , y ) = 0.282e - 4 /(T / 308.1) ^1.5* ε (x , y ) ^1.5 D eff H 2- H 2O (x , y ) = 0.915e - 4*(T / 307.1) ^1.5* ε (x , y ) ^1.5 Deff H 2O - N 2 (x , y ) = 0.256e - 4*(T / 307.5) ^1.5* ε (x , y ) ^1.5 (3.15) La résolution des équations se fait par la méthode des éléments finis avec le module COMSOL Multiphysics Chemical Engineering Module [Comsol 3.5a]. A chaque type de section de dent des PB et pour chaque pression de compression, des simulations sont faites pour différentes valeurs de la tension de la pile. La densité de courant local est calculée ce qui permet d’obtenir les courbes de polarisation et de puissance. La figure 3.5 montre le maillage pour une section de dent à 45° avec un rayon de courbure de 0.2 mm et pour une pression de 5 MPa. Il y a 1361 éléments et 9041 degrés de libertés. La taille maximale des éléments dans la zone raffinée du maillage est de l’ordre de 0.01 mm. 105 Fig. 3.5 – Exemple de maillage pour la détermination des champs de déformation, de porosité et de perméabilité dans la GDL (dent à 45° avec un rayon de courbure). 3.3 Les résultats de l'analyse numérique Comme rappel, le tableau 3 du paragraphe 3.1 de ce chapitre présente les valeurs des paramètres électrochimiques de transport du modèle étudié. Le tableau 4 du même paragraphe présente les paramètres géométriques et d’opération. 3.3.1 Les déformations mécaniques de la PB et de la GDL La figure 3.6 montre les déformations de la plaque bipolaire et de la GDL sous l’action d’une pression constante de 10 MPa exercée sur la face supérieure du modèle. Cette pression élevée permet de mettre en évidence la déformation de la surface de la GDL du côté de la plaque bipolaire et les zones de contact avec cette même PB (points de décollement). Ces résultats montrent aussi que dans le cas de PB en tôle, la GDL se déforme plus que dans le cas de PB en graphite et de façon non homogène. Par exemple, la déformation de la GDL pour des PB en graphite (section trapézoïdale avec congés de raccordement) sous une force de compression de 10 MPa est de 44.4 % alors que pour des PB en acier de même section, la déformation est de l’ordre de 50 %. D’autre part, les sections rectangulaires conduisent à un enfoncement de type coin dans la GDL, ce qui fait que le point de décollement se situe sur la partie verticale de la dent. Ces résultats sont utilisés pour la détermination des champs de porosité et de perméabilité de la GDL. 106 PB en Graphite – Section rectangulaire PB en Graphite – Rayon de courbure à 45 deg. PB en mince tôle d’acier – Section rectangulaire PB en mince tôle d’acier – Rayon de courbure à 45 deg. Fig. 3.6 – Déformations mécaniques pour différentes sections et matériaux (10 MPa). De plus, le tableau 5 résume les diminutions de l’épaisseur de la GDL sous la dent de graphite en fonction des pressions appliquées et des types de dent considérés. Pour une pression de 10 MPa et pour une dent à 90˚, l’épaisseur de la GDL sous la dent diminue considérablement (44 % environ) alors que pour une pression de 1 MPa et pour le même type de dent, l’épaisseur de la GDL sous la dent diminue d’environ 4.5 %. 107 Tableau 5 : Diminution de l’épaisseur de la GDL sous la dent en graphite Section de la PB Pression (MPa) 1 5 10 Rectangulaire 4.4% 22% 44% Trapézoïdale 4.5% 22.7% 45% (Fig. 3.7, figure supérieure) Trapézoïdale avec rayon 4.9% de courbure de 0.2 mm 3.3.2 23% 44.4% (Fig. 3.7, figure inférieure) Les champs de porosité et de perméabilité de la GDL La compression diminue localement l’épaisseur de la GDL, sa porosité et sa perméabilité varient en conséquent. La figure 3.7 est la GDL d’une PEMFC, constituée d’une moitié de dent et d’une moitié de canal côté cathode. Elle montre les champs de porosité associés à une pression de 5 MPa pour des plaques bipolaires en graphite. L’influence du rayon de courbure est marquée par une meilleure régularité dans la variation du champ local de porosité de la GDL. 108 Compression de 5 MPa– Section trapézoïdale Compression de 5 MPa– Section trapézoïde avec rayon de courbure Fig. 3.7 - Champs de porosité pour des dents en graphite. La figure 3.8 montre les champs locaux de porosité et perméabilité obtenus dans le cas d’une compression de 5 MPa et pour une section de dent à 45˚ des PB en graphite. Comme pour le champ de porosité, la force de compression appliquée fait varier localement la perméabilité de la GDL de façon proportionnelle à la porosité. 109 Champ de porosité Champ de perméabilité Fig. 3.8 - Comparaison des champs de porosité et de perméabilité. Compression de 5 MPa - PB en graphite - Dent à 45°. 3.3.3 Influence de la géométrie des plaques bipolaires en graphite Les courbes de polarisation et de densité de puissance ont été obtenues par des simulations successives relatives à chaque point de fonctionnement. Leur objectif est d’isoler l’influence des porosités et perméabilités locales. Elles ne tiennent donc pas compte de la variation des résistances de contact. Considérant des plaques bipolaires en graphite, les courbes pour les trois types différents de dents (à 90˚, à 45˚ et à 45˚ avec une courbure de 0.2 mm) et pour les valeurs de la pression de compression de 1 MPa, 5 MPa et 10 MPa ont été élaborées. La figure 3.9 met en évidence l’influence de la géométrie des PB. En effet, les paramètres géométriques (forme des dents des PB) ainsi que les forces de compression sur la pile influent sur son rendement à travers la variation locale des champs de porosité et de perméabilité. 110 La figure 3.9 montre qu’aux pressions élevées (10 MPa), la forme de la dent à 45˚ donne la meilleure performance alors qu’aux pressions inférieures, les courbes de polarisation et de densité de puissance sont très proches l’une de l’autre. 111 Courbes de polarisation Courbes de densité de puissance Fig. 3.9 – Courbes de polarisation et de densité de puissance pour les plaques bipolaires en graphite – Influence de la géométrie des plaques bipolaires. 112 3.3.4 Influence de la compression dans le cas de plaques bipolaires en graphite De plus, la compression de la pile rend la GDL plus fine sous la dent et plus gonflée sous le canal ; le volume de flux des gaz est alors plus restreint ce qui fait augmenter la vitesse des gaz et améliore le transfert de masse. Ainsi, les courbes de polarisation et de densités de puissance de la figure 3.10 montrent une amélioration de la performance de la pile aux pressions élevées. Ceci est valable pour toutes les formes de dent des PB. 113 Courbes de polarisation Courbes de densité de puissance Fig. 3.10 – Courbes de polarisation et de densité de puissance pour les plaques bipolaires en graphite – Influence de la compression. 114 3.3.5 Influence de la géométrie des plaques bipolaires en acier Le même travail a été effectué pour des plaques bipolaires en métal. La figure 3.11 montre les courbes de polarisation et de densité de puissance pour des PB en mince tôle d’acier, pour trois types de section de dents des PB et pour une charge de 5 MPa. Les tôles à 45° sont équivalentes par contre, les tôles avec des dents à section rectangulaire présentent des limites d’utilisation aux valeurs supérieures de la densité de courant. Ainsi, il existe une chute de performance à haute densité de courant pour les PB en tôle mince avec une section rectangulaire. Comme ce phénomène ne se produit pas pour les plaques en graphite de même section, nous supposons qu'il est dû aux grandes déformations mécaniques de la plaque. Quoiqu'il en soit, cela confirme le choix de section trapézoïdale pour les PB en mince tôle d’acier. Courbes de polarisation Courbes de densité de puissance Fig. 3.11 - Courbes de polarisation et de densité de puissance pour des PB en tôle mince d’acier – Influence de la géométrie des plaques bipolaires. 115 3.3.6 Comparaison des PB en graphite et en acier D’autre part, la figure 3.12 montre que, pour des plaques bipolaires à section trapézoïdale, au niveau des courbes de polarisation, il n’y a pas de différence entre les PB en graphite ou en métal du moins jusqu’à une compression de 5 MPa. Il en est de même pour les PB avec des congés de raccordement. Courbes de polarisation Courbes de densité de puissance Fig. 3.12 - Courbes de polarisation et de densité de puissance pour des PB en tôle d’acier et en graphite, à section trapézoïdale et pour une charge de 5MPa. A ce niveau, la distribution de la densité locale de courant peut être alors calculée. 3.3.7 Influence de la géométrie des PB en graphite sur la densité locale de courant La figure 3.13 montre, pour toutes les formes considérées de dents en graphite, pour une tension de 0.7 V et une pression de 5MPa, la distribution de la densité locale de courant sur deux interfaces : - l’interface entre la GDL et la PB - l’interface entre la GDL et la cathode. Le caractère non symétrique des courbes locales est dû aux gradients de pression existant entre les canaux à droite et à gauche de la dent. 116 Les résultats obtenus montrent que les régions de haute densité de courant sont les plus critiques à cause de la réduction de la diffusion d’oxygène suite à la compression de la GDL et car nous ne tenons pas compte de l’effet d’inondation à la cathode. De plus, les courbes inférieures de la figure 3.13 n'ont pas été lissées. Ce sont les résultats bruts du processus itératif de calcul dépendant de la taille du maillage et du seuil de convergence appliqué. Dans le cas d’une section de dent à 90° et à 5 MPa de compression, le décollement de la GDL se fait directement, le point de décollement étant le dernier point de contact entre la GDL non comprimée et la dent. Ainsi, les pics de courant apparaissent à ± 0.5 mm comme le montre la figure supérieure de la figure 3.13(c), à l’endroit où commence le canal. Toujours à 5 MPa, pour une dent à 45° avec ou sans rayon de courbure, la GDL reste collée sur le flanc du canal des PB et ne décolle pas directement. Ainsi, la dent paraît plus large et les pics de courant n’apparaissent qu’aux points de décollement de la GDL sur le flanc du canal comme le montrent les figures supérieures des figures 3.13(a) et (b). Ainsi une partie des flancs des canaux jouent le rôle de collecteurs de courant. Ces pics de courant montrent l’effet de la compression non homogène de la pile sur la distribution de la densité de courant. De plus, ces pics de courant peuvent avoir des conséquences sur la distribution de la température à l’intérieur de la pile. 117 118 a. Dent trapézoïdale b. Dent trapézoïdale avec rayon de courbure c. Dent rectangulaire Fig. 3.13 - Densité locale de courant – Dents en graphite – Compression de 5 MPa 119 120 3.4 Conclusion Dans ce chapitre, afin d’étudier les effets de la géométrie des dents des PB et de la compression de la GDL sur la performance de la pile, un modèle de déformation mécanique est d’abord présenté et ensuite un modèle multiphysique, prenant en compte des champs de porosité variables localement ainsi que des champs de perméabilité et des coefficients de diffusion effectifs variables localement eux aussi car ils sont fonction de la porosité. La performance de la pile est déduite d’après les courbes de polarisation et de densité de puissance pour différents types de dent et pour différentes pressions de compression. De plus, les courbes de densité locale de courant sont aussi étudiées. La méthode utilisée permet de prendre en compte la déformation des plaques bipolaires. Parmi les deux types de matériaux, seules les PB en tôles d’acier se déforment de façon visible. Par contre cette déformation n’a pas d’influence sur les courbes de polarisation, bien qu’il y ait un petit effet local de compression de la GDL. D’après les principaux résultats obtenus dans ce chapitre, nous déduisons l’influence du rayon de courbure de la dent de la PB qui est marquée par une meilleure régularité dans la variation du champ local de porosité et de perméabilité de la GDL. D’autre part, pour des PB en graphite, une amélioration de la performance de la pile est notée aux pressions élevées et la forme de la dent à 45˚ donne la meilleure performance. En ce qui concerne les PB en mince tôle d’acier, il a été démontré que les tôles avec des dents à section rectangulaire présentent des limites d’utilisation aux valeurs supérieures de la densité de courant. Ce phénomène ne se produit pas pour les PB en graphite ce qui confirme le choix de section trapézoïdale pour les PB en mince tôle d’acier. De plus, pour des PB à section trapézoïdale avec et sans congés de raccordement, il n’y a pas de différence, au niveau des courbes de polarisation, entre les PB en graphite ou en métal du moins jusqu’à une compression de 5 MPa. Pour ce qui est de la densité locale de courant, pour une dent trapézoïdale et une dent à 45° avec rayon de courbure, les pics de courant n’apparaissent qu’aux points de décollement de la GDL sur le flanc du canal alors que pour des dents rectangulaires, le point de décollement est le dernier point de contact entre la GDL non comprimée et la 121 dent. Ces pics de courant peuvent avoir des conséquences sur la distribution de la température à l’intérieur de la pile. Les résultats de ce chapitre ont fait l’objet d’un article [Akiki 2010] publié dans le journal Journal of Power Sources. Cependant, d’autres effets peuvent intervenir comme l’influence de la déformation des plaques bipolaires sur le comportement de la pile en général et sur sa performance et l’influence de la compression de la pile sur les flux des gaz dans les canaux et sur l’interface GDL/PB. Ce dernier point est étudié dans le chapitre suivant. 122 Chapitre 4 Influence sur les performances des PEMFC de la pression locale à l’interface entre la couche de diffusion des gaz et la plaque bipolaire 4.1 Approche générale 4.1.1 Le modèle partiel de représentation mécanique des fluides de l'écoulement de gaz 4.1.2 4.2 Modélisation multiphysique Les modélisations proposées 4.2.1 Le modèle mécanique des fluides pour la détermination de la pression locale à l'interface GDL/PB 4.2.2 Les équations du modèle multiphysique 4.3 Les résultats de l'analyse numérique 4.3.1 Le champ de pression locale à l'interface GDL/PB 4.3.2 Influence de la géométrie des plaques bipolaires en graphite 4.3.3 Influence de la compression dans le cas de plaques bipolaires en graphite 4.3.4 Influence d’une pression variable ou d’une pression moyenne à l’interface GDL /PB. Cas du graphite. 4.3.5 Influence de la géométrie des plaques bipolaires en acier 4.3.6 Influence d’une pression variable ou d’une pression moyenne à l’interface GDL /PB. Cas de l’acier. 4.3.7 4.4 Comparaison des PB en graphite et en acier Conclusion 123 Dans ce chapitre, une nouvelle approche numérique basée sur la mécanique des fluides est élaborée afin d’étudier l’effet de la compression mécanique de la pile sur la pression locale des gaz à l’interface entre la GDL et la PB. Nous démontrons au § 4.3.1 que cette pression locale d’interface n’est plus uniforme mais varie localement. Ensuite, le champ de pression locale d’interface GDL/PB est incorporé dans le modèle multiphysique de PEMFC (déjà présenté aux § 3.1.2 et 3.2.4) dans lequel les champs de porosité et de perméabilité varient localement afin de tirer les courbes de polarisation et de densité de puissance et de statuer sur la performance de la pile. Comme pour le chapitre précédent, cette étude s’intéresse aux PAC de type PEMFC pour le transport et pour trois types de section des canaux des PB (en graphite et en acier) considérés comme déformables. Dans les travaux suivants [P. Zhou 2006 et P. Zhou 2007], [Nitta 2007], [Fekrazad 2007], [Sun 2005] et [Al-Baghdadi 2007 et Al-Baghdadi 2009], des modèles de PEMFC sont développés afin d’étudier l’effet des forces de compression sur la pile mais la pression du gaz à l’entrée et à l’interface entre la GDL et la PB est gardée constante. Dans la suite, nous considérons un champ de pression d’interface GDL/PB variable localement afin de mieux simuler l’effet de compression de la PEMFC. 4.1 Approche générale 4.1.1 Le modèle partiel de représentation mécanique des fluides de l’écoulement de gaz On étudie l’écoulement fluide des gaz dans le serpentin de la PB à partir d’une structure mécanique 3D formée de la membrane, de la GDL déformable, d’une dent et de deux moitiés de canaux côté cathode d’une PEMFC. Ce modèle est développé pour tirer les valeurs du champ de pression locale à l’interface entre la GDL et la PB suite à la compression de la pile (figure 4.1). 124 Fig. 4.1 - Modèle 3D de pile pour la détermination de la pression locale à l’interface GDL/BP. En profondeur, nous avons considéré un tronçon de 1 cm et le domaine de calcul est ainsi limité car la pile présente une périodicité géométrique. La surface supérieure de la GDL (entre la GDL et la PB) est déformée suite aux différentes forces de compression appliquées et les côtés verticaux présentent des conditions de symétrie. Une approche de type mécanique des fluides est effectuée et le champ de pression d’interface GDL-PB obtenu est alors évalué et introduit dans le modèle multiphysique de pile afin de tirer les courbes de polarisation et de densité de puissance. Ceci est effectué dans le cadre du même plan d’expériences numériques formé des trois forces de compression appliquées (1 MPa, 5 MPa et 10 MPa) sur trois types de section des dents des PB (dent à section rectangulaire, dent à section trapézoïdale et dent à section trapézoïdale avec congés de raccordement) et pour deux types de PB (acier et graphite). 125 4.1.2 Modélisation multiphysique Les effets de la compression de la pile comme la déformation de la GDL, les champs de porosité et de perméabilité et le champ de pression locale d’interface GDL-PB variables localement sont pris en compte dans le modèle multiphysique 2D de la pile des § 3.1.2 et 3.2.4. La géométrie considérée est formée de la membrane et des GDL côté anode et côté cathode (figure 3.3). Les conditions et hypothèses du § 3.1.2 sont prises en compte également ici et les courbes de polarisation et de densité de puissance sont obtenues grâce au plan d’expériences numériques explicité dans le paragraphe précédent. 4.2 Les modélisations proposées 4.2.1 Le modèle mécanique des fluides pour la détermination de la pression locale à l’interface GDL/PB La déformation de la GDL suite à la compression exercée par la plaque bipolaire dont la face supérieure est chargée par une force répartie, a déjà été obtenue (voir § 3.3.1). Le champ des déformations est ensuite calculé et les champs de porosité et de perméabilité variant localement sont déduits (voir § 3.3.2 et [Akiki 2010]). Notre modèle est tridimensionnel et comporte trois solides obtenus par l’extrusion des profils d’une profondeur de 1 cm. Le premier solide est la PB qui se présente comme le profil d’une dent et de deux moitiés de canaux de part et d’autre de la dent, le second est la GDL déformée côté cathode et le troisième représente la membrane (figure 4.1). Trois géométries différentes de la dent de la PB sont étudiées : dent rectangulaire, dent à 45 degrés, dent à 45 degrés avec congés de raccordement pour deux types de PB : graphite et acier. Des raffinements de maillage sont pratiqués dans les zones critiques telles que l’angle sous le canal. Des conditions de symétrie sont appliquées sur les faces verticales latérales de la PB et de la GDL et une condition de pression sortante est appliquée sur une partie de la surface inférieure de la membrane pour simuler l’aspiration que crée la membrane. 126 Le gaz considéré ici est l’oxygène. Une différence de pression de 1000 Pa est prise en compte entre l’entrée et la sortie dans chaque partie de bras de serpentin. La simulation est réalisée avec le logiciel Fluent [Fluent 6.2] qui est un logiciel de mécanique des fluides basé sur la méthode des volumes finis pour simuler les écoulements de fluides. Il a été présenté au § 2.2. Le modèle 3D est statique et résolu par un algorithme à résolveur séparé. Dans les canaux de gaz, le transport se fait grâce à un gradient de pression ; les équations des fluides que constituent les équations de mouvement et de continuité de NavierStokes sont résolues : ρ .u .∇u = −∇p + η∇ 2u ∇.u = 0 (4.1) La membrane et la GDL sont des milieux poreux, l’équation de continuité est donnée par : ∇.( ρεu ) = 0 (4.2) Et la vitesse est obtenue par la loi de Darcy: u =− K η ∇p (4.3) Le champ de pression locale à l’interface entre la GDL et la PB est sauvegardé dans un fichier pour traitement ultérieur. Une grille régulière est construite à l’interface GDL/PB et les valeurs du champ de pression aux nœuds de la grille sont calculées par interpolation linéaire. Les fichiers d’interpolation de la pression locale d’interface obtenus ainsi que ceux de la porosité locale et de la perméabilité locale sont alors introduits dans le modèle multiphysique avec COMSOL. 4.2.2 Les équations du modèle multiphysique Les équations qui gouvernent les différents phénomènes dans la pile ont déjà été données aux § 3.2.3 et 3.2.4. Les symboles utilisés dans les équations sont définis dans la nomenclature en annexe A. 127 4.3 Les résultats de l'analyse numérique Les paramètres électrochimiques de transport et les paramètres géométriques et d’opération se trouvent dans les tableaux 3 et 4 au § 3.1. 4.3.1 Le champ de pression locale à l’interface GDL/PB Suite à la compression de la pile, la porosité et la perméabilité de la GDL varient localement et par conséquent la pression d’interface entre la GDL comprimée et la PB ne reste plus constante mais varie localement aussi. La figure 4.2 représente le modèle 3D formé de la membrane, de la GDL, d’une dent et de deux moitiés de canaux côté cathode et montre le champ de pression locale à l’interface GDL-PB associé à une force de compression de 5 MPa sur une dent trapézoïdale avec congés de raccordement d’une plaque bipolaire en graphite. Dans cette figure, la pression représente la pression relative (p-Patm). Fig. 4.2 - Champ de pression locale à l’interface GDL/PB sous une charge de 5 MPa et pour une section trapézoïdale avec un rayon de courbure de PB en graphite. L’échelle de la direction y est amplifiée 5 fois pour une bonne visualisation de la déformation de la GDL. 128 4.3.2 Influence de la géométrie des plaques bipolaires en graphite Les courbes de polarisation et de densités de puissance ont été obtenues suite à des simulations successives relatives à chaque point de fonctionnement où la tension de la pile varie de 0 à 1V avec des pas de 0.05V. L’objectif est d’isoler l’influence du champ de pression locale d’interface entre la GDL et la PB ainsi que celle des porosités et perméabilités locales. Considérant des plaques bipolaires en graphite, les courbes pour les trois types de dents (à 90˚, à 45˚ et à 45˚ avec une courbure de rayon 0.2 mm) et pour des forces de compression de 1 MPa, 5 MPa et 10 MPa sont données dans la figure 4.3. La pression locale d’interface GDL-PB est un champ variant localement tout comme la porosité et la perméabilité de la GDL. Nous notons une légère amélioration de la courbe de densité de puissance pour la dent trapézoïdale sous une forte charge. 129 Courbes de polarisation Courbes de densité de puissance Fig. 4.3 - Courbes de polarisation et de densités de puissance pour des PB en graphite – Influence de la géométrie des plaques bipolaires. 4.3.3 Influence de la compression dans le cas de plaques bipolaires en graphite Pour étudier les effets de la mise en charge, la figure 4.4 a été élaborée. Il s’agit de plaques bipolaires en graphite comprimées par 1 MPa, 5 MPa et 10 MPa. Ici aussi, la 130 pression locale à l’interface GDL-PB est un champ variant localement tout comme la porosité et la perméabilité de la GDL. La figure 4.4 montre une nette amélioration des courbes de polarisation et de densité de puissance aux pressions de compression élevées et ceci pour tous les types de dents. Ce résultat avait été déjà trouvé au chapitre précédent (voir § 3.3.4). 131 Courbes de polarisation Courbes de densité de puissance Fig. 4.4 - Courbes de polarisation et de densités de puissance pour des PB en graphite – Influence de la compression. Il peut donc être observé des figures 4.3 et 4.4 que les paramètres géométriques (forme des dents des PB en graphite) ainsi que les forces de compression influent sur le rendement de la pile à travers la variation locale du champ de pression d’interface GDLPB et des champs de porosité et de perméabilité. 132 4.3.4 Influence d’une pression variable ou d’une pression moyenne à l’interface GDL/PB. Cas du graphite. Il reste cependant à étudier l’effet d’un champ de pression d’interface GDL/PB variable ou d’une valeur moyenne de ce champ de pression d’interface. La figure 4.5 établit une comparaison entre deux modèles, l’un prend en compte un champ de pression d’interface variable localement et l’autre considère une pression constante de valeur égale à la moyenne du champ variable pour une dent trapézoïdale avec rayon de courbure des PB en graphite et sous une mise en charge de 1 MPa, 5 MPa et 10 MPa. Les profils sont similaires dans les deux cas d’étude ; c’est pourquoi, considérer la pression moyenne ou la pression variant localement n’influence pas de façon significative les courbes de polarisation et de densités de puissance. Cependant, il est évident que le calcul du champ de pression local est indispensable pour obtenir la valeur moyenne de la pression d’interface. 133 Courbes de polarisation Courbes de densité de puissance Fig. 4.5 - Courbes de polarisation et de densités de puissance pour des PB en graphite – Influence d’une pression variable ou d’une pression moyenne à l’interface GDL/PB. 4.3.5 Influence de la géométrie des plaques bipolaires en acier D’autre part, des simulations identiques ont été établies pour des PB en tôle mince d’acier. La figure 4.6 montre les courbes de polarisation et de densité de puissance pour 134 des PB en mince tôle d’acier, pour trois types de section et pour une mise en charge de 5 MPa avec la prise en compte d’un champ de pression d’interface GDL/PB variant localement. Les tôles à 45° sont équivalentes. Par contre, les tôles avec des dents à section rectangulaire présentent des limites d’utilisation aux valeurs supérieures de la densité de courant, ce qui diminue la performance des PB en tôle à section rectangulaire. Comme ce phénomène ne se produit pas pour les PB en graphite de même section, nous supposons qu'il est dû aux grandes déformations mécaniques de la plaque. Quoiqu'il en soit, cela confirme le choix de section trapézoïdale pour les PB en tôle mince. Courbes de polarisation Courbes de densité de puissance Fig. 4.6 - Courbes de polarisation et de densités de puissance pour des PB en tôle mince d’acier à 5 MPa. 4.3.6 Influence d’une pression variable ou d’une pression moyenne à l’interface GDL/PB. Cas de l’acier. De plus, l’effet de la considération, dans notre modèle multiphysique, d’un champ de pression d’interface GDL-PB variable ou d’une valeur moyenne de la pression d’interface sur des PB en tôle est aussi étudié. La figure 4.7 établit une comparaison entre deux modèles, l’un prend en compte un champ de pression d’interface variable localement et l’autre modèle considère une pression constante de valeur égale à la moyenne du champ variable pour les trois types de dents des PB en tôle et ceci sous une mise en charge de 5 MPa. 135 Là encore, comme pour les PB en graphite, une différence dans la performance de la pile n’est pas perceptible. Cependant, il est nécessaire de calculer le champ local de pression à l’interface GDL-PB afin d’en tirer sa valeur moyenne. 136 Courbes de polarisation Courbes de densité de puissance Fig. 4.7 - Courbes de polarisation et de densités de puissance pour des PB en tôle mince d’acier – Influence d’une pression variable ou d’une pression moyenne à l’interface GDL/PB. 137 4.3.7 Comparaison des PB en graphite et en acier Finalement, nous avons comparé les deux types de PB en graphite et en acier pour des dents rectangulaires, trapézoïdales et trapézoïdales avec congés de raccordement sous une force de compression de 5MPa. La pression locale d’interface entre la GDL et la PB est un champ variant localement d’où la figure 4.8. Les courbes sont presque identiques, au moins jusqu’à une compression de 5MPa. 138 Courbes de polarisation Courbes de densité de puissance Fig. 4.8 - Courbes de polarisation et de densités de puissance pour des PB en graphite et en acier à 5 MPa. 139 4.4 Conclusion Dans ce chapitre, différents modèles de pile dans lesquels la compression de la GDL est étudiée ainsi que les effets de cette compression sur la performance de la pile sont présentés. Les résultats typiques et les propriétés spécifiques à chaque modèle sont aussi donnés. Ainsi, la pression locale d’interface GDL/PB, la porosité, la perméabilité, la conductivité… sont influencées par des sollicitations mécaniques et ne peuvent plus être considérées constantes dans l’étude et la modélisation des PEMFC. Une modélisation complète de la PEMFC doit considérer des champs variant localement et non pas des valeurs constantes pour les propriétés ci-haut mentionnées. Afin d’étudier les effets de la géométrie des dents des PB et de la compression de la GDL sur la pression locale à l’interface GDL/PB et sur la performance de la pile, un modèle mécanique des fluides est d’abord présenté et ensuite un modèle multiphysique, prenant en compte des champs de pression locale d’interface GDL/PB, de porosité et de perméabilité variables localement, est développé. Le modèle considère aussi des coefficients de diffusion effectifs variables localement car ils sont fonction de la porosité et la méthode utilisée permet de prendre en compte la déformation des plaques bipolaires. Les modèles présentent des conditions de périodicité et de symétrie grâce à la géométrie du stack d’une PEMFC. La performance de la pile est déduite de par les courbes de polarisation et de densité de puissance pour différents types de dent et pour différentes pressions de compression. Ainsi, notre modèle multiphysique a été amélioré par l’introduction d’un champ de pression locale d’interface variant localement. Il a été démontré que cette distribution du champ de pression d’interface GDL/PB peut être remplacée par sa valeur moyenne sans changer de façon significative les courbes de polarisation et de densité de puissance cependant le calcul de ce champ local de pression d’interface est nécessaire pour tirer sa valeur moyenne. D’autres effets peuvent encore intervenir comme la modélisation, dans la résolution des équations du système ‘pile à combustible’, des électrodes et de la résistance de contact 140 sur l’interface entre la GDL et la PB pour approcher de plus en plus les conditions réelles d’utilisation de la pile. Le chapitre suivant traite l’effet de la compression de la pile sur la résistance de contact entre la GDL et la PB. 141 Chapitre 5 Influence sur les performances des PEMFC de la résistance de contact entre la couche de diffusion des gaz et la plaque bipolaire 5.1 Introduction 5.1.1 Phénomène de résistance de contact 5.1.2 Objectif du chapitre 5.2 Le modèle multiphysique global 5. 3 Les modélisations de la résistance de contact 5.3.1 Les types de modélisation 5.3.2 La modélisation choisie • Calcul des pressions • Calcul de la conductivité électrique • Introduction dans le modèle multiphysique global (distinction entre conductivité solide et conductivité de contact) 5.4 Les résultats de l'analyse numérique 5.4.1 Le champ de la résistivité électrique locale à l'interface GDL/PB 5.4.2 Influence de la compression dans le cas de PB en graphite 5.4.3 Comparaison de la résistivité de contact GDL/PB pour des PB en graphite et en acier 5.4.4 Comparaison de modèles avec ou sans le champ local de la résistivité de contact 5.5 • Cas du graphite • Cas de l’acier Conclusion 142 5.1 Introduction 5.1.1 Phénomène de résistance de contact L’interface entre deux matériaux est un contact plus ou moins important dépendamment de la pression de compression exercée. Plus la compression est grande, plus le contact est important et plus la résistance de contact diminue. Dans une PEMFC, plusieurs cellules sont connectées ensemble pour fournir une puissance suffisamment élevée. La GDL est en contact direct avec les PB. Cette résistance de contact à l’interface GDL/PB est un paramètre important pour le calcul des pertes et donc de la performance de la pile. 5.1.2 Objectif du chapitre Dans ce chapitre, les effets de la compression d’assemblage de la PAC sur la résistance de contact dans une cellule de pile sont étudiés à partir du modèle multiphysique, dans lequel les champs de porosité, de perméabilité et de pression à l’interface GDL/PB varient localement. Ce modèle est complété par l’introduction d’un champ local de la résistance de contact entre la GDL et la PB. Les effets de ces champs variant localement sur les courbes de polarisation et de densité de puissance sont rapportés et comparés avec deux modèles multiphysiques de pile dans lesquels le contact GDL/PB n’est pas pris en compte. Dans les travaux de [P. Zhou 2006 et P. Zhou 2007], [Y. Zhou 2009], [Zhang 2006], [Nitta 2007], [Fekrazad 2007], [Hottinen 2007], [Sun 2005], la résistance électrique de contact est modélisée suite à une compression inhomogène de la pile durant son assemblage. Il est démontré que cette résistance subit des modifications mais sa valeur est une valeur moyenne. Dans notre travail par contre, nous considérons un champ de résistance de contact variable localement afin de simuler correctement l’effet de compression inhomogène dans la pile. 143 5.2 Le modèle multiphysique global Les conditions de calcul et les hypothèses sont celles du chapitre 3. La géométrie considérée est celle de la figure 3.3. Une dent et deux moitiés de canaux de part et d’autre de la dent sont considérées et seule la déformation de la GDL du côté cathode est modélisée. Rappelons que les champs locaux de porosité et de perméabilité de la GDL, le champ local de pression à l’interface GDL/PB sont des grandeurs variables dans l’espace et non pas des valeurs moyennes. Il en est de même pour le champ de la résistivité électrique de la partie en contact entre la GDL et la PB. L’approche mécanique permettant d’obtenir les valeurs locales de cette résistivité est présentée dans le paragraphe suivant. Le même plan d’expériences décrit au chapitre 3 est appliqué. Différentes pressions uniformes (1 MPa, 5 MPa et 10 MPa) sont imposées aux plaques bipolaires déformables des figures 3.1 et 3.2. Les densités de courant sont obtenues par des simulations successives dans lesquelles la tension de la pile est fixée entre 0 et 1V avec des pas de 0.05V. Les équations qui gouvernent les différents phénomènes dans la pile ont déjà été données aux § 3.2.3 et 3.2.4. Les symboles utilisés dans les équations sont définis dans la nomenclature en annexe A. 5.3 Les modélisations de la résistance de contact 5.3.1 Les types de modélisation Dans la littérature, la résistance électrique de contact entre deux surfaces est considérée comme un circuit électrique formé de résistances en parallèle. La figure 5.1 montre la région de contact entre la GDL et la PB d’une cellule de PAC. RGDL est la résistance totale de la GDL, RPB celle des PB et RGDL / PB est la résistance de contact entre les deux domaines. 144 Fig. 5.1 – Région de contact à l’interface GDL/PB et résistance électrique équivalente [P. Zhou 2006]. Dans [P. Zhou 2006], l’auteur considère la résistivité électrique de contact ρ entre la GDL et la PB d’après la formule suivante : A⎛ B ⎞ ρ= ⎜ ⎟ L⎝T ⎠ C (5.1) avec : A, B et C des paramètres constants pour chaque type de PB et de GDL. L est la l’épaisseur de la région de contact et T (force divisée par une surface) est la contrainte locale de compression, grandeur pouvant être obtenue par la mise en œuvre du logiciel de calcul, ici COMSOL. 5.3.2 La modélisation choisie Il faut noter que COMSOL ne nous permet pas de donner directement la résistance de contact à l’interface entre deux corps. C’est pourquoi, il faut créer une paire de contact représentée par une couche très fine et lui affecter une conductivité électrique. Dans notre travail, nous avons choisi de modéliser le contact par l’introduction de la conductivité électrique (inverse de la résistivité électrique) d’une couche très fine. 145 • Calcul des pressions Nous avons défini au § 3.2.1 la déformation de la GDL comprimée par une force uniformément répartie sur la PB. Pour obtenir les contraintes locales à l’interface GDL/PB correspondantes, nous avions besoin d’un logiciel de calcul qui a été COMSOL dans notre cas. Les étapes donnant le champ de contraintes se résument par la définition d’une géométrie initiale (celle de la figure 3.7) puis par la définition de sa déformée (à l’aide de déplacements prédéfinis obtenus aussi au § 3.2.1). Le calcul numérique nous donne alors le champ local de pression à l’interface entre la GDL et la PB. • Calcul de la conductivité électrique Une fois que le champ local de contrainte est obtenu, la conductivité électrique σ est alors donnée par la formule suivante : σ= 1 A⎛ B ⎞ ⎜ ⎟ L⎝T ⎠ C (5.2) Pour des PB en graphite et une GDL de type GDL-10BA, on a d’après [P. Zhou 2006]: A = 3.32 mΩ.cm2 B = 1.01 MPa C = 0.534 Nous rappelons que L est la l’épaisseur de la région de contact et T le champ local de contrainte. Le fichier obtenu est ensuite traité avec MATLAB pour l’étendre au modèle de la figure 3.3 puis pour calculer, par interpolation linéaire, les valeurs sur chaque nœud de la grille régulière effectuée sur la GDL côté cathode. • Introduction dans le modèle multiphysique global (distinction entre conductivité solide et conductivité de contact) Le but est maintenant d’introduire les valeurs obtenues de conductivité dans le modèle multiphysique 2D de pile (voir figure 3.3). 146 Pour tenir compte de la conductivité solide de la GDL et de la conductivité du contact GDL/PB, nous considérons la dernière tranche de la grille régulière de la GDL (côté cathode) en contact avec la PB comme étant la fine couche qui modélise le contact GDL/PB. Il faut attribuer alors à cette fine couche les valeurs de résistivité telles que : ρ solide + ρcontact = • * ρ solide + ρ contact 2 * ρcontact (5.3) étant la résistivité de contact recherchée. * Ainsi donc, la résistivité de contact recherchée ρcontact vaut : * ρcontact = ρsolide + 2 ρcontact (5.4) ρsolide est l’inverse de la conductivité électrique de la GDL et ρcontact est obtenue d’après l’équation (5.1). La conductivité recherchée de la fine couche n’est autre alors que l’inverse * de ρcontact . 5.4 Les résultats de l’analyse numérique 5.4.1 Le champ de la résistivité électrique locale à l'interface GDL/PB La figure 5.2 représente la résistivité électrique de la partie de contact entre la GDL et la PB en graphite, obtenue d’après la formule (5.1) et associée à une force de compression de 10 MPa pour tous les types de dents des PB. L’axe des abscisses représente la largeur de la dent qui n’est autre que la région de contact entre la PB et la GDL comprimée (de x = 0 jusqu’aux points de décollement). La figure 5.2 montre une résistivité électrique plus petite sous la dent que sous le canal, pour toutes les sections des PB en graphite. Ceci est plus clair dans la figure 5.3 qui montre la courbe de la résistivité électrique sur toute la frontière entre la GDL et la PB pour une dent trapézoïdale en graphite sous 10 MPa. La résistivité électrique continue d’augmenter car sous le canal, le contact se perd. 147 Dent trapézoïdale – Dent rectangulaire – Dent trapézoïdale Compression de 10 MPa Compression de 10 MPa avec courbure – Compression de 10 MPa Fig. 5.2 - Courbes de la résistivité électrique de contact GDL/PB en graphite sous 10 MPa. 148 Fig. 5.3 - Courbe de la résistivité électrique sur toute la frontière GDL/PB pour une dent trapézoïdale en graphite sous 10MPa. 5.4.2 Influence de la compression dans le cas de PB en graphite Les courbes de la résistivité électrique de contact GDL/PB pour des dents trapézoïdales en graphite sous 1 MPa, 5 MPa et 10 MPa sont données à la figure 5.4. Cette figure montre la diminution de la résistivité électrique de la région de contact entre la GDL et la dent de la PB quand la pression de compression augmente. En effet, la compression de la pile améliore le contact entre la GDL et la dent de la PB et la résistivité électrique diminue alors. 149 Fig. 5.4 – Courbes de la résistivité électrique de contact GDL/PB pour des dents trapézoïdales en graphite – Influence de la compression. 5.4.3 Comparaison de la résistivité de contact GDL/PB pour des PB en graphite et en acier La figure 5.5 établit une comparaison entre la courbe de la résistivité électrique de contact pour une dent trapézoïdale en graphite sous 5 MPa et celle pour une dent trapézoïdale en tôle mince d’acier sous 5 MPa. La résistivité électrique de contact commence par être plus élevée dans le cas de l’acier mais devient par la suite plus petite aux points où les valeurs des abscisses sont proches du décollement. Ceci est dû au fait que la déformation des PB en tôle mince d’acier est plus grande que celle des PB en graphite (voir § 3.3.5). 150 Fig. 5.5 – Courbes de la résistivité électrique de contact GDL/PB pour des dents trapézoïdales en graphite et en acier sous 5 MPa. 5.4.4 Comparaison de modèles avec ou sans le champ local de la résistivité de contact. Les courbes de polarisation et de densités de puissance sont obtenues, comme dans les chapitres précédents, suite à des simulations successives relatives à chaque point de fonctionnement où la tension de la pile varie de 0 à 1 V avec des pas de 0.05 V. Ces courbes tiennent compte de la variation de la résistance de contact, de la porosité et de la perméabilité de la GDL ainsi que de la pression à l’interface GDL/PB. Nous remarquons que les courbes de polarisation et de densité de puissance obtenues pour un modèle multiphysique incluant le champ local de la résistivité électrique de contact GDL/PB et pour des modèles multiphysiques n’incluant pas ce champ variable sont quasi identiques. 151 En effet, les valeurs des résistivités de contact sont nettement plus petites que la résistivité électrique de la GDL. De plus, cette résistance de contact diminue localement si la pression locale de contact augmente. Nous avons étudié trois modèles multiphysiques : - Le premier modèle prend en considération des champs locaux de porosité, de perméabilité, de pression d’interface GDL/PB et de résistivité électrique de contact GDL/PB variables localement, - Le deuxième modèle considère des champs locaux de porosité, de perméabilité et de pression d’interface GDL/PB variables localement, - Le troisième modèle considère uniquement des champs locaux de porosité et de perméabilité variables localement. • Cas du graphite. La figure 5.6 montre les courbes de polarisation et de densité de puissance pour des PB en graphite sous une charge de 5 MPa pour les trois modèles explicités plus haut. Les courbes sont presque identiques. 152 Courbes de polarisation Courbes de densité de puissance Dent trapézoïdale Dent rectangulaire Dent trapézoïdale avec rayon de courbure Fig. 5.6 - Courbes de polarisation et de densités de puissance pour des PB en graphite sous une charge de 5 MPa. 153 • Cas de l’acier. Le même travail est fait pour des PB en acier. Les courbes obtenues pour les trois types de modèles étudiés sont identiques. Les sections rectangulaires présentent cependant des limites d’utilisation dans le cas des PB en mince tôle d’acier. Comme ce phénomène ne se produit pas pour les plaques en graphite de même section, nous supposons qu'il est dû aux larges déformations mécaniques des plaques en acier. 154 Courbes de polarisation Courbes de densité de puissance Dent trapézoïdale Dent rectangulaire Dent trapézoïdale avec rayon de courbure Figure 5.7 - Courbes de polarisation et de densités de puissance pour des PB en acier sous une charge de 5 MPa. 155 5.5 Conclusion Dans ce chapitre, nous avons étudié l’effet de la compression d’une PEMFC sur la résistance électrique de contact entre la PB et la GDL. Dans notre modélisation de la pile, les champs locaux de porosité et de perméabilité de la GDL, le champ de pression et de résistance de contact à l’interface GDL/PB varient localement. Les résultats des simulations ont montré que la compression de la pile fait diminuer la résistance de contact entre la GDL et la PB (voir figure 5.4). D’autre part, les courbes de polarisation et de densité de puissance obtenues pour les trois modèles suivants sont presque identiques : - Le premier modèle prend en considération des champs locaux de porosité, de perméabilité, de pression d’interface GDL/PB et de résistivité électrique de contact GDL/PB variables localement, - Le deuxième modèle considère des champs locaux de porosité, de perméabilité et de pression d’interface GDL/PB variables localement, - Le troisième modèle considère uniquement des champs locaux de porosité et de perméabilité variables localement. En effet, les pertes dues à la résistance de contact font partie des pertes ohmiques de la pile. Mais elles sont négligeables devant d’autres pertes comme la perte due à la résistance de l’électrolyte par exemple et qui est, elle, une source de limitation de courant (voir §1.2.3). De plus, nous n’avons étudié que les pressions de compression de 1 MPa, 5 MPa et 10 MPa. En réalité, il faudrait vraisemblablement descendre en dessous de 1 MPa de pression pour que la résistance de contact joue un rôle. Localement, cela peut se produire lors de vibrations. 156 Chapitre 6 Etude de sensibilité et de variabilité 6.1 Objectifs du chapitre 6.2 Présentation des résultats dans un plan d'expériences 6.3 Analyse de la variance 6.4 6.3.1 Influence sur la densité de puissance 6.3.2 Influence sur la densité de courant Régression multilinéaire 6.4.1 6.4.2 6.4.3 Coefficients des polynômes de la régression multilinéaire • Pour la densité de puissance • Pour la densité de courant • Pour la porosité Comparaison entre le modèle trouvé et l’expérience • Pour la densité de puissance • Pour la densité de courant • Pour la porosité Exemple de calcul pour la porosité 6.5 Incertitudes sur les modèles étudiés 6.6 Conclusion 157 6.1 Objectifs du chapitre Nous avons vu dans les chapitres 3, 4 et 5 l’influence de la compression d’une PEMFC sur les propriétés physiques de la GDL (à savoir sa porosité et sa perméabilité), sur la pression d’interface GDL/PB et sur la résistance de contact entre la GDL et la PB. La variation de ces grandeurs modifie les courbes de polarisation, les courbes de densité de puissance et la performance de la pile. La compression de la pile a été étudiée en considérant deux types de matériau des PB (le graphite et l’acier), trois géométries (rectangulaire et trapézoïdale avec et sans congés de raccordement) et trois pressions de compression (1 MPa, 5 MPa et 10 MPa). Apparaît alors l’importance d’une analyse d’influence et de sensibilité afin d’identifier les paramètres qui auront le plus d'influence sur les variables simulées. Cette étude peut s'avérer un outil fort utile à la prise de décision concernant la géométrie de la dent des PB, la nature des PB, … Elle est effectuée selon les méthodes des plans d’expériences. Ces méthodes nous donnent des surfaces de réponse sur lesquelles nous introduisons des incertitudes. Ainsi se fait la compréhension de l'influence des paramètres de conception (rayons de courbure, épaisseurs, paramètres liés au matériau, ...) et l’optimisation des paramètres physiques et géométriques de la pile pour en tirer la meilleure rentabilité possible tout en respectant les contraintes existantes. Ce chapitre présente donc une nouvelle vision des résultats de notre travail précédent de modélisation d’une cellule de pile. Les résultats obtenus seront rassemblés et analysés de façon à mettre en évidence l’influence et l’interaction des paramètres considérés. D’abord une analyse de la variance14 pour la détermination du pourcentage d’influence des variables considérées sur les paramètres analysés est effectuée. Ensuite, le calcul des coefficients de la régression multilinéaire15 est fait pour enfin introduire des incertitudes sur les variables étudiées. 14 L’analyse de la variance est une méthode statistique pour déterminer la variance de la moyenne de plusieurs ensembles d’observations. 15 La régression multilinéaire consiste à modéliser la relation entre deux ou plusieurs variables et une variable réponse en interpolant une équation linéaire aux données observées. 158 6.2 Présentation des résultats dans un plan d'expériences Pour cette analyse de sensibilité, nous avons considéré les résultats obtenus dans le chapitre 4. Il s’agissait de l’étude de l’influence sur les performances des PEMFC des champs locaux de porosité et de perméabilité de la GDL et d’une pression locale à l’interface entre la GDL et la PB. Le plan d’expériences numérique adopté est formé par les combinaisons des variables suivantes : • Le type du matériau des PB (Module d’Young) noté A • La géométrie des PB notée B • La pression d’assemblage notée C Les PB prises en compte dans nos modélisations sont soit en graphite soit en acier ; le paramètre A possède donc deux modalités. Les PB considérées présentent trois choix de géométrie ; dent à section rectangulaire, dent trapézoïdale avec ou sans congés de raccordement. Le paramètre B possède donc trois modalités. La pression d’assemblage exercée sur les PB est de 1 MPa, 5 MPa ou 10 MPa. Le paramètre C possède aussi trois modalités. La table suivante résume les modalités des différentes variables considérées. Modalité Module d’Young des PB Géométrie des PB Pression d’assemblage A B C 1 10 GPa (graphite) Rectangulaire 1 MPa 2 200 GPa (acier) Trapézoïdale 5 MPa 3 - Trapézoïdale avec 10 MPa courbure Ainsi, notre plan d’expériences consiste des 18 expériences suivantes : 159 Tableau 6 : Les 18 expériences du plan d’expériences numériques adopté No. de Module d’Young des PB Géométrie des PB Pression d’assemblage A B C 1 1 1 1 2 1 1 2 3 1 1 3 4 1 2 1 5 1 2 2 6 1 2 3 7 1 3 1 8 1 3 2 9 1 3 3 10 2 1 1 11 2 1 2 12 2 1 3 13 2 2 1 14 2 2 2 15 2 2 3 16 2 3 1 17 2 3 2 18 2 3 3 l’expérience 160 6.3 Analyse de la variance Afin d’étudier l’influence des différentes variables A, B et C, déjà définies, sur la densité de puissance et sur la densité de courant, nous avons établi les figures 6.1 et 6.2. 6.3.1 Influence sur la densité de puissance La figure 6.1 montre l’influence sur la densité de puissance. Les séries 1, 2, 3, 4, 5, 6 et 7 correspondent respectivement aux pourcentages d’influence du type de matériau A, de la géométrie des plaques bipolaires B, de la pression d’assemblage C, des interactions AB, AC et BC, et des résidus qui correspondent à l’interaction ABC. Fig. 6.1 – Pourcentage d’influence des variables et de leurs interactions sur la densité de puissance Nous notons l’influence de la pression d’assemblage sur la densité de puissance de la pile comme étant la plus importante et ceci pour toutes les tensions de la pile. Elle correspond à la série 3 (verte et triangulaire). Nous retrouvons le résultat de la figure 4.4 qui montrait l’effet de la compression sur les courbes de polarisation et de densité de 161 puissance. Cette influence de la pression d’assemblage est plus importante aux valeurs élevées de la tension de la pile qu’aux faibles valeurs de V. De plus, l’influence de toutes les autres variables (voir le reste des séries) est bien loin de celle de la pression d’assemblage C. 6.3.2 Influence sur la densité de courant Les séries 1 à 7 sont les mêmes que celles du paragraphe 6.3.1. Elles représentent les variables de type de matériau, de géométrie de la dent et de pression d’assemblage ainsi que leurs interactions. Fig. 6.2 - Pourcentage d’influence des variables et de leurs interactions sur la densité de courant Ici aussi, la figure 6.2 montre que l’effet de la pression d’assemblage (série 3) sur la densité de courant est le plus important. Résultat prévisible à partir de la figure 4.4. Aux tensions élevées de la pile, l’influence de la pression d’assemblage est plus grande. 162 6.4 Régression multilinéaire 6.4.1 Coefficients des polynômes de la régression multilinéaire Nous avons trois modalités pour chacune des variables B (géométrie de la dent) et C (pressions d’assemblage). Leurs effets sont donc non linéaires. Les modalités considérées sont telles que : Module d’Young des PB Géométrie des PB Pression d’assemblage A B C -1 10 GPa (graphite) Rectangulaire 1 MPa 0 - Trapézoïdale 5 MPa 1 200 GPa (acier) Trapéz. avec courbure 10 MPa Modalité Le plan d’expériences sera le suivant: Tableau 7 : Les modalités des 18 expériences du plan d’expériences numériques adopté No. de Module d’Young des PB Géométrie des PB Pression d’assemblage A B C 1 -1 -1 -1 2 -1 -1 0 3 -1 -1 1 4 -1 0 -1 5 -1 0 0 6 -1 0 1 7 -1 1 -1 8 -1 1 0 9 -1 1 1 10 1 -1 -1 l’expérience 163 11 1 -1 0 12 1 -1 1 13 1 0 -1 14 1 0 0 15 1 0 1 16 1 1 -1 17 1 1 0 18 1 1 1 La régression multilinéaire est donc de la forme : y = α 0 + α1A + α 2 B + α 3 B 2 + α 4C + α 5C 2 + α 6 AB + α 7 AB 2 +α 8 AC + α 9 AC 2 + α10 BC + α11BC 2 + α12 B 2C + α13 B 2C 2 • (6.1) Pour la densité de puissance La figure 6.3 montre les coefficients de la régression multilinéaire qui donnent la densité de puissance de la pile. Ainsi, nous pouvons calculer la densité de puissance de la pile pour n’importe quel module d’Young (compris entre 10 et 200 GPa) et n’importe quelle pression de compression (comprise entre 1 et 10 MPa). La géométrie des dents de la PB doit être obligatoirement choisie entre des dents à section rectangulaire et trapézoïdale avec ou sans congés de raccordement. 164 Fig. 6.3 – Coefficients des polynômes de la régression multilinéaire pour la densité de puissance de la pile Nous notons que le facteur de compression c fera augmenter la densité de puissance surtout aux faibles valeurs de la tension de la pile. • Pour la densité de courant La figure 6.4 donne les coefficients de la régression multilinéaire qui donnent la densité de courant de la pile en fonction des variables du plan d’expériences. 165 Fig. 6.4 – Coefficients des polynômes de la régression multilinéaire pour la densité de courant de la pile Nous notons que le facteur de compression c fera augmenter ici aussi la densité de courant surtout aux faibles valeurs de la tension de la pile. • Pour la porosité La figure 6.5 montre les coefficients de la régression multilinéaire qui donnent la porosité de la GDL. Nous avons considéré trois valeurs de la porosité qui correspondent à une coupe de la GDL à x = 0.5 mm, pour y = 0, y = -0,1 mm et y = -0,2 mm. Donc, la série ‘losange’ correspond à la porosité au point (0,5 ; 0), la série ‘carré’ à la porosité au point (0,5 ; -0.1) et la série ‘triangle’ à la porosité au point (0,5 ; -0,2). L’unité étant ici le millimètre. 166 Fig. 6.5 - Coefficients des polynômes de la régression multilinéaire pour la porosité de la GDL Dans les trois courbes, nous remarquons que la pression de compression c fera diminuer la porosité mais beaucoup plus pour y = 0 que pour y = -0.2 mm. 6.4.2 Comparaison entre le modèle trouvé et l’expérience • Pour la densité de puissance La figure 6.6 compare la densité de puissance obtenue par simulation et le modèle obtenu par régression multilinéaire pour une tension V = 0,5 V. 167 Fig. 6.6 – Comparaison de la densité de puissance entre le modèle et la simulation Le modèle linéaire est proche des valeurs obtenues par la simulation et les valeurs du modèle obtenues par régression linéaire ne sont pas dispersées. Les coordonnées x et y sont les valeurs de la densité de puissance. 168 • Pour la densité de courant Fig. 6.7 - Comparaison de la densité de courant entre le modèle et la simulation Ici aussi, les valeurs du modèle obtenues par régression linéaire sont autour des valeurs obtenues par simulation. Les coordonnées x et y sont les valeurs de la densité de courant. 169 • Pour la porosité Fig. 6.8 - Comparaison de la porosité entre le modèle et la simulation Les valeurs du modèle obtenues par régression linéaire ne sont pas dispersées et restent voisines des valeurs de porosité obtenues par simulation. Les coordonnées x et y sont les valeurs de la porosité. 6.4.3 Exemple de calcul pour la porosité Nous considérons l’exemple suivant : Pour le paramètre de type de matériau A, nous prenons une modalité a = 0.8. Cette modalité a = 0.8 correspond à un module d’Young E = 181.5 GPa environ. En effet, la relation linéaire pour la modalité a est : a = 0.0105*E - 1.105 (6.2) Cette relation est obtenue par la résolution d’un système de deux équations à 2 inconnues tel que : a = -1 pour E = 10 GPa et a = 1 pour E = 200 GPa. Pour la géométrie de la dent B, nous considérons la dent trapézoïdale donc b = 0. Pour la pression d’assemblage C, nous prenons c variant linéairement entre -1 et 1. Nous rappelons que la modalité c = -1, par exemple, correspond à la pression d’assemblage de 1 MPa. 170 La porosité est calculée sur la coupe x = 0.5 mm (aux alentours du point de décollement de la GDL) et pour différentes valeurs de y (y = 0, y = -0.1 mm et y = -0.2 mm). Fig. 6.9 – Exemples de régression linéaire pour la porosité en fonction de la pression d’assemblage La figure 6.9 montre comment la porosité évolue avec la pression de compression. La porosité est le plus influencée par la pression d’assemblage pour le modèle correspondant à la coupe y = 0. 6.5 Incertitudes sur les modèles étudiés Nous appliquons sur les paramètres a (type du matériau) et c (pression d’assemblage) des incertitudes qui modélisent les phénomènes aléatoires existant et qui influent sur le module d’Young et sur la pression d’assemblage. Ces incertitudes peuvent se produire par exemple lors de vibrations. Nous supposons que les paramètres a et c suivent une loi normale de moyenne 0 et d’écart-type 0,1. 171 Le paramètre b ne peut prendre que trois valeurs -1, 0 et 1 car nous avons uniquement étudié trois types de géométries bien discrètes (rectangulaire, trapézoïdale avec et sans congés de raccordement). Nous avons considéré un échantillon de 10 simulations pour la porosité à y = 0. La figure 6.10 montre la porosité en fonction de la géométrie de la dent, sachant que nous avons introduit des incertitudes sur le module d’Young et sur la pression d’assemblage. Les séries de 1 à 10 sont les 10 simulations de l’échantillon considéré. Fig 6.10 – Effets des incertitudes sur la porosité de la pression d’assemblage et du type de matériau des PB, à géométrie constante C’est pour la géométrie trapézoïdale avec congés de raccordement (b=1) que la porosité est la plus importante, ensuite la géométrie rectangulaire (b=-1) et finalement la géométrie trapézoïdale (b=0). Pour chaque type de géométrie, les valeurs de la porosité sont proches l’une de l’autre, sans grande dispersion. 172 6.6 Conclusion Ce chapitre a été consacré pour l’étude de sensibilité des paramètres de conception d’une pile PEMFC. Un plan d’expériences numériques formé de 18 expériences est établi. L’analyse de la variance, la régression multilinéaire et les incertitudes introduites sur les variables étudiées ont montré l’influence primordiale de la pression de compression sur la densité de puissance et la densité de courant de la pile. Ce résultat était déjà obtenu au paragraphe 4.3.3. Aussi, il a été démontré que le facteur de compression fera augmenter la densité de courant et de puissance surtout aux faibles valeurs de la tension de la pile. De plus, la porosité est plus influencée par la pression d’assemblage pour le modèle correspondant à la coupe y = 0 que pour y = -0.2 mm. Finalement, nous avons montré que la porosité est la plus importante pour la géométrie trapézoïdale avec congés de raccordement. Ce travail n’est en fait qu’une étude préliminaire et l’exploitation des résultats par les méthodes statistiques est vaste et ne peut être limitée à ce chapitre. Des perspectives de travaux ultérieurs sont explicitées dans la conclusion générale de ce mémoire de thèse. 173 Conclusion générale et perspectives L’objectif de ce travail a été d’étudier l’influence des sollicitations mécaniques sur les propriétés physiques des parties constituantes de la PAC et de montrer comment la variation des propriétés physiques (comme la porosité et la perméabilité de la GDL) et comportementales de la pile (comme la pression à l’interface GDL/PB et la résistance de contact entre la GDL et la dent de la PB) pouvait modifier la production de puissance de la pile et donc sa performance. De plus, l’étude numérique présentée dans ce mémoire a permis d’améliorer la connaissance et la compréhension du fonctionnement d’une pile à combustible de type PEMFC. D’une part, elle a apporté des précisions aux travaux antérieurs, et par la prise en compte des effets 3D, elle a permis d’évaluer la performance de la PAC par des simulations plus réalistes. Une analyse des différents paramètres de modèles dépendant des aspects mécaniques a été effectuée et les principaux paramètres à étudier dans le cadre de cette thèse ont été sélectionnés : porosité, perméabilité, coefficients de diffusion de la GDL, conductivité électrique du contact GDL/PB et volume des canaux après compression d’assemblage de la cellule. Les sections des canaux des PB les plus courantes ont été sélectionnées, associées à leur matériau : section rectangulaire pour les PB en graphite et section trapézoïdale avec et sans congés de raccordement pour les PB métalliques fines embouties. Un modèle partiel pour la représentation mécanique de la couche de diffusion d’une PEMFC du côté cathode a d’abord été mis en œuvre afin de déterminer la déformation de la GDL comprimée par une force répartie sur la PB. Les profils de déformation ont été obtenus pour les différents types de PB. Par la suite et à partir de ces profils, le champ local de contrainte mécanique sur la GDL a pu être déterminé ainsi que les zones de contact entre la GDL et la PB et les zones de décollement. Sur la base des contraintes mécaniques calculées dans la GDL, le champ local de porosité est obtenu. Nous avons procédé de même avec le champ local de résistance électrique de contact entre la GDL et la PB, lui aussi déterminé à partir du champ local de contrainte mécanique. 174 Enfin, une modélisation 3D de type volumes finis pour l’étude de la pression du fluide à l’interface entre la GDL et la PB a été élaborée. Elle est fondée sur les équations de mouvement et de continuité de Navier-Stokes dans les canaux de gaz de la pile, et dans la membrane et la GDL poreuses. La vitesse est obtenue par la loi de Darcy. L’analyse a permis de déterminer le champ local de pression d’oxygène sur la GDL du côté cathode. D’autre part, une modélisation multiphysique 2D d’une cellule de pile était disponible sur la base d’équations bien connues. La mise en charge de pressions de gaz est obtenue à partir d’une coupe effectuée sur les résultats de l’étude fluidique 3D. Les champs locaux de porosité de la GDL et de résistance électrique de contact GDL/PB sont aussi introduits dans ce modèle multiphysique. Une étude détaillée du comportement de la pile et de la modification de sa performance a pu être réalisée. Les résultats ont été présentés sous forme de courbes de polarisation et de densité de puissance. Concernant la porosité, on a pu observer que les PB en tôle se déforment plus que les PB en graphite. De plus, une nette amélioration des courbes de densité de puissance de la pile aux pressions élevées pour toutes les formes de dents des PB est visible et c’est la forme de la dent trapézoïdale qui donne la meilleure performance. Dans le cas de tôles avec des dents à section rectangulaire, elles présentent des limites d’utilisation aux valeurs supérieures de la densité de courant. D’où l’importance de considérer un champ local de porosité et non une valeur moyenne. Aussi, nous avons pu déduire l’influence du rayon de courbure de la dent de la PB qui est marquée par une meilleure régularité dans la variation du champ local de porosité et de perméabilité de la GDL. D’autre part, dans le cas d’une section de dent à 90°, le décollement de la GDL se fait directement sur le dernier point de contact entre la GDL non comprimée et la dent. Pour une dent à 45° avec ou sans congés de raccordement, la GDL reste collée sur le flanc du canal des PB et ne décolle pas directement. Ainsi, la dent paraît plus large et une partie des flancs des canaux jouent le rôle de collecteurs de courant. Concernant la pression de l’oxygène sur l’interface GDL/PB, comme pour la porosité variable, il existe une nette amélioration des courbes de polarisation et de densité de puissance aux pressions de compression élevées et ceci pour tous les types de dents avec une légère amélioration de la courbe de densité de puissance pour la dent trapézoïdale. Cependant, pour des PB en graphite ou en acier, les profils sont similaires 175 lorsque nous considérons une pression moyenne ou un champ local de pression à l’interface GDL/PB. Concernant la résistance de contact, il a été démontré que la compression de la pile fait diminuer la résistance de contact GDL/PB et les courbes de polarisation et de densité de puissance obtenues pour les trois modèles considérés (1- champs locaux de porosité, de perméabilité, de pression d’interface GDL/PB et de résistivité électrique de contact GDL/PB variables localement, 2- champs locaux de porosité, de perméabilité et de pression d’interface GDL/PB variables localement, 3- champs locaux de porosité et de perméabilité variables localement) sont presque identiques. Ceci est dû au fait que les pertes sont négligeables devant d’autres pertes limitatives du courant. Finalement tous les résultats ont été rassemblés pour une analyse de l’importance relative des différents paramètres sur la puissance de la pile. Il a été démontré que la pression de compression a le plus d’influence sur la densité de puissance et la densité de courant de la pile. Aussi, le facteur de compression fait augmenter la densité de courant et de puissance surtout aux faibles valeurs de la tension de la pile. De plus, la porosité de la GDL est plus influencée par la pression d’assemblage au niveau de la surface de contact GDL/PB qu’aux niveaux proches de la membrane. Finalement, nous avons montré que la porosité a les valeurs les plus élevées pour la géométrie trapézoïdale avec congés de raccordement. En perspective et suite à l’analyse des couplages multiphysiques incluant la mécanique (Ch. 1), l’influence de nombreux paramètres doit être analysée comme nous avons étudié l’influence de la porosité, de la pression à l’interface GDL/PB et de la résistance de contact GDL/PB. Il s’agit par exemple de : • la variation d’épaisseur de la membrane due à l’hydratation, • l’influence de la dilatation thermique des éléments de pile, … Dans cette thèse, nous avons étudié les variations de l’assemblage GDL/PB cathodique dans la cellule complète. Il faut encore lui ajouter l’étude des variations des couches électrolytiques, de la membrane et des éléments du côté anodique. Bien évidemment, il faudra étendre les études de la cellule à la pile complète. 176 Les champs obtenus dans cette étude sont généralement des champs 2D. Or, les modèles système utilisés à FCLAB sont des modèles 1D. La direction géométrique étant l’axe de la pile. Les coefficients de ces modèles peuvent donc dépendre seulement de cette direction. Il s’agira donc d’élaborer des procédures permettant d’identifier des champs 1D équivalents à l’action des champs 2D (voir figure I). Fig. I – Méthode d’identification du champ de porosité 1D de la GDL équivalent à l’action du champ 2D Il reste également beaucoup de travail à effectuer sur l’exploitation des résultats par les méthodes statistiques ou de fiabilité pour trouver des relations entre les différents paramètres en ce qui concerne leur influence sur les performances de la pile. Ces études sont des préliminaires aux études de durabilité faisant intervenir une modélisation élaborée du vieillissement sur les paramètres primaires de type surface 177 active, augmentation des résistances de contact par corrosion … Des études de type ACP16 (Analyse en Composantes Principales) doivent permettre de déterminer le niveau des interactions entre variables et paramètres pour piloter les modélisations système en introduisant des paramètres globaux et leur dépendance effective des conditions mécaniques. Par exemple, dans l’étude que nous avons faite, la puissance de la pile dépend des variables de pression d’assemblage, de géométrie et de type du matériau des PB. Ces variables influencent la puissance au travers des champs de porosité, de résistance de contact et de pression des gaz dans les canaux. Il s’agirait de voir l’importance relative des différents champs dans leur contribution à la variation de puissance. S’il s’avérait par exemple que le champ de pression de gaz ne contribuait pas à la variation de puissance suite à une modification de la géométrie des PB, il ne serait pas nécessaire de tenir compte de la pression des gaz dans la modélisation système de la fonction de transfert entre géométrie des PB et puissance. Par contre, si on remarquait une variation de puissance de 80% due à la contribution du champ de résistance de contact suite à une variation de la pression d’assemblage, il faudrait naturellement en tenir compte dans la modélisation de la fonction de transfert entre pression d’assemblage et puissance (voir figure II). 16 L’Analyse en Composantes principales (ACP) fait partie du groupe des méthodes descriptives multidimensionnelles appelées méthodes factorielles. L’ACP servira à mieux connaître les données sur lesquelles on travaille, à détecter éventuellement des valeurs suspectes, et aidera à formuler des hypothèses qu’il faudra étudier à l’aide de modèles et d’études statistiques [ACP 2006]. 178 Porosité Pression d'assemblage Pression des gaz Géométrie des PB Puissance de la pile Résistance de contact Matériau des PB Fig. II - Interaction entre les différentes variables et paramètres simulés sur la puissance et la densité de courant de la pile 179 Annexes A Nomenclature Alphabets c concentration (mol/l) Dij coefficient de diffusion (m2/s) drag H2O coefficient de la traînée d’eau : water drag ev dilatation volumique Ethermo potentiel thermodynamique (V) F constante de Faraday F = 96,487 (C/mole) I courant (A) ia densité de courant local à l’anode (A/m2) ref i 0,a densité de courant d’échange à l’anode (A/m2) ref i 0,c densité de courant d’échange à la cathode (A/m2) ic densité de courant local à la cathode (A/m2) ke conductivité électronique effective de l’électrode (S/m) km conductivité ionique de la membrane (S/m) k mem conductivité thermique de la membrane (W/mK) K perméabilité de la couche de diffusion (m2) M H2 masse moléculaire de l’hydrogène (kg/mole) M O2 masse moléculaire de l’oxygène (kg/mole) M H2O masse moléculaire de l’eau (kg/mole) N O2 flux massique d’oxygène N H2 flux massique d’hydrogène N H2O flux massique d’eau p pression (Pa) Patm pression atmosphérique (Pa) rcont résistance de contact 180 R constante universelle R = 8,314 (J/mole K) T température (K) u vitesse du flux (m/s) Vcell potentiel d’opération de la pile (V) Ve potentiel des électrodes Vm potentiel de la membrane Symboles grecs δGDL épaisseur de la couche de diffusion des gaz ε porosité de la GDL ρ masse volumique (kg/ m3) ωi fraction massique de l’espèce i η viscosité dynamique (Pa.s) ηact surtension d’activation (V) ηconc surtension de concentration (V) ηohm surtension ohmique due aux ions H+ dans la membrane et aux e- dans le circuit extérieur (V) α coefficient de transfert 181 B Liste des figures Figure 1.1 Eléments de modélisation d’une cellule. Figure 1.2 Caractéristique d’une PAC. Figure 1.3 Circuit équivalent et courbe de Nyquist de l’élément de Warburg. Figure 1.4 Circuit équivalent et courbe de Nyquist de l’élément de Warburg poreux. Figure 1.5 Circuit équivalent simple et courbe de Nyquist d’une PAC. Figure 1.6 Technique de mesure par interruption de courant. Figure 1.7 Technique de voltamétrie cyclique. Figure 1.8 Méthode de scellage. Figure 1.9 Localisation spatiale des phénomènes étudiés dans ce mémoire. Figure 3.1 Profils de dent de plaques bipolaires en graphite en contact avec la GDL. Figure 3.2 Profils de dent de plaques bipolaires en acier en contact avec la GDL. Figure 3.3 Conditions aux limites du modèle étudié. Figure 3.4 Les modèles éléments finis relatifs aux différentes sections. Figure 3.5 Exemple de maillage pour la détermination des champs de déformation, de porosité et de perméabilité dans la GDL (dent à 45° avec un rayon de courbure). Figure 3.6 Déformations mécaniques pour différentes sections et matériaux (10 MPa). Figure 3.7 Champs de porosité pour des dents en graphite. Figure 3.8 Comparaison des champs de porosité et de perméabilité. Compression de 5 MPa – PB en graphite - Dent à 45°. Figure 3.9 Courbes de polarisation et de densité de puissance pour les plaques bipolaires en graphite – Influence de la géométrie des plaques bipolaires. Figure 3.10 Courbes de polarisation et de densité de puissance pour les plaques bipolaires en graphite – Influence de la compression. Figure 3.11 Courbes de polarisation et de densité de puissance pour des PB en tôle mince d’acier – Influence de la géométrie des plaques bipolaires. Figure 3.12 Courbes de polarisation et de densité de puissance pour des PB en tôle d’acier et en graphite, à section trapézoïdale et pour une charge de 5MPa. 182 Figure 3.13 Densité locale de courant – Dents en graphite – Compression de 5 MPa. Figure 4.1 Modèle 3D de pile pour la détermination de la pression locale à l’interface GDL/BP. Figure 4.2 Champ de pression locale à l’interface GDL/PB sous une charge de 5 MPa et pour une section trapézoïdale avec un rayon de courbure de PB en graphite. L’échelle de la direction y est amplifiée 5 fois pour une bonne visualisation de la déformation de la GDL. Figure 4.3 Courbes de polarisation et de densités de puissance pour des PB en graphite – Influence de la géométrie des plaques bipolaires. Figure 4.4 Courbes de polarisation et de densités de puissance pour des PB en graphite – Influence de la compression. Figure 4.5 Courbes de polarisation et de densités de puissance pour des PB en graphite – Influence d’une pression variable ou d’une pression moyenne à l’interface GDL/PB. Figure 4.6 Courbes de polarisation et de densités de puissance pour des PB en tôle mince d’acier à 5 MPa. Figure 4.7 Courbes de polarisation et de densités de puissance pour des PB en tôle mince d’acier – Influence d’une pression variable ou d’une pression moyenne à l’interface GDL/PB. Figure 4.8 Courbes de polarisation et de densités de puissance pour des PB en graphite et en acier à 5 MPa. Figure 5.1 Région de contact à l’interface GDL/PB et résistance électrique équivalente. Figure 5.2 Courbes de la résistivité électrique de contact GDL/PB en graphite sous 10 MPa Figure 5.3 Courbe de la résistivité électrique sur toute la frontière GDL/PB pour une dent trapézoïdale en graphite sous 10MPa. Figure 5.4 Courbes de la résistivité électrique de contact GDL/PB pour des dents trapézoïdales en graphite – Influence de la compression. Figure 5.5 Courbes de la résistivité électrique de contact GDL/PB pour des dents trapézoïdales en graphite et en acier sous 5 MPa. Figure 5.6 Courbes de polarisation et de densités de puissance pour des PB en graphite sous une charge de 5 MPa. Figure 5.7 Courbes de polarisation et de densités de puissance pour des PB en acier sous une charge de 5 MPa. 183 Figure 6.1 Pourcentage d’influence des variables et de leurs interactions sur la densité de puissance Figure 6.2 Pourcentage d’influence des variables et de leurs interactions sur la densité de courant Figure 6.3 Coefficients des polynômes de la régression multilinéaire pour la densité de puissance de la pile Figure 6.4 Coefficients des polynômes de la régression multilinéaire pour la densité de courant de la pile Figure 6.5 Coefficients des polynômes de la régression multilinéaire pour la porosité de la GDL Figure 6.6 Comparaison de la densité de puissance entre le modèle et la simulation Figure 6.7 Comparaison de la densité de courant entre le modèle et la simulation Figure 6.8 Comparaison de la porosité entre le modèle et la simulation Figure 6.9 Exemples de régression linéaire pour la porosité en fonction de la pression d’assemblage Figure 6.10 Effets des incertitudes sur la porosité de la pression d’assemblage et du type de matériau des PB, à géométrie constante Figure I Méthode d’identification du champ de porosité 1D de la GDL équivalent à l’action du champ 2D Figure II Interaction entre les différentes variables et paramètres simulés sur la puissance et la densité de courant de la pile 184 C Liste des tableaux Tableau 1 Tableau récapitulatif des différents modèles étudiés Tableau 2 Tableau récapitulatif des différents phénomènes physiques pouvant influer sur la production de puissance d’une pile Tableau 3 Paramètres de transport électrochimiques du modèle Tableau 4 Paramètres géométriques et d’opération Tableau 5 Diminution de l’épaisseur de la GDL sous la dent en graphite Tableau 6 Les 18 expériences du plan d’expériences numériques adopté Tableau 7 Les modalités des 18 expériences du plan d’expériences numériques adopté 185 D Calcul du champ de porosité La porosité initiale ε 0 est égale au volume initial de vide V 0,vide sur le volume • total V total : ε0 = V 0,vide V total c.à.d. 1 = V 0,vide V 0,vide +V matiere = 1+ ε0 V matiere V 0,vide (1) (2) ce qui implique : V 0,vide V matiere • ε= = ε0 1− ε0 (3) Après compression, la porosité devient : V vide V vide = V total V vide +V matiere c.à.d. 1 ε = 1+ V matiere V vide et V vide =V 0,vide + ΔV vide (4) (5) (6) D’autre part, V total =V vide +V matiere (7) En supposant que V matiere est constant, on aura : ΔV total = ΔV vide (8) Alors l’équation 6 devient : V vide =V 0,vide + ΔV total • ev = (9) Comme la dilatation volumique e v vaut : ΔV total V total (10) L’équation 9 devient alors : V vide =V 0,vide +V total .ev =V 0,vide + (V 0,vide +V matiere ).ev =V 0,vide (1 + ev ) +V matiere .ev (11) ce qui donne d’après l’équation 3 : 186 V ε ε +e V vide = 0,vide (1 + ev ) + ev = 0 (1 + ev ) + ev = 0 v 1− ε0 1− ε0 V matiere V matiere • 1 ε = 1+ (12) Finalement d’après les équations 5 et 12, on obtient: V matiere 1− ε0 = 1+ ε 0 + ev V vide c.à.d. ε = ε 0 + ev 1 + ev (13) (14) Ainsi donc, le champ de porosité local est calculé selon la formule suivante : ε (x , y ) = ε 0 + ev (x , y ) 1 + e v (x , y ) (15) 187 Références [Fowler 2002] M. Fowler, R. Mann, J. Amphlett, B. Peppley, P. Rpberge, Incorporation of voltage degradation into a generalized steady state electrochemical model for a PEM fuel cell, J. Power Sources 106 (2002) 274-283. [Kundu 2006] S. Kundu, M. Fowler, L.C. Simon, S. Grot, Morphological features (defects) in fuel cell membrane electrode assemblies, J. Power Sources 157 (2006) 650-656. [P. Zhou 2006] P. Zhou, C.W. Wu, G.J. Ma, Contact resistance prediction and structure optimization of bipolar plates, J. Power Sources 159 (2006) 1115-1122. [P. Zhou 2007] P. Zhou, C.W. Wu, G.J. Ma, Influence of clamping force on the performance of PEMFCs, J. Power Sources 163 (2007) 874-881. [Y. Zhou 2009] Y. Zhou, G. Lin, A. J. Shih, S. J. Hu, Multiphysics Modeling of Assembly Pressure Effects on Proton Exchange Membrane Fuel Cell Performance, J. Fuel Cell Science and Technology 6 (2009) 041005 (7 pages). [Zhang 2006] L. Zhang, Y. Liu, H. Song, S. Wang, Y. Zhou, S. Jack Hu, Estimation of contact resistance in proton exchange membrane fuel cells, J. Power Sources 162 (2006) 11651171. [Nitta 2007] I. Nitta, T. Hottinen, O. Himanen M. Mikkola, Inhomogeneous compression of PEMFC gas diffusion layer, Part I. Experimental, J. Power Sources 171 (2007) 26-36. [Fekrazad 2007] N. Fekrazad, T.L. Bergman, The Effect of Compressive Load on Proton Exchange Membrane Fuel Cell Stack Performance and Behavior, J. Heat Transfer 129 (2007) 1004-1013. [Hottinen 2007] T. Hottinen, O. Himanen, S. Karvonen, I. Nitta, Inhomogeneous compression of PEMFC gas diffusion layer, Part II. Modeling the effect, J. Power Sources 171 (2007) 113-121. [Sun 2005] W. Sun, B. Peppley, K. Karan, Modeling the influence of GDL and flow-field plate parameters on the reaction distribution in the PEMFC cathode catalyst layer, J. Power Sources 144 (2005) 42-53. [Al-Baghdadi 2007] M. Sadiq Al-Baghdadi, H. Shahad Al-Janabi, Parametric and optimization study of a PEM fuel cell performance using three-dimensional computational fluid dynamics 188 model, J. Renewable Energy 32 (2007) 1077-1101. [Al-Baghdadi 2009] M. Sadiq Al-Baghdadi, A CFD study of hygro-thermal stresses distribution in PEM fuel cell during regular cell operation, J. Renewable Energy 34 (2009) 674-682. [Biyikoglu 2005] A. Biyikoglu, Review of proton exchange membrane fuel cell models, J. Hydrogen energy 30 (2005) 1181-1212. [O’Hayre 2006] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006. [O’Hayre 2006 - 0] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006, pp. 14. [Fuel Cell Handbook - 1] EG&G Services, Fuel Cell Handbook, Fifth Edition, Science Applications International Corporation, Parsons, 2000, pp. 2-19. [Fuel Cell Technology Handbook] Fuel Cell Technology Handbook, Hoogers, 2003. [O’Hayre 2006 - 1] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006, pp. 43. [Barbir 2005 - 1] F. Barbir, PEM Fuel Cells, Theory and Practice, Elsevier Academic Press, 2005, pp.18. [Noren 2005] D.A. Noren, M.A. Hoffman, Clarifying the Butler-Volmer equation and related approximations for calculating activation losses in solid oxide fuel cell models, J. Power Sources 152 (2005) 175-181. [O’Hayre 2006 - 2] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006, pp. 77. [Seong 2006] J.Y. Seong, Y.Ch. Bae, Y.K Sun, Water activities of polymeric membrane/water systems in fuel cells, SHORT COMMUNICATION, J. Power Sources 157 (2006) 733738. [Trabold 2006] T.A. Trabold, J.P. Owejan, D.L. Jacobson, M. Arif, P.R. Huffman, In situ investigation of water transport in an operating PEM fuel cell using neutron radiography: Part 1 – Experimental method and serpentine flow field results, International J. Heat and Mass Transfer 49 (2006) 47124720. [Chen 2005] D. Chen, H. Peng, A thermodynamic model of membrane humidifiers for PEM fuel cell humidification control, J. Dynamic Systems, Measurement and Control, 127 (2005) 189 424-432. [Barbir 2005 - 2] F. Barbir, PEM Fuel Cells, Theory and Practice, Elsevier Academic Press, 2005, pp.187. [Busquet 2004] S. Busquet, C.E. Hubert, J. Labbé, D. Mayer, R. Metkemeijer, A new approach to empirical modeling of a fuel cell, an electrolyser or a regenerative fuel cell, J. Power Sources 134 (2004) 41-48. [Xue 2006] X.D. Xue, K.W.E. Cheng, D. Sutanto, Unified mathematical modeling of steady-state and dynamic voltage-current characteristics for PEM fuel cells, Electrochimica Acta 52 (2006) 1135-1144. [Tao 2005] S. Tao, C. Guang-yi, Z. Xin-jian, Non linear modeling of PEMFC based on neural networks identification, J Zhejiang University Science (2005) 6A(5) 365-370. [Saisset 2006] R. Saisset, G. Fontes, Ch. Turpin, S. Astier, Bond graph model of a PEM fuel cell, J. Power Sources 156 (2006) 100107. [Haddad 2006] A. Haddad, R. Bouyekhf, A. El Moudni, M. Wack, Nonlinear dynamic modeling of proton exchange membrane fuel cell, J. Power Sources 163 (2006) 420-432. [O’Hayre 2006 - 3] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006, pp. 219. [O’Hayre 2006 - 4] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006, pp. 219. [O’Hayre 2006 - 5] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006, pp. 221. [O’Hayre 2006 - 6] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006, pp. 225. [O’Hayre 2006 - 7] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006, pp. 227. [Liu 2006] D. Liu, S. Case, Durability study of proton exchange membrane fuel cells under dynamic testing conditions with cyclic current profile, J. Power Sources 162 (2006) 521-531. [Kim 2005] S. Kim, S. Shimpalee, J. W. Van Zee, Effect of flow field design and voltage change range on the dynamic behavior of PEMFCs, J. Electrochemical Society 152 (6) (2005) A1265-A1271. 190 [Friede 2004] W. Friede, S. Raël, B. Davat, Mathematical Model and Characterization of the Transient Behavior of a PEM Fuel Cell, IEEE Transactions on power electronics 19 (5) (2004). [O’Hayre 2006 - 8] R. O’Hayre, S.W. Cha, W. Colella, F. Prinz, Fuel cell fundamentals, Wiley, 2006, pp. 254. [Baschuk 2004] J.J. Baschuk, Xianguo Li, Modelling of PEM fuel cells and stacks, Department of Mechanical Engineering, University of Waterloo, Waterloo, Ontario, Canada, April 2004. [Barreras 2005] F. Barreras, A. Lozano, L. Valiño, C. Marín, A. Pascau, Flow distribution in a bipolar plate of a proton exchange membrane fuel cell: experiments and numerical simulation studies, J. Power Sources 144 (1) (2005) 54-66. [Yen 2006] Ch. Yen, Sh. Liao, Y. Lin, Ch. Hung, Y. Lin, Ch. Ma, Preparation and properties of high performance nanocomposite bipolar plate for fuel cell, J. Power Sources 162 (2006) 309-315. [Blunk 2006 - 1] R. Blunk, F. Zhong, J. Owens, Automotive composite fuel cell bipolar plates: Hydrogen permeation concerns, J. Power Sources 159 (1) (2006) 533-542. [Blunk 2006 - 2] R. Blunk, M.H. Abd Elhamid, D. Lisi, Y. Mikhail, Polymeric composite bipolar plates for vehicle applications, J. Power Sources 156 (2) (2006) 151-157. [Wolf 2006] H. Wolf, M. Willert-Porada, Electrically conductive LCPcarbon composite with low carbon content for bipolar plate application in polymer electrolyte membrane fuel cell, SHORT COMMUNICATION, J. Power Sources 153 (2006) 41-46. [Antoni 2007] L. Antoni, J.P. Poirot-Crouvezier, F. Roy, X. GLIPA, Pile à combustible GENEPAC, 2007. [Fuel Cell Handbook - 2] EG&G Services, Fuel Cell Handbook, Fifth Edition, Science Applications International Corporation, Parsons, 2000, pp.10-18. [Blunier 2007] B. Blunier, A. Miraoui, Air Management in PEM Fuel Cells: State-of-the-Art and Prospectives, 2007. [Feser 2006] J.P. Feser, A.K. Prasad, S.G. Advani, Experimental characterization of in-plane permeability of gas diffusion layers, J. Power Sources 162 (2006) 1226-1231 [Comsol 3.5a] Comsol Multiphysics, Version 3.5a. Electrochemical Engineering: The Proton Exchange Membrane Fuel Cell. 191 [Akiki 2010] T. Akiki, W. Charon, M.C. Iltchev, G. Accary, R. Kouta, Influence of local porosity and local permeability on the performances of a polymer electrolyte membrane fuel cell, J. Power Sources 195(16) (2010) 5258-5268. [Fluent 6.2] Fluent 6.2, User’s Guide, January 2005. [ACP 2006] C. Duby, S. Robin, Analyse en Composantes Principales, Institut National Agronomique Paris – Grignon, Département O.M.I.P., 2006. 192