# Modélisation et simulation numérique des transferts fluidique et thermique dans le canal et couches cathodiques d'une PEMFC

Hocine Ben Moussa¹, Djamel Haddad¹, Kafia Oulmi²,
Bariza Zitouni¹, Bouziane Mahmah³ et Maiouf Belhamel³

¹Département de Mécanique, Faculté des Sciences de l'Ingénieur
Université Hadj Lakhdar, Rue Chahid M. El Hadi Boukhlouf, Batna, Algérie

²Département de Chimie, Faculté des Sciences
Université Hadj Lakhdar, Rue Chahid M. El Hadi Boukhlouf, Batna, Algérie

³ Division Bioénergie et Environnement, Centre de Développement des Energies Renouvelables,
B.P. 62, Route de l'Observatoire, Bouzareah, Alger, Algérie

(reçu le 12 Mars 2007 – accepté le 25 Mars 2007)

**Résumé** - A l'issue de cette étude, des phénomènes ayant lieu lors du fonctionnement d'une pile à combustible à membrane échangeuse de protons (PEMFC) ont pu être analysés. Un modèle, bidimensionnel transitoire de l'écoulement de gaz, et de transfert de chaleur dans la partie cathodique: canal, couche de diffusion, et couche de catalyseur, est présenté. Une approche bidimensionnelle a été choisie selon deux plans, longitudinal (PL) et transversal (P2). Les équations de conservation utilisées sont l'équation de continuité, l'équation de mouvement, et l'équation d'énergie. Ces équations ont été résolues numériquement en utilisant la méthode des volumes finis, sous un programme informatique (Fortran) développé localement. L'objectif de ce travail est l'obtention des profils et champs de vitesse des gaz dans le canal, la couche de diffusion (GDL) et la couche de catalyseur, ainsi que le profil et champ de température. On s'est intéressé à l'étude de l'effet des paramètres, perméabilité et nombre de Reynolds, sur l'écoulement, et l'effet des paramètres, perméabilité et la surtension cathodique sur les profils de température.

**Abstract** – A the exit of this study, the phenomena taking place during the operation of a fuel cell with exchanging membrane of protons (PEMFC) could be analyzed. A model, two-dimensional transitory of the flow of gas, and transfer of heat in the cathodic part: channel, layer of diffusion and lay down catalyst, is presented. A two-dimensional approach was selected according to two plans, longitudinal (PL) and transverse (P2). The conservation equations used are the equation of continuity, the equation of motion and the equation of energy. These equations were solved numerically by using the method of finished volumes, under a data-processing program (Fortran) developed locally. The objective of this work is obtaining the profiles and fields speed of gases in the channel, the layer of diffusion (GDL) and lay down it catalyst, as well as the profile and field of temperature. One was interested in the study of the effect of the parameters, permeability and Reynolds number, on the flow, and the effect of the parameters, permeability and the cathodic overpressure on the profiles of temperature.

**Mots clés:** PEMFC – Modèle – Bidimensionnel – Vitesse, température.

## 1 INTRODUCTION

Les piles à combustible à membrane échangeuse de protons (PEMFC) concentrent actuellement les plus gros efforts de recherche dans le monde. Elles fonctionnent à des températures relativement basses entre 80°C et 100°C et à des pressions de 1 à 5 bars. Les électrodes sont exposées à un flux de gaz, soit l'hydrogène (ou autre combustible) et l'oxygène (ou air ou autre comburant), elles sont généralement poreuses permettant la diffusion des gaz pour approvisionner en réactifs les zones actives où le catalyseur est en contact avec le conducteur ionique. Les couches de diffusion entourant les électrodes diffusent les gaz, permettent le transfert des électrons et assurent la gestion de l'eau pour une humidification optimale de la membrane et éventuellement une évacuation. Les plaques bipolaires permettent la distribution des gaz, l'évacuation de l'eau, la collecte du courant et assurent la connexion entre les piles élémentaires
qui composent la pile à combustible. Elles sont constituées de matériaux possédant un niveau de conductivité électrique suffisant et une imperméabilité au gaz supérieure à la membrane et doivent contenir des canaux de distribution de gaz généralement obtenus par usinage ou directement par moulage du polymère chargé sur une empreinte sans phase ultérieure d'usinage.

Cependant selon les applications visées, de nombreux problèmes de PEMFC restent à résoudre avant leur commercialisation. Ayant une compréhension des éléments qui affectent la performance de ce type de pile et de la façon dont les changements de conception affectent ces éléments, il est important de concevoir de meilleures piles à combustible. Un modèle de simulation est nécessaire qui peut être employé pour examiner des conceptions optimales. Actuellement, beaucoup de modèles existent et qui s'intéressent à cet effet. Mais, pour créer un outil qui est utile aux concepteurs, le modèle de simulation doit tenir compte de beaucoup de phénomènes physiques comme les transferts de chaleur et de masse, aussi bien les effets électriques et électrochimiques qui se produisent et doit tenir compte des composants importants de la pile à combustible, comme les plaques collectrices et la membrane. Le champ d'écoulement de réactifs dans les électrodes, les couches de catalyseur, et la membrane doivent être modelés ensemble de sorte que les effets des changements de conception soient reflétés à travers le modèle. Un modèle complet pourrait encore identifier les distributions uniformes de densité de courant et pourrait aider à identifier les propriétés matérielles qui sont importantes pour la performance de la pile à combustible.

Pour modéliser les piles à combustible de type PEMFC, il faut satisfaire les critères suivants du modèle qui doit être:

- basé sur une description mathématique de la physique de la pile, puis résolu en une, deux, et trois dimensions,

- utilisé comme outil d'étude parce qu'il facilite la compréhension des phénomènes physiques qui régissent le comportement de la pile,

- utile comme outil de conception parce qu'il tient compte de l'évaluation de la conception de la pile à combustible.

Ceci exige que le modèle soit flexible en ce qui concerne la géométrie et le matériel approprié. Les résultats de la modélisation mathématique de la pile à combustible à membrane échangeuse de protons (PEMFC) dépendent d'une série d'hypothèses. Ces dernières sont faites de telles façons à simplifier la résolution numérique et non pas pour l'exactitude de l'étude du modèle. Cependant, dans les modèles complets où on cherche à modéliser la pile entière, les hypothèses qui négligent le transport de l'eau liquide et de traiter les couches de catalyseur comme interfaces plutôt que comme volumes, diminuent considérablement l'exactitude et l'utilité du modèle.

Le modèle complet devrait inclure le transport le long de la direction d'écoulement afin d'étudier des effets importants comme l'épuisement de réactifs et la non uniformité de la densité de courant et la température de pile. L'utilisation des propriétés physiques correctes est aussi importante pour une nécessité mathématique, pour un modèle donné, plusieurs propriétés physiques sont exigées qui peuvent être déterminées expérimentalement. Les propriétés qui ne peuvent pas être déterminées expérimentalement sont typiquement traitées comme propriétés estimées, certaines sont ajustées pour avoir les résultats modélisés en accord avec les données expérimentales. Quand, le nombre et l'incertitude de propriétés estimées augmentent, l'exactitude des modèles diminue.

L'objectif de cette étude est la modélisation numérique bidimensionnelle transitoire de transfert de chaleur et de l'écoulement turbulent de l'oxygène humide afin d'obtenir des profils et champs de vitesse des gaz, et de température dans la partie cathodique (canal, couche de diffusion de gaz, couche de catalyseur). L'approche bidimensionnelle a été choisie selon deux plans, longitudinal (P1) et transversal (P2) (Fig. 1). Cette modélisation permet de visualiser la transition du milieu continu (canal) et le milieu poreux (couche de diffusion de gaz et couche de catalyseur), ainsi que l'effet de perméabilité, le nombre de Reynolds et la surtension cathodique sur l'écoulement, et la température.
## 2. MODELE MATHEMATIQUE

Dans ce modèle, les équations mathématiques de transport, les termes sources et les relations de fermeture sont présentés. La forme des équations peut changer d'un composant à un autre de la pile à combustible, c'est-à-dire chaque élément possède sa propre équation. Les canaux de gaz sont considérés comme un milieu continu, le régime turbulent d'écoulement est considéré et régi par l'équation de Navier Stockes. La couche de diffusion, la couche de catalyseur sont des milieux poreux; l'effet de la force exercée par le fluide sur le solide est caractérisé par un nouveau terme additionnel dans l'équation de mouvement. Le phénomène de turbulence est négligé dans ces milieux, car les vitesses d'écoulement sont faibles.

<figure><img src="image_1.png" /><figcaption>Fig. 1: Domaine d'étude</figcaption></figure>

### 2.1 Canaux de gaz

Les réactifs entrent et sortent de la pile à combustible par les canaux de gaz. Généralement, les canaux sont usinés dans les plaques collectrices (bipolaires ou mono polaires) et peuvent également être fabriqués par emboutissage ou moulage par injection. Le réseau entier du canal de gaz, est conçu de façon à distribuer uniformément les réactifs à travers la surface du GDL. Dans les canaux de gaz, l'équation de continuité, l'équation de Navier Stockes et celle de l'énergie, régissant les processus physiques s'écrivent:

**Equation de continuité**

$$ \frac{\partial \rho}{\partial t} + \frac{\partial \rho U}{\partial x} + \frac{\partial \rho V}{\partial y} = 0 \qquad (1) $$

**Equation de mouvement:** (projection sur x et sur y)

$$ \frac{\partial \rho U}{\partial t} + U \frac{\partial \rho U}{\partial x} + V \frac{\partial \rho U}{\partial y} = - \frac{\partial P}{\partial x} + \frac{\partial}{\partial x} \left( (\mu + \mu_t) \frac{\partial U}{\partial x} \right) + \frac{\partial}{\partial y} \left( (\mu + \mu_t) \frac{\partial U}{\partial y} \right) \quad (2) $$
$$ \frac{\partial \rho V}{\partial t} + U \frac{\partial \rho V}{\partial x} + V \frac{\partial \rho V}{\partial y} = - \frac{\partial P}{\partial y} + \frac{\partial}{\partial x} \left( (\mu + \mu_t) \frac{\partial V}{\partial x} \right) + \frac{\partial}{\partial y} \left( (\mu + \mu_t) \frac{\partial V}{\partial y} \right) \quad (3) $$

### Equation d'énergie

$$ \frac{\partial T}{\partial t} + U \frac{\partial T}{\partial x} + V \frac{\partial T}{\partial y} = \frac{\partial}{\partial x} \left( \left( \alpha + \frac{\mu_t}{Pr_t} \right) \frac{\partial T}{\partial x} \right) + \frac{\partial}{\partial y} \left( \left( \alpha + \frac{\mu_t}{Pr_t} \right) \frac{\partial T}{\partial y} \right) \quad (4) $$

Pour le calcul de la viscosité turbulente, on a choisi deux modèles de fermeture: modèle à zéro équation (longueur de mélange) et le modèle à deux équations ($k-\epsilon$).

- modèle à 0 équation (longueur de mélange)

$$ \mu_t = \rho l_m^2 \left| \frac{\partial U}{\partial y} \right| \quad (5) $$

avec $l_m$: longueur de mélange

- modèle à 2 équations ($k-\epsilon$)

$$ \mu_t = \frac{\rho C_{\mu} k^2}{\epsilon} \quad (6) $$

### Equation de bilan de l'énergie cinétique: k

$$ \frac{\partial \rho k}{\partial t} + U \frac{\partial \rho k}{\partial x} + V \frac{\partial \rho k}{\partial y} = - \frac{\partial}{\partial x} \left( \left( \mu + \frac{\mu_t}{\sigma_k} \right) \frac{\partial k}{\partial x} \right) + \frac{\partial}{\partial y} \left( \left( \mu + \frac{\mu_t}{\sigma_k} \right) \frac{\partial k}{\partial y} \right) + \mu_t \left( \frac{\partial U}{\partial y} \right)^2 - \rho \cdot \epsilon \quad (7) $$

### Equation de bilan de dissipation: $\epsilon$

$$ \frac{\partial \rho \epsilon}{\partial t} + U \frac{\partial \rho \epsilon}{\partial x} + V \frac{\partial \rho \epsilon}{\partial y} = \frac{\partial}{\partial x} \left( \left( \mu + \frac{\mu_t}{\sigma_\epsilon} \right) \frac{\partial \epsilon}{\partial x} \right) + \frac{\partial}{\partial y} \left( \left( \mu + \frac{\mu_t}{\sigma_\epsilon} \right) \frac{\partial \epsilon}{\partial y} \right) + \left( C_{1\epsilon} \mu_t \left( \frac{\partial U}{\partial y} \right)^2 - C_{2\epsilon} \rho \cdot \epsilon \right) \frac{\epsilon}{k} \quad (8) $$

Les constantes empiriques prennent les valeurs suivantes pour les écoulements simples (écoulement en canal) [14]: $C_\mu = 0.09$, $C_{1\epsilon} = 1.44$, $C_{2\epsilon} = 1.92$, $\sigma_k = 1$, $\sigma_\epsilon = 1.3$.

## 2.2 Canaux de diffusion de gaz (GDL)

Dans les diffuseurs de gaz, afin de tenir compte de la nature poreuse, et d'employer des formes semblables d'équations, U et V sont les composantes de la vitesse dans les pores.

### Equation de continuité

$$ \frac{\partial \epsilon_d \rho}{\partial t} + \frac{\partial \rho \epsilon_d U}{\partial x} + \frac{\partial \rho \epsilon_d V}{\partial y} = 0 \quad (9) $$

### Equation de mouvement: (projection sur x et sur y)

$$ \frac{\partial \rho \epsilon_d U}{\partial t} + U \frac{\partial \rho \epsilon_d U}{\partial x} + V \frac{\partial \rho \epsilon_d U}{\partial y} = - \epsilon_d \frac{\partial P}{\partial x} + \frac{\partial}{\partial x} \left( r \cdot \mu \frac{\partial U}{\partial x} \right) + \frac{\partial}{\partial y} \left( r \cdot \mu \frac{\partial U}{\partial y} \right) - \epsilon_d^2 \frac{\mu}{k_h} U \quad (10) $$

$$ \frac{\partial \rho \epsilon_d V}{\partial t} + U \frac{\partial \rho \epsilon_d V}{\partial x} + V \frac{\partial \rho \epsilon_d V}{\partial y} = - \epsilon_d \frac{\partial P}{\partial y} + \frac{\partial}{\partial x} \left( r \cdot \mu \frac{\partial V}{\partial x} \right) + \frac{\partial}{\partial y} \left( r \cdot \mu \frac{\partial V}{\partial y} \right) - \epsilon_d^2 \frac{\mu}{k_h} V \quad (11) $$
où les termes sources de l'équation de mouvement expliquent la résistance de l'écoulement (c'est-à-dire le terme de Darcy).

## Equation d'énergie

$$ \epsilon_d \rho C_p \left( \frac{\partial T}{\partial t} + U \frac{\partial T}{\partial x} + V \frac{\partial T}{\partial y} \right) = k^{\text{eff}} \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) \quad (12) $$

Les nombres adimensionnels Prandtl et Schmitt sont calculés avec des valeurs effectives:

$$ Pr_t = \frac{\mu C_p}{k^{\text{eff}}}, $$

$$ Sc_t = \frac{\mu}{\rho D_k^{\text{eff}}} $$

La conductivité thermique effective du milieu poreux est [3]:

$$ k^{\text{eff}} = -2k_{\text{gr}} + \frac{1}{\frac{\epsilon_g}{2k_{\text{gr}} + k_{\text{gas}}} + \frac{1 - \epsilon_g}{3k_{\text{gr}}}} $$

La diffusivité effective du milieu poreux est [6]:

$$ D_k^{\text{eff}} = D_k \epsilon_d^{1.5} $$

Le facteur de correction du milieu poreux r dans l'équation de quantité de mouvement a comme expression [3]:

$$ r = 2.25 \frac{(1 - \epsilon_d)^2}{\epsilon_d^2} $$

## 2.3 Couche de catalyseur

Dans cette zone, on définit les termes sources dans les équations d'énergie et des espèces qui expliquent le phénomène de la réaction électrochimique. L'équation de continuité et l'équation de mouvement restent applicables, dans cette couche de catalyseur.

La porosité effective de la couche de catalyseur peut être calculée en fonction de la porosité de la membrane $\epsilon_m$ et la fraction de volume de la membrane dans la couche de catalyseur $\epsilon_{mc}$.

### Equation de continuité

$$ \frac{\partial \epsilon_c \rho}{\partial t} + \frac{\partial \rho \epsilon_c U}{\partial x} + \frac{\partial \rho \epsilon_c V}{\partial y} = 0 \quad (13) $$

### Equation de mouvement: (projection sur x et sur y)

$$ \frac{\partial \rho \epsilon_c U}{\partial t} + U \frac{\partial \rho \epsilon_c U}{\partial x} + V \frac{\partial \rho \epsilon_c U}{\partial y} = -\epsilon_c \frac{\partial P}{\partial x} + \frac{\partial}{\partial x} \left( r \cdot \mu \frac{\partial U}{\partial x} \right) + \frac{\partial}{\partial y} \left( r \cdot \mu \frac{\partial U}{\partial y} \right) - \epsilon_c^2 \frac{\mu}{k_h} U \quad (14) $$

$$ \frac{\partial \rho \epsilon_c V}{\partial t} + U \frac{\partial \rho \epsilon_c V}{\partial x} + V \frac{\partial \rho \epsilon_c V}{\partial y} = -\epsilon_c \frac{\partial P}{\partial y} + \frac{\partial}{\partial x} \left( r \cdot \mu \frac{\partial V}{\partial x} \right) + \frac{\partial}{\partial y} \left( r \cdot \mu \frac{\partial V}{\partial y} \right) - \epsilon_c^2 \frac{\mu}{k_h} V \quad (15) $$

$$ \frac{\partial \epsilon_c C_k}{\partial t} + U \frac{\partial \epsilon_c C_k}{\partial x} + V \frac{\partial \epsilon_c C_k}{\partial y} = \frac{\partial}{\partial x} \left( D_k^{\text{eff}} \frac{\partial C_k}{\partial x} \right) + \frac{\partial}{\partial y} \left( D_k^{\text{eff}} \frac{\partial C_k}{\partial y} \right) + S_k \quad (16) $$

### Equation d'énergie

$$ \epsilon_c \rho C_p \left( \frac{\partial T}{\partial t} + U \frac{\partial T}{\partial x} + V \frac{\partial T}{\partial y} \right) = k^{\text{eff}} \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) + S_T \quad (17) $$
Le terme source de chaleur dans l'équation d'énergie est:

$$S_T = i \eta + \frac{i^2}{\sigma_{ct}}$$

où le premier terme est la chaleur produite par la surtension d'activation et le second terme est la chaleur produite par effet joule.

### 3. RESOLUTION NUMERIQUE

L'outil numérique utilisé pour résoudre nos systèmes d'équations est la méthode des volumes finis. La technique de cette méthode consiste à discrétiser la forme intégrale des équations de conservations plutôt que leurs formes différentielles. La formulation différences finies adoptée ici, consiste en une approche du type volume de contrôle. Son principe est basé sur le fait de proposer un champ de pression qui permet d'obtenir une première approximation du champ de vitesse de l'écoulement. Puis, on applique des corrections du champ de pression dans le but de rendre le champ de vitesse résultant conforme avec l'équation de continuité. Cette procédure a été réexaminée par la suite et adaptée à de multiples problèmes de transport-diffusion comme le notre. Les équations peuvent s'écrire sous la forme générale suivante:

$$\frac{\partial}{\partial t}(\varepsilon \rho \phi) + \text{div}(\varepsilon \mu u \phi) = \text{div}(\Gamma \cdot \text{grad}(\phi)) + S_{\phi} \quad (18)$$

L'évaluation des différentes intégrales de l'équation de transport-diffusion ci-dessous demande que l'on fasse un certain nombre d'approximations de façon à former un système d'équations linéaires liant entre elles les valeurs de la fonction $\phi$ aux points du maillage $P_{IJ}$.

$$\int_{VC}^{t+\Delta t} \int_{t} \left( \frac{\partial}{\partial t}(\varepsilon \rho \phi) + \frac{\partial}{\partial x}(\varepsilon \rho U \phi) + \frac{\partial}{\partial y}(\varepsilon \rho V \phi) \right) dv \cdot dt = \\ \int_{VC}^{t+\Delta t} \int_{t} \left( \frac{\partial}{\partial x} \left( \Gamma \frac{\partial}{\partial x}(\phi) \right) + \frac{\partial}{\partial y} \left( \Gamma \frac{\partial}{\partial y}(\phi) \right) + S_{\phi} \right) dv \cdot dt \quad (19)$$

avec $dv = dx \cdot dy$ l.

#### • Discrétisation des différents termes de l'équation ci-dessus

##### a- Terme de stockage

$$\int_{VC}^{t+\Delta t} \int_{t} \left( \frac{\partial}{\partial t}(\varepsilon \rho \phi) \right) dv \cdot dt = (\varepsilon \rho \phi)_{t}^{t+\Delta t} \cdot \Delta x \cdot \Delta y = ((\varepsilon \rho \phi_{p})^{t+\Delta t} - (\varepsilon \rho \phi_{p})^{t}) \cdot \Delta x \cdot \Delta y$$

##### b- Terme convectif longitudinal

$$\int_{t}^{t+\Delta t} \int_{VC} \left( \frac{\partial}{\partial t}(\varepsilon \rho U \phi) \right) dv \cdot dt = (\varepsilon \rho U \phi)_{w}^{e} \cdot \Delta y \cdot \Delta t = ((\varepsilon \rho U \phi)_{e} - (\varepsilon \rho U \phi)_{w}) \cdot \Delta y \cdot \Delta t = (F_{e} \phi_{e} - F_{w} \phi_{w}) \cdot \Delta t$$

avec: $F_e = \varepsilon \rho U_e \cdot \Delta y$; $F_w = \varepsilon \rho U_w \cdot \Delta y$.

##### c- Terme convectif vertical

$$\int_{t}^{t+\Delta t} \int_{VC} \left( \frac{\partial}{\partial t}(\varepsilon \rho V \phi) \right) dv \cdot dt = (\varepsilon \rho V \phi)_{s}^{n} \cdot \Delta x \cdot \Delta t = ((\varepsilon \rho V \phi)_{n} - (\varepsilon \rho V \phi)_{s}) \cdot \Delta x \cdot \Delta t = (F_{n} \phi_{n} - F_{s} \phi_{s}) \cdot \Delta t$$

avec: $F_n = \varepsilon \rho V_n \cdot \Delta x$; $F_s = \varepsilon \rho V_s \cdot \Delta x$.
### e- Terme diffusif suivant x

$$ 
\begin{aligned}
& \int_{t}^{t+\Delta t} \int_{VC} \left( \frac{\partial}{\partial x} \left( \Gamma \frac{\partial \Phi}{\partial x} \right) \right) \cdot dx \cdot dy \cdot dt = \left( \Gamma \frac{\partial \phi}{\partial x} \right)_{\text{w}}^{\text{e}} \cdot \Delta y \cdot \Delta t = \left( \Gamma \frac{\partial \phi}{\partial x} \Big|_{\text{e}} - \Gamma \frac{\partial \phi}{\partial x} \Big|_{\text{w}} \right) \cdot \Delta y \cdot \Delta t \\
&= \left( \Gamma_e \frac{\phi_E - \phi_P}{\delta x_{PE}} - \Gamma_w \frac{\phi_P - \phi_W}{\delta x_{WP}} \right) \cdot \Delta y \cdot \Delta t = \left( \frac{\Gamma_e \cdot \Delta y}{\delta x_{PE}} \cdot (\phi_E - \phi_P) - \frac{\Gamma_w \cdot \Delta y}{\delta x_{WP}} \cdot (\phi_P - \phi_W) \right) \cdot \Delta t \\
&= (D_e \cdot (\phi_E - \phi_P) - D_w \cdot (\phi_P - \phi_W)) \cdot \Delta t
\end{aligned}
$$

$$ \text{avec: } D_e = \frac{\Gamma_e \cdot \Delta y}{\delta x_{PE}}, $$

$$ D_w = \frac{\Gamma_w \cdot \Delta y}{\delta x_{WP}}. $$

### f- Terme diffusif suivant y

$$ 
\begin{aligned}
& \int_{t}^{t+\Delta t} \int_{VC} \left( \frac{\partial}{\partial y} \left( \Gamma \frac{\partial \Phi}{\partial y} \right) \right) \cdot dx \cdot dy \cdot dt = \left( \Gamma \frac{\partial \phi}{\partial x} \right)_{\text{s}}^{\text{n}} \cdot \Delta x \cdot \Delta t = \left( \Gamma \frac{\partial \phi}{\partial y} \Big|_{\text{n}} - \Gamma \frac{\partial \phi}{\partial y} \Big|_{\text{s}} \right) \cdot \Delta x \cdot \Delta t \\
&= \left( \Gamma_n \frac{\phi_N - \phi_P}{\delta y_{PN}} - \Gamma_s \frac{\phi_P - \phi_S}{\delta y_{SP}} \right) \cdot \Delta x \cdot \Delta t = \left( \frac{\Gamma_n \cdot \Delta x}{\delta x_{PN}} \cdot (\phi_N - \phi_P) - \frac{\Gamma_s \cdot \Delta x}{\delta y_{SP}} \cdot (\phi_P - \phi_S) \right) \cdot \Delta t \\
&= (D_n \cdot (\phi_N - \phi_P) - D_s \cdot (\phi_P - \phi_S)) \cdot \Delta t
\end{aligned}
$$

$$ \text{avec: } D_n = \frac{\Gamma_n \cdot \Delta x}{\delta y_{PN}}, $$

$$ D_s = \frac{\Gamma_s \cdot \Delta x}{\delta y_{SP}}. $$

### g- Terme source

$$ \int_{t}^{t+\Delta t} \int_{VC} (S_{\phi}) dv \cdot dt = S_{\phi} \cdot \Delta x \cdot \Delta y \cdot \Delta t = (S_{p} \cdot \phi_{p} + S_{u}) \cdot \Delta x \cdot \Delta y \cdot \Delta t $$

(On considère que le terme source $S_\phi$ est linaire par rapport à la grandeur $\phi$).

Après une substitution dans l'équation transport-diffusion et en appliquant un schéma implicite pour le temps, on obtient la forme discrétisée suivante:

$$ a_p^0 (\phi_p^{t+\Delta t} - \phi_p^t) + (F_e \phi_e - F_w \phi_w) + (F_n \phi_n - F_s \phi_s) = (D_e (\phi_E - \phi_P) - D_w (\phi_P - \phi_W)) \\ + (D_n (\phi_N - \phi_P) - D_s (\phi_P - \phi_S)) \cdot \Delta x \cdot \Delta y $$

$$ \text{avec : } a_p^0 = \varepsilon \rho \frac{\Delta x \cdot \Delta y}{\Delta t}. $$

Pour arriver à une forme finale, on doit choisir un schéma de discrétisation des termes convectifs, le schéma employé pour la discrétisation, c'est le schéma hybride. Ce schéma, dont la précision spatiale est réduite au premier ordre, possède les propriétés de transport et de conservation. Il se réduit à l'approximation du schéma centrée CSD lorsque $|P_{e_m}| < 2$, mais il devient identique au schéma Upwind UDS1 lorsque $|P_{e_m}| > 2$.

Après un réarrangement, on aura:

$$ a_p \cdot \phi_p^1 = a_w \cdot \phi_w^1 + a_e \cdot \phi_e^1 + a_s \cdot \phi_s^1 + a_n \cdot \phi_n^1 + a_p^0 \cdot \phi_p^0 + S_u \cdot \Delta x \cdot \Delta y $$

Les coefficients $a_m$ ($m = w, e, s, n$) sont donnés par:

$$ a_w = \max \left( F_w, D_w + \frac{F_w}{2}, 0 \right); a_e = \max \left( -F_e, D_e - \frac{F_e}{2}, 0 \right); $$
$$a_s = \max\left(-F_s, D_s - \frac{F_s}{2}, 0\right); a_n = \max\left(-F_n, D_n - \frac{F_n}{2}, 0\right)$$

Il est à noter que cette formulation est valable pour n'importe quelle localisation des interfaces entre les mailles. Elle permet de fournir des solutions physiquement acceptables aux problèmes sans toutefois garantir la précision de la solution obtenue.

## Equation de pression

Par la suite on va voir que l'équation de pression a une forme générale de l'équation de poisson:

$$\nabla(\Gamma \cdot \nabla(\phi)) + S_{\phi} = 0$$

L'intégration sur un volume s'écrit:

$$\int_{VC} \left( \frac{\partial}{\partial x} \left( \Gamma \frac{\partial}{\partial x}(\phi) \right) + \frac{\partial}{\partial y} \left( \Gamma \frac{\partial}{\partial y}(\phi) \right) + S_{\phi} \right) \cdot dv = 0$$

Sa forme discrétisée est la suivante:

$$a_p \cdot \phi_p = a_w \cdot \phi_w + a_e \cdot \phi_e + a_s \cdot \phi_s + a_n \cdot \phi_n + S_{\phi} \cdot \Delta x \cdot \Delta y$$

où les différents coefficients: $a_w = D_w$, $a_e = D_e$, $a_s = D_s$, $a_n = D_n$.

$$a_p = D_w + D_e + D_s + D_n - S_p \Delta x \cdot \Delta y$$

La difficulté principale, lors de la résolution des équations de Navier-Stokes est la vérification de l'équation de continuité à chaque pas de temps. En effet, en se basant sur les équations de quantité de mouvement discrétisées, celles-ci étant pleinement explicites, on pourra se contenter à chaque itération de déduire le champ de vitesse à partir du champ de pression.

Cependant, il y a fort à croire que le champ ainsi calculé ne respecte pas la condition de conservation de la masse. On serait d'autre part dans l'incapacité de déduire l'évolution du champ de pression à partir de la seule équation de continuité. En incompressible, le champ de pression n'est pas connu, et on a besoin des algorithmes de type couplage pression-vitesse. L'astuce numérique que nous allons mettre en œuvre pour pallier à cette difficulté repose sur une méthode basée sur un fractionnement du pas de temps. Réécrivons l'équation de conservation de la quantité de mouvement comme suit:

$$\left\{ \begin{aligned} \frac{\partial}{\partial t}(\varepsilon \rho U) &= -\frac{\partial(\varepsilon \rho U U)}{\partial x} - \frac{\partial(\varepsilon \rho V U)}{\partial y} + \frac{\partial}{\partial x}\left(\Gamma \frac{\partial U}{\partial x}\right) + \frac{\partial}{\partial y}\left(\Gamma \frac{\partial U}{\partial y}\right) + S_U - \varepsilon \frac{\partial P}{\partial x} \\ \frac{\partial}{\partial t}(\varepsilon \rho V) &= -\frac{\partial(\varepsilon \rho U V)}{\partial x} - \frac{\partial(\varepsilon \rho V V)}{\partial y} + \frac{\partial}{\partial x}\left(\Gamma \frac{\partial V}{\partial x}\right) + \frac{\partial}{\partial y}\left(\Gamma \frac{\partial V}{\partial y}\right) + S_V - \varepsilon \frac{\partial P}{\partial y} \end{aligned} \right.$$

On décompose le terme de stockage $\left(\frac{\partial}{\partial t}\right)$ du temps $t$ à $t+\frac{\Delta t}{2}$ et $t+\frac{\Delta t}{2}$ à $t+\Delta t$, on aura:
$$
\left\{
\begin{aligned}
& \frac{\epsilon \rho U^{t+\Delta t/2} - \epsilon \rho U^t}{\Delta t} + \frac{\epsilon \rho U^{t+\Delta t} - \epsilon \rho U^{t+\Delta t/2}}{\Delta t} = -\frac{\partial(\epsilon \rho U U)}{\partial x} - \frac{\partial(\epsilon \rho V U)}{\partial y} + \\
& \qquad \frac{\partial}{\partial x} \left( \Gamma \frac{\partial U}{\partial x} \right) + \frac{\partial}{\partial y} \left( \Gamma \frac{\partial U}{\partial y} \right) + S_U - \epsilon \frac{\partial P}{\partial x} \\[2ex]
& \frac{\epsilon \rho V^{t+\Delta t/2} - \epsilon \rho V^t}{\Delta t} + \frac{\epsilon \rho V^{t+\Delta t} - \epsilon \rho V^{t+\Delta t/2}}{\Delta t} = -\frac{\partial(\epsilon \rho U V)}{\partial x} - \frac{\partial(\epsilon \rho V V)}{\partial y} + \\
& \qquad \frac{\partial}{\partial x} \left( \Gamma \frac{\partial V}{\partial x} \right) + \frac{\partial}{\partial y} \left( \Gamma \frac{\partial V}{\partial y} \right) + S_V - \epsilon \frac{\partial P}{\partial y}
\end{aligned}
\right.
$$

On pose:

$$
\left\{
\begin{aligned}
\frac{\epsilon \rho U^{t+\Delta t/2} - \epsilon \rho U^t}{\Delta t} &= -\frac{\partial(\epsilon \rho U U)}{\partial x} - \frac{\partial(\epsilon \rho V U)}{\partial y} + \frac{\partial}{\partial x} \left( \Gamma \frac{\partial U}{\partial x} \right) + \frac{\partial}{\partial y} \left( \Gamma \frac{\partial U}{\partial y} \right) + S_U \\
\frac{\epsilon \rho V^{t+\Delta t/2} - \epsilon \rho V^t}{\Delta t} &= -\frac{\partial(\epsilon \rho U V)}{\partial x} - \frac{\partial(\epsilon \rho V V)}{\partial y} + \frac{\partial}{\partial x} \left( \Gamma \frac{\partial V}{\partial x} \right) + \frac{\partial}{\partial y} \left( \Gamma \frac{\partial V}{\partial y} \right) + S_V \\
\frac{\epsilon \rho U^{t+\Delta t} - \epsilon \rho U^{t+\Delta t/2}}{\Delta t} &= -\epsilon \frac{\partial P}{\partial x} && (a) \\
\frac{\epsilon \rho V^{t+\Delta t} - \epsilon \rho V^{t+\Delta t/2}}{\Delta t} &= -\epsilon \frac{\partial P}{\partial y} && (b)
\end{aligned}
\right.
$$

La résolution des deux dernières équations (a) et (b) représente l'étape dite de prédiction, elle
permet, l'obtention d'un champ provisoire de vitesse $U^{t+\frac{1}{2}\Delta t}$ et $V^{t+\frac{1}{2}\Delta t}$ calculé uniquement à
partir du champ $U^t$ et $V^t$. Ce champ provisoire ne vérifie a priori pas l'équation de continuité.
On dérive les deux dernières équations respectivement par rapport à $x$ et par rapport à $y$, d'où:

$$
\frac{\partial}{\partial x} \left( \frac{\rho U^{t+\Delta T} - \rho U^{t+\frac{1}{2}\Delta t}}{\Delta t} \right) = - \frac{\partial}{\partial x} \left( \frac{\partial P}{\partial x} \right) , \quad \frac{\partial}{\partial y} \left( \frac{\rho V^{t+\Delta T} - \rho V^{t+\frac{1}{2}\Delta t}}{\Delta t} \right) = - \frac{\partial}{\partial y} \left( \frac{\partial P}{\partial y} \right)
$$

On fait la somme et on obtient :

$$
\frac{\partial}{\partial x} \left( \frac{\rho U^{t+\Delta T} - \rho U^{t+\frac{1}{2}\Delta t}}{\Delta t} \right) + \frac{\partial}{\partial y} \left( \frac{\rho V^{t+\Delta T} - \rho V^{t+\frac{1}{2}\Delta t}}{\Delta t} \right) = - \frac{\partial}{\partial x} \left( \frac{\partial P}{\partial x} \right) - \frac{\partial}{\partial y} \left( \frac{\partial P}{\partial y} \right) \\
\left[ \left( \frac{\partial(\rho U^{t+\Delta T})}{\partial x} + \frac{\partial(\rho V^{t+\Delta T})}{\partial y} \right) - \left( \frac{\partial(\rho U^{t+\frac{1}{2}\Delta t})}{\partial x} + \frac{\partial(\rho V^{t+\frac{1}{2}\Delta t})}{\partial y} \right) \right] / \Delta t = - \frac{\partial^2 P}{\partial x^2} - \frac{\partial^2 P}{\partial y^2}
$$
$$ \left( \frac{\partial (\rho U^{t+\Delta t})}{\partial x} + \frac{\partial (\rho V^{t+\Delta t})}{\partial y} \right) = 0 \quad (\text{équation de continuité}), $$

$$ \text{on aura: } \frac{\partial^2 P}{\partial x^2} - \frac{\partial^2 P}{\partial y^2} = \left( \frac{\partial \left( \rho U^{t+\frac{1}{2}\Delta t} \right)}{\partial x} + \frac{\partial \left( \rho V^{t+\frac{1}{2}\Delta t} \right)}{\partial y} \right) / \Delta t $$

La résolution de l'équation de pression (ci-dessus), nécessite évidemment de fixer des conditions aux limites. On choisira une condition de type Neumann sur le gradient de pression au niveau des frontières du domaine $\left(\frac{\partial P}{\partial n} = 0\right)$.

Après convergence, on recalcule les pressions aux frontières du domaine, en les extrapolant à partir du champ en aval. La dernière étape de la méthode consiste évidemment à réinjecter le champ de pression, calculé ci-dessus, dans les équations (a) et (b), ce qui permet d'obtenir le champ de vitesse $U^{t+\Delta t}$, $V^{t+\Delta t}$ celui-ci satisfaisant aussi bien à l'équation de quantité de mouvement qu'à l'équation de continuité:

$$ U_{IJ}^{t+\Delta t} = U_{IJ}^{t+\Delta t/2} - \Delta t \cdot \frac{(P_{I+1,J} - P_{IJ})}{(\rho \cdot \Delta x)} $$

$$ V_{IJ}^{t+\Delta t} = V_{IJ}^{t+\Delta t/2} - \Delta t \cdot \frac{(P_{I,J+1} - P_{IJ})}{(\rho \cdot \Delta y)} $$

Il existe d'autres algorithmes qui résolvent le problème de couplage pression-vitesse comme Simple, Simplec, etc...

Dans ce qui suit, on s'intéresse à la résolution du modèle mathématique évoqué, qui est relatif aux différentes équations de transport dans les parties de la pile à combustible du côte cathodique (canal, couche de diffusion, couche de catalyseur).

Ces équations sont:

- Equation de mouvement pour déterminer les composantes de vitesses U, V.

- Equation de l'énergie pour déterminer la température T.

- Deux équations de transport pour la fermeture du modèle de turbulence: l'énergie cinétique et la dissipation.

Ces équations peuvent s'écrire sous la forme générale suivante:

$$ \frac{\partial}{\partial t}(\varepsilon \rho \phi) + \frac{\partial}{\partial x}(\varepsilon \rho U \phi) + \frac{\partial}{\partial y}(\varepsilon \rho V \phi) = \frac{\partial}{\partial x}\left(\Gamma \frac{\partial}{\partial x}(\phi)\right) + \frac{\partial}{\partial y}\left(\Gamma \frac{\partial}{\partial y}(\phi)\right) + S_{\phi} $$

On retrouve les équations de chaque composant de la pile et ceci en affectant les quantités appropriées au grandeur $\phi$, au coefficient de diffusion $\Gamma$ et au terme source $S_{\phi}$ (voir tableaux ci-dessous). La porosité du canal, couche de diffusion et couche de catalyseur est égale respectivement à 1, $\varepsilon_d$, $\varepsilon_c$.

**Equation de mouvement**, $\phi = U_i$ (selon $x$ $U_i = U$, $x_i = x$, selon $y$ $U_i = V$, $x_i = y$)
<table><thead><tr><th></th><th>Γ</th><th>S<sub>φ</sub></th></tr></thead><tbody><tr><td>Canal</td><td>μ + μ<sub>t</sub></td><td>- ∂P/∂x<sub>i</sub></td></tr><tr><td>Couche de diffusion</td><td>(ε<sub>d</sub><sup>1.5</sup>) . μ</td><td>- μ/k<sub>d</sub> . ε<sub>d</sub><sup>2</sup> . U<sub>i</sub> - ε<sub>d</sub> . ∂P/∂x<sub>i</sub></td></tr><tr><td>Couche de catalyseur</td><td>(ε<sub>c</sub><sup>1.5</sup>) . μ</td><td>- μ/k<sub>c</sub> . ε<sub>c</sub><sup>2</sup> . U<sub>i</sub> - ε<sub>c</sub> . ∂P/∂x<sub>i</sub></td></tr></tbody></table>

Equation d'énergie, $\phi = T$

<table><thead><tr><th></th><th>Γ</th><th>S<sub>φ</sub></th></tr></thead><tbody><tr><td>Canal</td><td>α + μ<sub>t</sub>/Pr</td><td>0</td></tr><tr><td>Couche de diffusion</td><td>(λ<sub>d</sub>/ρ·C<sub>p</sub>)<sup>eff</sup></td><td>0</td></tr><tr><td>Couche de catalyseur</td><td>(λ<sub>c</sub>/ρ·C<sub>p</sub>)<sup>eff</sup></td><td>i<sup>2</sup>/σ<sub>c</sub> + i·η</td></tr></tbody></table>

Equation d'énergie cinétique ($\phi = k$) et de dissipation ($\phi = \varepsilon$)

<table><thead><tr><th></th><th>Γ</th><th>S<sub>φ</sub> (cinétique)</th><th>S<sub>φ</sub> (dissipation)</th></tr></thead><tbody><tr><td>Canal</td><td>μ + μ<sub>t</sub></td><td>G - ε</td><td>(c<sub>1ε</sub>·G - c<sub>2ε</sub>·ε)·ε/k</td></tr><tr><td>Couche de diffusion</td><td>(ε<sub>d</sub><sup>1.5</sup>) . μ</td><td>0</td><td>0</td></tr><tr><td>Couche de catalyseur</td><td>(ε<sub>c</sub><sup>1.5</sup>) . μ</td><td>0</td><td>0</td></tr></tbody></table>

$$ \text{avec: } G = \mu_t \left( \frac{u_{ij+1} - u_{ij-1}}{2 \Delta y} \right)^2 . $$

Equation de pression, $\phi = P$, $\Gamma = 1$

$$ S_{\phi} = - \left( \frac{\partial(\rho.U)}{\partial x} + \frac{\partial(\rho.V)}{\partial y} \right) / \Delta t $$

## 4 RESULTATS ET DISCUSSIONS

Un programme informatique (Fortran) est développé pour l'obtention des profils de température, et le champ de vitesse. Parmi les paramètres importants qui ont une grande influence sur le fonctionnement de la pile et qui ont été analysés sont: la densité de courant (ou la surtension cathodique), la perméabilité de la couche de diffusion (GDL) et de la couche de catalyseur. Le régime d'écoulement est considéré turbulent dans les canaux des gaz et laminaire dans les autres composants. Le modèle de turbulence étudié dans notre cas et celui de $k_\epsilon$. Les résultats illustrent certains phénomènes de transport qui se produisent dans la pile à combustible et démontrent les capacités du modèle mathématique utilisé qui répondent aux différentes conditions aux limites. Le comportement transitoire de la pile à combustible PEMFC est aussi discuté.
## a- Champ de vitesse

Dans cette étape, on s'intéresse à étudier le champ dynamique (profil de vitesse) pour différents paramètres, tels que le régime d'écoulement (influence du nombre de Reynolds) et la nature du milieu continu ou discontinu (perméabilité). La fermeture du modèle de turbulence est testée sur le modèle de longueur de mélange et le modèle $k-\epsilon$.

La figure 2 (a)-(b) montre une comparaison entre les deux modèles de turbulence pour un nombre de Reynolds $Re = 10^4$ et pour deux valeurs différentes de perméabilité de la couche de diffusion. Pour une valeur de perméabilité de $10^{-15}$, la différence entre les deux modèles est nette, surtout à l'interface de la couche de diffusion, car le modèle de longueur de mélange demande des informations sur l'épaisseur de la couche limite pour calculer la viscosité turbulente. Les profils de la vitesse représentés dans cette figure sont ceux de la composante horizontale. A cet effet, le modèle $k-\epsilon$ est opté, ce modèle demande aussi des informations sur les conditions d'entrée du canal, surtout la valeur de l'intensité turbulente.

<figure><img src="image_2.png" /><figcaption>Fig. 2: Profil de vitesse dans le canal, GDL et couche de catalyseur</figcaption></figure>

La figure 3 représente l'effet de la perméabilité $k$ de la couche de diffusion sur le comportement de l'écoulement, pour différentes valeurs du nombre de Reynolds. On remarque que le profil de vitesse a une forme parabolique symétrique pour les faibles valeurs de perméabilités (cas d'écoulement entre deux parois solides). La valeur de la vitesse augmente lorsque la perméabilité est grande, et devient plus grande en augmentant le nombre de Reynolds. On distingue que pour une perméabilité très grande donc le terme de Darcy devient négligeable, et pour une perméabilité se rapprochant de 0 le terme de Darcy est prépondérant. La représentation en deux dimensions des champs de vitesses se trouve dans la figure 4 (a)-(b) Ces contours confirment ce qui a été dit précédemment concernant l'effet de la perméabilité.

<figure><img src="image_3.png" /><figcaption>Fig. 3: Effet du nombre de Reynolds et de la perméabilité sur les profils de vitesse (modèle $k-\epsilon$)</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 4: Répartition de champs de vitesses</figcaption></figure>

b- Champ de température

L'effet de la réaction électrochimique qui peut être représentée par la surtension cathodique η est analysé. Le régime d'écoulement turbulent pour un nombre de Reynolds de Re = 10⁴ est choisi. Les résultats prouvent et montrent qu'à la cathode la réaction électrochimique fait augmenter la température. Les champs de température sont présentés sous forme de contours dans la figure 5. L'influence de l'écoulement et de la réaction électrochimique η sur la répartition de ces paramètres dans les trois parties cathodique (canal, couche de diffusion (GDL), couche de catalyseur) dans le plan OXY est analysée. La description des profils montre qu'il existe deux zones importantes: la première zone représente la couche limite massique et thermique à la paroi du canal où on observe l'augmentation de l'épaisseur de cette couche en se déplaçant le long du canal. On peut voir dans la deuxième zone, l'effet de la surtension d'activation quand on atteint la zone de réaction électrochimique (couche de catalyseur). Pour deux valeurs de la surtension d'activation (η = 0 et η = 0.3), on remarque que la diffusion de masse et de la chaleur dans le milieu poreux est ralentie à cause de l'effet de la perméabilité de ce milieu, dans le cas où η = 0 (surtensions nulles). Par contre, si on tient compte de la réaction, la répartition de température est influencée.

<figure><img src="image_5.png" /><figcaption>Fig. 5: Distribution de température dans le canal, GDL et couche de catalyseur
Pour différentes valeurs de surtension d'activation</figcaption></figure>
### c- Evolution transitoire de température

La figure 6 représente, le comportement transitoire de la température pour différentes valeurs de Perméabilité et de la surtension d'activation. On remarque une augmentation de l'évolution temporelle par rapport à l'état initial.

<figure><img src="image_6.png" /><figcaption>Fig. 6: Effet de la perméabilité et de la surtension d'activation sur la température</figcaption></figure>

### d-Evolution dans le Plan OYZ (P2)

#### -Contour de température

Pour plus d'éclaircissement des phénomènes, une représentation des résultats dans le plan P2 (OYZ) s'impose. L'effet des mêmes paramètres, réaction et perméabilité a été présenté sous forme de contours dans le coté cathodique; couche du catalyseur et couche de diffusion. La perméabilité et la surtension d'activation agissent de la même manière sur la distribution de température. La visibilité de l'effet de la perméabilité est claire dans ce plan contrairement au plan P1 (OXY). Pour une température fixe à l'entrée du canal (température d'oxygène humide), la diffusion de la chaleur est influencée par la perméabilité avec ou sans réaction.

Dans le cas d'une perméabilité importante et sans réaction tous les composants sont presque à la même température (Fig. 7.b), la chaleur transportée par les gaz entrant occupe environ la totalité de la cathode. Cette chaleur se localise aux entrées des canaux. Pour des faibles perméabilités (Fig. 7.a), cette chaleur est limitée aux entrées des gaz dans la GDL (la diffusion de la chaleur décroît du canal vers la couche catalyseur).

Par contre avec réaction la propagation de la chaleur se fait dans le sens opposé. Elle est plus importante entre les deux canaux et ceci s'explique par les flux de gaz entrant qui sont froid. Pour une faible perméabilité la cathode s'échauffe presque dans sa totalité et les deux entrées sont moins chaudes (Fig. 7.c).

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7: Distribution de température dans la couche de diffusion et de la couche de catalyseur pour différentes perméabilités et surtensions d'activation</figcaption></figure>

## 5. CONCLUSION

Notre travail est basé sur l'étude de la performance de la pile à combustible à membrane échangeuse de protons. Pour une meilleure compréhension des phénomènes ayant lieu dans une PEMFC, des équations mathématiques ainsi que des relations empiriques ont été développées. L'outil numérique choisi pour la résolution de ces équations est la méthode des volumes finis. Des arrangements d'approximation numériques sont abordés pour atteindre la solution du problème étudié. Des difficultés sont rencontrées en raison de la présence de dimensions spatiales extrêmement petites aussi bien que des coefficients variables qui ont exigé beaucoup de changements et d'approximations numériques.

D'après les analyses des résultats obtenus, il est apparu que le modèle développé pour la pile peut être un outil de conception futur très puissant (milieu poreux, écoulement turbulent, réactions chimiques, présence des termes sources, etc....).

Les solutions du modèle ont éclairci beaucoup de phénomènes dans les composants de la pile: canal, couche de diffusion, couche de catalyseur et ont permis des analyses plus détaillées de chacun.

Les effets suivants ont été observés :

- Le champ de vitesses est obtenu pour un écoulement turbulent dans un domaine de dimensions faibles. La comparaison des modèles de turbulence nous a permis de faire le choix adéquat.

- Une bonne compréhension de la différence entre l'écoulement dans un milieu continu et un milieu poreux, avec et sans réaction est acquise.

- La température et la pression d'entrée du canal sont des paramètres d'une grande influence sur la performance de la pile.

- L'obtention de l'évolution des grandeurs physiques dans deux plans différents (transversal et longitudinal) permet de bien cerner les phénomènes.
<figure><img src="image_9.png" /><figcaption>Fig. 8: Organigramme de calcul</figcaption></figure>

## NOMENCLATURE

<table><tr><td>a</td><td>Activité de l'eau</td><td>M<sub>H+</sub></td><td>Masse molaire du proton d'hydrogène, (kg/mole)</td></tr><tr><td>a. J<sub>0,c</sub><sup>ref</sup></td><td>Densité de courant volumique de référence, (A/m<sup>3</sup>)</td><td>P</td><td>Pression, (N/m<sup>2</sup>)</td></tr><tr><td>C<sub>O2</sub></td><td>Concentration d'oxygène, (mole/m<sup>3</sup>)</td><td>Pe</td><td>Nombre de Peclet</td></tr><tr><td>C<sub>H2O</sub></td><td>Concentration de l'eau, (mole/m<sup>3</sup>)</td><td>Pr</td><td>Nombre de Prandtl</td></tr><tr><td>C<sub>H+</sub></td><td>Concentration du proton d'hydrogène, (mole/m<sup>3</sup>)</td><td>Re</td><td>Nombre de Reynolds</td></tr><tr><td>C<sub>pH2O</sub></td><td>Chaleur spécifique de l'eau, (J/kg.K)</td><td>R<sub>gaz</sub></td><td>Constante des gaz, (J/mole.K)</td></tr><tr><td>C<sub>pH+</sub></td><td>Chaleur spécifique du proton d'hydrogène, (J/kg.K)</td><td>r</td><td>Rapport de correction du milieu poreux</td></tr><tr><td>D</td><td>Coefficient de diffusion, (m<sup>2</sup>/s)</td><td>S<sub>U</sub></td><td>Terme source de l'équation de mouvement</td></tr><tr><td>D<sub>k</sub></td><td>Coefficient de diffusion de l'espèce k, (m<sup>2</sup>/s)</td><td>Sc</td><td>Nombre de Schmidt</td></tr></table>
<table><tr><td>D<sub>k</sub><sup>eff</sup></td><td>Coefficient de diffusion effectif de l'espèce k, (m<sup>2</sup>/s)</td><td>T</td><td>Température, (K)</td></tr><tr><td>F</td><td>Constante de Faraday, (C/mole)</td><td>t</td><td>Temps, (s)</td></tr><tr><td>h<sub>c</sub></td><td>Hauteur du canal, (m)</td><td>U</td><td>Composante longitudinale de la vitesse, (m/s)</td></tr><tr><td>h<sub>d</sub></td><td>Hauteur de couche de diffusion, (m)</td><td>V</td><td>Composante transversale de la vitesse, (m/s)</td></tr><tr><td>h<sub>ct</sub></td><td>Hauteur de couche de catalyseur, (m)</td><td>x</td><td>Abscisse, (m)</td></tr><tr><td>h<sub>m</sub></td><td>Hauteur de la membrane, (m)</td><td>y</td><td>Ordonnée, (m)</td></tr><tr><td>K</td><td>Courant surfacique dans la couche de catalyseur, (A/m<sup>2</sup>)</td><td>ε</td><td>Dissipation de l'énergie cinétique ou porosité</td></tr><tr><td>J<sub>c</sub></td><td>Courant volumique cathodique, (A/m<sup>3</sup>)</td><td>ε<sub>d</sub></td><td>Porosité de couche de diffusion</td></tr><tr><td>K<sub>d</sub></td><td>Perméabilité de la couche de diffusion, (m<sup>2</sup>)</td><td>ε<sub>c</sub></td><td>Porosité de couche du catalyseur</td></tr><tr><td>K<sub>c</sub></td><td>Perméabilité de la couche de catalyseur, (m<sup>2</sup>)</td><td>ρ</td><td>Masse volumique, (kg/m<sup>3</sup>)</td></tr><tr><td>k</td><td>Energie cinétique turbulente, (m<sup>2</sup>/s<sup>2</sup>)</td><td>σ<sub>c</sub></td><td>Conductivité électrique du catalyseur, (Ω<sup>-1</sup>.m<sup>-1</sup>)</td></tr><tr><td>k<sup>eff</sup></td><td>Conductivité thermique effective, (W/m.K)</td><td>η</td><td>Surtension d'activation de la cathode, (V)</td></tr><tr><td>L<sub>x</sub></td><td>Longueur de la pile suivant x, (m)</td><td>μ</td><td>Viscosité dynamique, (kg/m.s)</td></tr><tr><td>L<sub>z</sub></td><td>Longueur de la pile suivant z, (m)</td><td>μ<sub>H+</sub></td><td>Viscosité dynamique du proton de l'hydrogène, (kg/m.s)</td></tr><tr><td>l<sub>m</sub></td><td>Longueur de mélange, (m)</td><td>μ<sub>t</sub></td><td>Viscosité turbulente, (kg/m.s)</td></tr><tr><td>M<sub>H2O</sub></td><td>Masse molaire de l'eau, (kg/mole)</td><td>μ<sub>H2O</sub></td><td>Viscosité dynamique de l'eau, (kg/m.s)</td></tr><tr><td>M<sub>O2</sub></td><td>Masse molaire de l'oxygène, (kg/mole)</td><td>MEA</td><td>Assemblage membrane électrode</td></tr></table>

REFERENCES

[1] B. Eaton, M.R. Von Spakovsky, M.W. Ellis, D.J. Nelson, B. Olsommer and N. Siegel, ‘One-Dimensional, Transient Model of Heat, Mass, and Charge Transfer in a Proton Exchange Membrane’, Proceedings of IMECE, 2001.

[2] S. Um, C.Y. Wang and K.S. Chen, ‘Computational Fluid Dynamics Modeling of Proton Exchange Membrane Fuel Cells’, J. Electrochemical Soc., Vol. 147, N°12, pp. 485 - 493, 2000.

[3] Z.H. Wang and C.Y. Wang, ‘Two-Phase Flow and Transport in the Interdigitated Air Cathode of Proton Exchange Membrane Fuel Cells’, Proceedings ASME Heat Transfer Division, HTD-Vol. 366-1, pp. 27 - 33, 2000.

[4] D. Natarajan and T. Van Nguyen, ‘A Two-Dimensional, Two-Phase, Multicomponent, Transient Model of the Cathode of a Proton Exchange Membrane Fuel Cell Using Conventional Gas Distributors’, J. Electrochemical Soc., Vol. 148, N°12, pp. A1324 - A1335, 2001.

[5] D. Genevey, M.R. Von Spakovsky, M.W. Ellis, D.J. Nelson, B. Olsommer, F. Topin, N. Montel and N.P. Siegel, ‘Transient Model of the Heat, Mass and Charge Transfer as Well as Electrochemistry in the Catalyst Layer of a PEMFC’, International Mechanical Engineering Congress and Exposition – IMECE’2002, ASME IMECE Paper N°33322, N.Y., November 2002.

[6] D.M. Bernardi and M.W. Verbrugge, ‘A Mathematical Model of the Solid-Polymer-Electrolyte Fuel Cell’, J. Electrochemical Soc., Vol. 139, N°9, pp. 2477 - 2490, 1992.

[7] J.S. Yi and T.V. Nguyen, ‘An Along-the-Channel Model for Proton Exchange Membrane Fuel Cells’, J. Electrochemical Soc., Vol. 145 (4), pp. 1149 - 1159, 1998.

[8] N.P. Siegel, M.W. Ellis, D.J. Nelson and M.R. Von Spakovsky, ‘Single Domain PEMFC Model Based on Agglomerate Catalyst Geometry’, J. Power Sources, Vol. 115, pp. 81 - 89, 2003.

[9] K. Oulmi, H. Ben Moussa and D. Haddad, ‘One-Dimensional, Transient Numerical Study of Mass, Heat and Charge Transfer in a Proton Exchange Membrane for PEMFC’, Electrochemical Energy
Conversion and Storage (Fuel Cells), The 56th International Society of Electrochemistry, Bexco, Busan, Korea, Sept. 25-30, 2005.

[10] H. Ben Moussa, D. Haddad et K. Oulmi. 'Simulation Numérique par Eléments Finis de Transfert Thermique, d'Humidité et de Charges dans une Membrane Echangeuse de Proton d'une Pile à Combustible', 12èmes Journées Internationales de Thermique (JITH 2005), pp. 85, Tanger, Maroc, 15-17 Nov. 2005.

[11] T.E. Springer and S. Gottesfeld, 'Pseudohomogeneous Catalyst Layer for Polymer Electrolyte Fuel Cells', Proceedings of the Symposium on Modelling of Batteries and Fuel Cells, R.E. White, M.W. Verbrugge and J.F. Stockel, Editors, PV 91-10, pp. 197 - 208, The Electrochemical Society Softbound Proceedings Series, Pennington, N.J. (b), 1991.

[12] T.E. Springer, M.S. Wilson and S. Gottesfeld, 'Modelling and Experimental Diagnostics in Polymer Electrolyte Fuel Cells', J. Electrochemical Soc., Vol 140, N°12, pp. 3513-3526, 1993.

[13] J. Bard et E Faulkner, 'Electrochimie, Principes, Méthodes et Applications', Ed. Masson, 1980. Peter W. Atkins, Chimie Physique, Bock Université, Chap. 24, pp. 723-753, 2000.

[14] J. Cousteix, 'Turbulence et Couche Limite', Cepadues-Editions, CEPAD, 1989.

[15] D.B. Genevey, 'Transient Model of Heat, Mass, and Charge Transfer as Well as Electrochemistry in the Cathode Catalyst Layer of a PEMFC', Thesis Master of Science in Mechanical Engineering, Blacksburg, Virginia, 2001.