Repositorio Institucional de Documentos

Obtenci´on por simulaciones num´ericas de la curva de polarizaci´on de una pila de combustible de membrana de intercambio de protones mediante software libre Trabajo Fin de M´aster presentado en la Universidad de Zaragoza para la obtenci´on del grado de M´aster en Mec´anica Aplicada por Leonard Efr´en Due˜ nas Guti´ errez Director Dr. Radu MUSTATA CoDirector Dr. Luis VALI ˜ NO POP en Ingen´ıeria Mec´anica y de Materiales Curso acad´emico 2010-2011 Centro Polit´ecnico Superior Universidad de Zaragoza Septiembre, 2011 En esta pagina quiero dedicar este trabajo a dos seres muy especiales. Primer que todo a mi abuelita Cruz Ana Moreno (cachanita) la cual llevo en mi coraz´on como un sello imborrable, y su recuerdo sigue tan intacto como si nunca hubiese partido. No podr´ıa dejar pasar y mencionar su nombre, porque cada paso que doy siempre va presente su recuerdo, por eso como lo dije un d´ıa no vean la vida que dejo, sino que miren lo que comenc´e y por donde sigo avanzando, el camino es largo pero llegar´e a lo mas alto. Y en segundo lugar al que quiero dedicar este trabajo al ser que me tiene en este punto de mi proceso, mi Dios creador, ya que siempre me ha concedido todo lo que le he pedido y d´ıa a d´ıa me protege y bendice en todos lo sentidos en mi vida. Nunca me cansare de agradecer a estos dos seres ....... los amo con todo mi Coraz´on Agradecimientos En el momento de realizar esta trabajo quiero extender un agradecimiento especial al Dr. Radu Mustata (director de tesis) y tambi´en al Dr. Luis Vali˜no (codirector) por su gran colaboraci´on y por ser una invaluable fuente de conocimiento y apoyo incondicional. Quiero agradecer, como siempre, el apoyo y los ´animos constantes de mi familia, en especial a mi madre Mar´ıa Roc´ıo Guti´errez, la cual d´ıa a d´ıa me brindaba su apoyo y ´animo para seguir adelante, y tambi´en a mi hermana Lina Rocio Bravo por su apoyo incondicional. Y por ´ultimo quiero agradecer a mis mejores amigos este trabajo, ya que su apoyo para seguir adelante fue de mucha ayuda. Gracias a todos. Resumen En ese trabajo se propone obtener num´ericamente la curva de polarizaci´on de una pila de combustible PEM (proton exchange membrane) mediante una discretizaci´on por vol´umenes finitos, utilizando software de dominio p´ublico, en concreto OpenFOAM. Se considera una geometr´ıa sencilla, pero con todos los elementos de la pila (no se resuelve todo el ”stack”, sino solo una monocelda). Para este trabajo se realizan algunas simplificaciones como las siguientes: 1) Geometr´ıa 2D simplificada, pero se incluyen todos los componentes de la pila. 2) Caso isotermo. 3) Caso estacionario. 4) Densidad constante. 5) Difusi´on de Fick. 6) PEM de alta temperatura (no hay agua l´ıquida). 7) Hidr´ogeno y vapor de agua en la entrada an´odica y ox´ıgeno y vapor de agua en la cat´odica. 8) Se consideran p´erdidas por activaci´on, masa y ´ohmicas en la membrana. 9) Ecuaciones de transporte completas, incluyendo medio poroso. Este trabajo requiere primero la correcta modelizaci´on de todos los fen´omenos f´ısicos que tienen lugar en la pila (fluidodin´amica, flujo en medios porosos, incapacidad de transporte de hidr´ogeno y ox´ıgeno, reacciones electroqu´ımicas, transporte de cargas el´ectricas) y que est´an fuertemente acoplados. Una vez obtenidas de la literatura las ecuaciones de transporte adecuadas al problema, con las contribuciones personales necesarias y el modelo de acoplo, se proceder´a a su discretizaci´on por vol´umenes finitos en una malla sencilla bidimensional que representa una monocelda. Los par´ametros f´ısicos que caracterizan los distintos materiales y componentes de la pila se obtendr´an de los valores proporcionados por los fabricantes, y que est´an a disposicion del LITEC-CSIC. Una vez obtenido un c´odigo fiable, se obtendr´a la curva de polarizaci´on de la pila mediante la variaci´on en la demanda de corriente el´ectrica a la monocelda, que se impone como una condici´on de contorno al problema. ´ Indice general ´ Indice general V ´ Indice de figuras VII ´ Indice de cuadros IX 1. Introducci´on 1 1.0.1. Curva Polarizaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.0.2. Consideraciones te´oricas del modelado de los fen´omenos de transporte en el interior de una pila PEM . . . . . . . . . . . . . . . . 3 1.0.3. Entornos de simulaci´on . . . . . . . . . . . . . . . . . . . . . . . 7 1.0.4. Objetivos del trabajo . . . . . . . . . . . . . . . . . . . . . . . . 9 2. Modelo Matem´atico para un pila de tipo PEM 11 2.1. Ecuaciones que describen el movimiento fluido. . . . . . . . . . . . . . . 12 2.1.1. Ecuaciones para el medio poroso . . . . . . . . . . . . . . . . . . 13 2.1.2. ´ Anodo................................. 14 2.1.3. C´atodo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2. Ecuaciones Campo El´ectrico, Membrana . . . . . . . . . . . . . . . . . . 16 2.2.1. Condiciones de contorno . . . . . . . . . . . . . . . . . . . . . . . 18 2.3. Acoplamiento de las ecuaciones y m´etodo num´erico . . . . . . . . . . . . 19 2.3.1. M´etodo iterativo caso real simplificado . . . . . . . . . . . . . . . 20 2.3.2. Comentarios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3. Simulaci´on Num´erica 23 3.1. Geometr´ıa del Dominio . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2. Malla computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3. Simulaci´on num´erica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 v ´ INDICE GENERAL 3.3.1. Par´ametros de simulaci´on . . . . . . . . . . . . . . . . . . . . . . 25 4. Resultados Num´ericos 27 4.1. Velocidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1.1. Hidr´ogeno y oxigeno . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2. Curva de polarizaci´on y otras . . . . . . . . . . . . . . . . . . . . . . . . 33 5. Conclusiones 39 Bibliograf´ıa 41 vi ´ Indice de figuras 1.1. Esquema de una pila de combustible de tipo PEM. . . . . . . . . . . . . 1 1.2. Curva de Polarizaci´on de una pila de combustible de membrana de intercambio de protones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.1. Geometr´ıa del dominio. . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2. Malla del domino computacional. . . . . . . . . . . . . . . . . . . . . . . 24 4.1. Modulo de la Velocidad en el ´anodo(izquierdo) y c´atodo(derecho) . . . . 28 4.2. Perfil de velocidad en corte transversal del ´andodo a media altura. . . . 29 4.3. Vectores de velocidad ´anodo y c´atodo. . . . . . . . . . . . . . . . . . . . 30 4.4. Perfil de velocidad salida capa catal´ıtica del ´anodo. . . . . . . . . . . . . 30 4.5. Fracciones m´asicas ´anodo y c´atodo . . . . . . . . . . . . . . . . . . . . . 31 4.6. Variaci´on de la fracci´on m´asica H2yO2.................. 32 4.7. Variaci´on de la densidad de corriente en las capas catal´ıticas. . . . . . . 32 4.8. Perfiles de η´anodo y c´atodo. . . . . . . . . . . . . . . . . . . . . . . . . 33 4.9. Comparativo curva de p´erdidas por activaci´on. . . . . . . . . . . . . . . 34 4.10. Comparativo curva de p´erdidas ´ohmicas . . . . . . . . . . . . . . . . . . 35 4.11. Comparativo curva de p´erdidas por concentraci´on de masa . . . . . . . . 35 4.12. Comparativo Curva de polarizacion . . . . . . . . . . . . . . . . . . . . . 36 4.13. Comparativo densidad potencia frente a densidad de corriente . . . . . . 37 4.14. Comparativo de curva potencial-eficiencia frente a la densidad de potencia 37 vii ´ INDICE DE FIGURAS viii En cualquier caso el flujo se considera incompresible. Adem´as en los canales puede existir agua liquida H2O(l), que provenga de las capas difusoras an´odica y cat´odica. Se supondr´a que la cantidad de agua que llega a los canales es muy peque˜na, de tal modo que se evacua por gravedad y no interfiere en el flujo de los reactantes en los canales. Evidentemente es necesario comprobar este supuesto (la obstrucci´on de los canales indicar´ıa un r´egimen de funcionamiento poco deseable) a posteriori. Las ecuaciones que se van a considerar (Navier-Stokes) son •Continuidad de la velocidad (de la mezcla de gases) (conservaci´on de masa) •Conservaci´on de cantidad de movimiento •Conservaci´on de los componentes gaseosos (ecuaciones de convecci´on-difusi´on para los gases, incluyendo difusi´on de Maxwell-Stefan) Capas difusoras: en esta zona de los dos electrodos constituida por material poroso y conductor el´ectrico (grafito), los gases se difunden en su camino hacia la zona catal´ıtica de los electrodos. En el flujo de medios porosos, la descarga espec´ıfica desempe˜na el papel de la velocidad. Proviene del promediado espacial de ´esta. Es la cantidad que proporciona el caudal. Se trabajan con las siguientes ecuaciones: •Conservaci´on de masa (Continuidad de la descarga espec´ıfica) •Conservaci´on cantidad de movimiento (Leyes de Darcy generalizadas) •Conservaci´on de los componentes gaseosos (ecuaciones de convecci´on-difusi´on para los gases, ley de difusi´on adaptada al medio poroso). Son necesarios, en las ecuaciones anteriores, distintos par´ametros que caracterizan los medios f´ısicos por los que se realizan los transportes y su interacci´on con los fluidos (permeabilidades, porosidad, coeficientes de difusi´on). Esto vale tambi´en para las restantes zonas de la pila. Membrana: en esta zona est´a el electrolito consistente en Nafi´on, con H2O(l), por el que circulan los iones H+ del ´anodo hacia el c´atodo. La naturaleza del electrolito impide el movimiento de los gases a trav´es de ´el. Se consideran las siguientes ecuaciones de transporte: •Continuidad de la descarga espec´ıfica (velocidad) 5 1. INTRODUCCI ´ ON •Conservaci´on cantidad de movimiento (Leyes de Sch¨ogl generalizadas, como las de Darcy, incluyendo el efecto del arrastre electrosm´otico del H2O(l) por los iones H+) •Conservaci´on iones H+ (ley de Darcy generalizada, incluyendo el t´ermino de arrastre debido a las diferencias de potencial el´ectrico en la membrana) •Conservaci´on de la energ´ıa (ecuaci´on de convecci´on-difusi´on para la temperatura, adaptada al medio poroso, m´as termino fuente ´ohmico debido a la densidad de corriente el´ectrica generada por los iones H+ ) •Ecuaci´on del potencial el´ectrico en la membrana •Ecuaci´on de la densidad de corriente en funci´on del potencial el´ectrico. De nuevo, son necesarios, en las ecuaciones anteriores, distintos par´ametros que caracterizan los medios f´ısicos por los que se realizan los transportes y su interacci´on con el H2O(l) y con los iones H+ (permeabilidades, porosidad, coeficientes de difusi´on, permeabilidades electrocin´eticas, conductividades el´ectricas). Capas catal´ıticas: en estas regiones tienen lugar las reacciones qu´ımicas de oxidaci´on y reducci´on que caracterizan la pila de combustible. Aqu´ı se encuentra el Pt finamente dividido que cataliza las reacciones de oxidaci´on y reducci´on. Una parte de Nafi´on que forma la membrana envuelve el catalizador para permitir la llegada de los reactantes. La zona es muy delgada en comparaci´on con los restantes dominios de la pila. Las ecuaciones, en principio an´alogas a las anteriores, incluir´an ahora t´erminos fuente (distintos en ´anodo y c´atodo) para las sustancias reactivas (H+, H2O, O2, H2). Estos t´erminos fuente se expresan como funciones de la densidad de corriente de transferencia del ´anodo (para las reacciones an´odicas) y del c´atodo (para las reacciones cat´odicas). Las densidades de corriente de transferencia se relacionan con el sobrepotencial mediante las ecuaciones de Butler-Volmer. Tambi´en se relacionan con el potencial el´ectrico mediante la conductividad i´onica. Como en la membrana y en las capas difusoras, es necesario conocer los distintos par´ametros que caracterizan los materiales por los que se realizan los transportes. Ecuaciones auxiliares: adem´as de la ecuaci´on de estado de los gases, son necesarias ecuaciones auxiliares adicionales para evaluar algunos de los mencionados par´ametros que caracterizan los materiales, cuyos valores pueden ser funci´on de variables termodin´amicas, como la temperatura, la presi´on y las concentraciones de alguna especie. Suelen considerarse ajustes funcionales a medidas experimentales. Son necesarias para evaluar, entre otras, la conductividad i´onica del Nafi´on, 6 la presi´on de vapor del agua o los coeficientes de difusi´on binarios (necesarios tambi´en para evaluar difusiones no fickianas). Por ´ultimo, tambi´en son necesarias expresiones para evaluar permeabilidades y coeficientes de difusi´on y conductividades t´ermicas en los medios porosos. Condiciones de contorno: es necesario en primer lugar prescribir los flujos y composici´on de los gases a la entrada de los canales. En las paredes de las placas bipolares que cierran los canales, se considera velocidad nula, flujo nulo de masa. Por ´ultimo, se consideran flujos nulos para las cantidades que no pueden atravesar entrefases en dichas entrefases. Consideraciones adicionales A la vista de lo expuesto, se comprende la dificultad de la descripci´on de los fen´omenos de transporte en una pila de combustible polim´erica. Existen en la literatura cient´ıfica diversos niveles de aproximaci´on a esta descripci´on, en ocasiones con hip´otesis diferentes para describir los mismos fen´omenos f´ısicos. Es fundamental ser equilibrado y consistente con el nivel de detalle escogido en la simulaci´on de todos los procesos en la pila. En las siguientes secciones se describen en detalle una propuesta de modelado con las correspondientes ecuaciones. 1.0.3. Entornos de simulaci´on Las ecuaciones que se han descrito a lo largo de este secci´on forman un sistema de ecuaciones diferenciales en derivadas parciales acopladas. Como se ha mencionado, dichas ecuaciones contienen t´erminos convectivos, difusivos y fuente, adem´as de la presi´on para el transporte de cantidad de movimiento. Este tipo de ecuaciones son habituales en mec´anica de fluidos y por tanto parece adecuada su simulaci´on num´erica mediante alg´un c´odigo de los utilizados en esa disciplina, adaptado a las caracter´ısticas propias de las ecuaciones que describen las pilas de combustible. En primer lugar es menester ser consistente con el nivel de detalle escogido en la simulaci´on de la pila y equilibrado en dicho nivel de detalle para los distintos procesos f´ısicos y qu´ımicos que tienen lugar en ella. Dado el elevado n´umero de ecuaciones existentes, incluso en las formulaciones m´as sencillas, es fundamental la comprobaci´on de que no hay ning´un “hueco” en la descripci´on escogida, que no se olvida ning´un par´ametro y que el sistema de ecuaciones es consistente y completo. Para que el m´etodo de discretizaci´on escogido sobre el mallado de la pila polim´erica sea de utilidad, es preciso que se satisfagan una serie de condiciones que se describen a 7 1. INTRODUCCI ´ ON continuaci´on. Consistencia: la discretizaci´on debe tender a ser exacta conforme el espaciado de la malla tiende a cero. Estabilidad: no se ampl´ıan los errores que aparecen en el curso de los procesos de soluci´on. Convergencia: la soluci´on de las ecuaciones discretizadas tiende a la soluci´on exacta conforme el espaciado de malla tiende a cero. Conservaci´on: en ausencia de fuentes y en estado estacionario, la cantidad de una magnitud conservada que entra en un volumen es igual a la que sale. Acotaci´on: las soluciones num´ericas deben tener valores comprendidos entre sus l´ımites f´ısicos. Realizabilidad: las soluciones deben tener sentido f´ısico Las soluciones son evidentemente aproximadas. Los errores provienen del modelado, discretizaci´on y de los proceso iterativos. Hay que ponderar la importancia de estos errores frente al tiempo de c´alculo extra necesario para minimizarlos. El problema de simular una pila polim´erica es, desde el punto de vista num´erico y para modelos realistas, extraordinariamente complicado, debido a la enorme cantidad de procesos que ocurren en el interior de la pila y que han de ser simulados. Hay que ser cuidadoso a la hora de escoger el c´odigo computacional sobre el que se va a resolver el sistema de ecuaciones resultante, atendiendo a los criterios expuestos en los p´arrafos anteriores. Idealmente se deber´ıa tener acceso a todas las l´ıneas del c´odigo. Esto sucede al utilizar c´odigos propios o cedidos con c´odigos fuente, probablemente la mejor soluci´on para el simulador con experiencia en m´etodos num´ericos. El uso de programas comerciales tiene el inconveniente de desconocer en muchos casos la naturaleza de los errores cometidos. Adem´as, en general se prima en exceso la estabilidad frente a la precisi´on. En este trabajo se ha escogido como c´odigo de trabajo OpenFOAM (http://www.openfoam.com/#openfoam), codigo CFD en v´olumenes finitos, de elementos poliedrales, no estructurado, que oferece una gran flexibilidad a la hora de dise˜nar un solver proprio para incluir los f´enomenos que ocurren dentro de una pila PEM. 8 1.0.4. Objetivos del trabajo El resultado de la soluci´on num´erica del conjunto de ecuaciones elegido para modelar la pila polim´erica es un conjunto de valores que representan la distribuci´on espacial (y temporal en su caso) para cada una de las magnitudes simuladas. Estos resultados han de ser validados compar´andolos con medidas experimentales existentes para condiciones de funcionamiento similares a las supuestas en el modelado. Tambi´en pueden servir para verificar la validez de esas suposiciones. Aunque hay muchos trabajos valiosos que enfocan aspectos particulares del comportamiento de la pila polim´erica, por ejemplo centr´andose en la capa difusiva o los electrodos, etc., el objetivo final es l´ogicamente la simulaci´on global del funcionamiento de la pila. Existe una intensa actividad cient´ıfica en este campo, que se beneficia de las continuas mejoras en los modelados, la caracterizaci´on param´etrica de los materiales y m´etodos num´ericos. La ventaja fundamental de las simulaciones num´ericas es que permiten cambiar par´ametros importantes como la permeabilidad o porosidad simplemente cambiando un n´umero en un archivo de datos, mientras que en un experimento esto mismo supondr´ıa el reemplazo de materiales por otros en algunos casos imposibles de adquirir. Una vez validado un modelo num´erico, se puede estudiar f´acilmente el efecto de la variaci´on de algunos par´ametros f´ısicos en la pila polim´erica. Tambi´en se mostrar´an los valores de algunas magnitudes importantes en la descripci´on del estado estacionario del funcionamiento de la pila, como las fracciones molares de vapor de agua, hidr´ogeno, ox´ıgeno, etc., as´ı como la distribuci´on espacial de la densidad de corriente y de los reactantes. 9 1. INTRODUCCI ´ ON 10 Cap´ıtulo 2 Modelo Matem´atico para un pila de tipo PEM Como se ha mencionado en la introducci´on, en todo modelo hay que llegar a un compromiso entre el nivel de detalle y el esfuerzo requerido para obtener resultados a partir de ´el. En el presente trabajo se pretenden demostrar las ventajas de una nueva aproximaci´on, para lo cual es conveniente la utilizaci´on de configuraciones geom´etricas e hip´otesis f´ısicas relativamente sencillas, como las que se muestran a continuaci´on: 1. Geometr´ıa 2D simplificada, pero se incluyen todos los componentes de la pila. 2. Caso isotermo. 3. Caso estacionario. 4. Densidad constante. 5. Difusi´on de Fick. 6. PEM de alta temperatura (no hay agua l´ıquida). 7. Hidr´ogeno y vapor de agua en la entrada an´odica y ox´ıgeno y vapor de agua en la cat´odica. 8. Se consideran p´erdidas por activaci´on, masa y ´ohmicas en la membrana. 9. Ecuaciones de transporte completas, incluyendo medio poroso. Este cap´ıtulo consta de dos partes, en la primera parte se exponen las ecuaciones matem´aticas que describen el movimiento de un fluido. En la segunda parte se exponen 11 2. MODELO MATEM ´ ATICO PARA UN PILA DE TIPO PEM las ecuaciones de campo el´ectrico y se explica como se acoplan estas ecuaciones en la fronteras para poder resolver el problema multif´ısico presente en este trabajo. La idea general es considerar las reacciones electroqu´ımicas que tienen lugar en las capas catal´ıticas como flujos salientes del dominio. La pila se dividir´ıa entonces en tres dominios, ´anodo, membrana y c´atodo. La informaci´on entre ellos vendr´a dada precisamente por el acoplo de los flujos en las fronteras. 2.1. Ecuaciones que describen el movimiento fluido. El ingeniero franc´es Claude Navier y el matem´atico ingl´es George Stokes escribieron las ecuaciones b´asicas que describen el movimiento de un fluido, a las cuales se les conoce como ¨ecuaciones de Navier-Stokes¨. Estas ecuaciones expresan en el lenguaje del medio continuo las tres leyes de conservaci´on b´asicas de la f´ısica: ecuaci´on de continuidad o conservaci´on de la masa, ecuaci´on de conservaci´on del movimiento y la ecuaci´on de conservaci´on de la energ´ıa. La ecuaci´on de continuidad se basa en la ley de conservaci´on de la masa. Aplicado al concepto de movimiento de un fluido, la tasa de variaci´on de la masa en un volumen de control es equivalente a la diferencia de la masa que entra y sale a trav´es de sus fronteras. La ecuaci´on de conservaci´on del movimiento se deriva de la aplicaci´on del concepto de la segunda ley de Newton a un fluido en movimiento. La ecuaci´on de movimiento se expresa en t´erminos de la presi´on y los esfuerzos debido a la viscosidad actuando sobre una part´ıcula fluida. La tasa de variaci´on de movimiento en una part´ıcula fluida es la diferencia de las fuerzas totales debido a los esfuerzos de la superficie y las fuerzas volum´etricas que act´uan sobre ella. Combinando estos principios fundamentales, el movimiento de un fluido se describe mediante un conjunto de ecuaciones en derivadas parciales conocidas como las ecuaciones de Navier-Stokes. Estas ecuaciones son una representaci´on matem´atica de las leyes de conservaci´on de la f´ısica. En el caso de un flujo laminar, estacionario y incompresible las ecuaciones tienen la siguiente forma: Ecuaci´on de continuidad. ∂ui ∂xi = 0 (2.1) 12 Ecuaci´on de Conservaci´on de movimiento. ui ∂ui ∂xj =−1 ρ ∂p ∂xi +ν∂2ui ∂xj∂xj (2.2) donde pes la presi´on, νla viscosidad cinem´atica, ρes la densidad. 2.1.1. Ecuaciones para el medio poroso Para un medio poroso (como en nuestro caso, por ejemplo, las capas difusoras) es conveniente distinguir entre dos tipos de promedios Ochoa-Tapia and Whitaker [1995]: Promedio superficial: h•i =1 VZV•dV,(2.3) Promedio intr´ınseco: h•iβ=1 VβZVβ •dVβ,(2.4) donde βindica la zona disponible para el que flujo se mueva libremente dentro del medio poroso (la parte ”vac´ıa”del medio poroso), y V indica un volumen lo suficientemente peque˜no que se va a usar para calcular el promedio. De acuerdo con la notaci´on Vβindica las zonas vac´ıas dentro del volumen V, su fracci´on esta dada por la porosidad ε, por definici´on. Por lo tanto los promedios superficiales e intr´ınsecos est´an relacionados por h•i =εh•iβ. Se puede demostrar que el promedio superficial de la velocidad es la cantidad que se acopla a la velocidad del flujo en el medio libre (canales) y que el promedio intr´ınseco de la presi´on dentro del medio poroso es la cantidad que se acopla a presi´on del flujo libre. Simplificando la notaci´on, uypvan a representar estas cantidades para el medio poroso. Para un aproximaci´on estacionaria, las ecuaciones de conservaci´on van a tener la forma: 1. Continuidad dentro medio poroso ∂uj ∂xj = 0 (2.5) 2. Cantidad de movimiento dentro medio poroso 1 ε2uj ∂ui ∂xj =−1 ρ ∂p ∂xi +ν ε ∂2ui ∂xj∂xj −ν Kui,(2.6) 13 2. MODELO MATEM ´ ATICO PARA UN PILA DE TIPO PEM donde pes la presi´on, νla viscosidad cinem´atica, εla porosidad y Kla permeabilidad, supuesto un medio poroso homog´eneo e is´otropo. El segundo t´ermino de la ecuaci´on 2.6 se conoce como la aproximaci´on de Brinkman, y el tercero refleja la contribuci´on de la ley de Darcy para el medio poroso. Se puede observar que las ecuaciones 2.1 y 2.2 son un caso particular de las ecuaciones 2.5 y2.6 cuando la porosidad tiene el valor ε= 1 y la permeabilidad K=∞. Como tal, las ecuaciones 2.5 y2.6 se pueden usar, Ochoa-Tapia and Whitaker [1995], para todo el dominio computacional sin la necesidad de imponer condiciones de contorno “internas” para las entrefases entre el medio poroso y el medio libre. 2.1.2. ´ Anodo Para el dominio an´odico, tenemos que resolver el campo de velocidades ui, la pre- si´on py la fracci´on m´asica de hidr´ogeno CH2. La fracci´on m´asica de vapor de agua ser´a CH2O= 1 −CH2. Las ecuaciones para el campo fluido utilizadas ser´an las de Navier-Stokes en el canal, promediadas para medio poroso en la capa difusora, vea la secci´on 2.1.1. En cuanto a la fracci´on m´asica del hidr´ogeno, CH2, su ecuaci´on de transporte es uj ∂CH2 ∂xj =∂ ∂xj γef ∂CH2 ∂xj ,(2.7) donde γef es el coeficiente de difusi´on efectivo, que depende del medio poroso. Seg´un Fishman and Bazylak [2011], γef =εε−εp 1−εpα γ(2.8) donde γes el coeficiente de difusi´on en medio libre, εpes el umbral de porosidad para que exista filtrado y αes un par´ametro de ajuste. Las condiciones de contorno son 1. Paredes u= 0 ∂CH2p ∂xn = 0 2. Entrada ρuea CH2e 14 habr´ıa que tomar el valor mayor de los encontrados, pues ser´ıa el potencial requerido para absorber el flujo de protones en el punto de mayor flujo. En el resto de puntos habr´ıa, por decirlo de manera coloquial, un “exceso” de concentraci´on de ox´ıgeno sobre el valor requerido para generar la corriente que no se aprovechar´ıa. 2.3.2. Comentarios En realidad la soluci´on anterior es un poco disipativa de m´as. La raz´on es que las p´erdidas por sobrepotencial an´odico son muy peque˜nas comparadas con las del sobrepotencial cat´odico, lo que para situaciones alejadas del “starving” de hidr´ogeno har´ıa m´as eficientes energ´eticamente y por tanto m´as realistas, soluciones en las que ηate´orico no fuera igual en todo el ´anodo, a cambio de que si fuera m´as uniforme en el c´atodo. Una posible aproximaci´on ser´ıa considerar flujo uniforme en el ´anodo (y c´atodo), de tal forma que, al ser abundante el ox´ıgeno en el c´atodo, el sobrepotencial cat´odico ser´ıa de valor muy uniforme, y por tanto habr´ıa menos disipaci´on energ´etica. Se deja como tema abierto, as´ı como la soluci´on en situaciones no extremas. 21 2. MODELO MATEM ´ ATICO PARA UN PILA DE TIPO PEM 22 Cap´ıtulo 3 Simulaci´on Num´erica 3.1. Geometr´ıa del Dominio En los cap´ıtulos anteriores se explic´o que para este trabajo se utilizar´a un modelo simpificado 2D. A continuaci´on se puede observar la geometr´ıa del dominio: Figura 3.1: Geometr´ıa del dominio. Y en el cuadro a continuaci´on se muestran las caracter´ısticas del domino: 23 3. SIMULACI ´ ON NUM´ ERICA Dimensi´on Valor Longitud total del domino 0,0712m Anchura de los canales Anodo y C´atodo 3,18 ·10−3m Anchura de la capa difusora (GDL) (2) y (6) 3,10 ·10−4m Anchura capa catal´ıtica (3) y (5) 1,0·10−5m Anchura de la membrana (4) 5,1·10−5m Cuadro 3.1: Dimensiones del dominio computacional 3.2. Malla computacional La malla que se ha utilizado para este trabajo es una malla de elementos cuadril´ateros y de un tama˜no aproximado de 8.600 elementos. Tambi´en se observa en la figura 3.3 que la raz´on de aspecto de la malla decrece conforme nos movemos hacia la capa catal´ıtica. De esta manera se garantiza que las variaciones de la magnitudes de inter´es son captadas correctamente, a la vez que se ahorran nodos computacionales. En la zona central el mallado es ya uniforme. Figura 3.2: Malla del domino computacional. 24 3.3. Simulaci´on num´erica Seg´un el modelo que se ha mostrado en el cap´ıtulo anterior, se distinguen b´asicamente dos tipos de regiones en la monocelda: regiones donde circulan fluidos (los canales m´as las capas difusoras -GDL-) y regiones donde circulan los protones (la membrana ´unicamente). Las capas catal´ıticas marcan la frontera entre ambos tipo de regiones, ya que en este trabajo se consideran infinitamente delgadas (l´ıneas en el caso 2D). En los canales y las capas difusoras (GDL) se aplican las ecuaciones de Navier-Stokes, con la formulaci´on de Ochoa en medios porosos (GDL) que evita la introducci´on de condiciones de contorno internas entre el canal y la capa difusora. En las capas catal´ıticas se generan las reacciones electrol´ıticas, que se expresan como flujos entrantes o salientes de los reactantes. En la membrana el flujo de protones en las fronteras se adec´ua a los anteriores mediante la correspondiente ley de conservaci´on y las ecuaciones de Butler- Volmer. El modelo desarrollado para la membrana bajo la hip´otesis de Onsager implica para los protones un movimiento unidireccional normal a la membrana, lo que se refleja como una simple ley de Ohm. Las subrutinas necesarias se han escrito como m´odulos para OpenFOAM, un c´odigo num´erico basado en vol´umenes finitos de libre distribuci´on, OpenFOAM. 3.3.1. Par´ametros de simulaci´on ´ Anodo C´atodo Combustible Hidr´ogeno Ox´ıgeno Permeabilidad K(m2) 2,584 x 10−13 2,584 x 10−13 Porosidad (ε) 0,517 0,517 Cref 0,909657 1,09013 * jref (A/m2) 10000 0,032 Constante Faraday (C/mol) 96485 96485 Temperatura (Kelvin) 353 353 Constante de gases R (J/molK) 8,314 8,314 Masa Molar (Kg/mol) 0,002 0,032 Presi´on kPa 101,1 101,1 N´umero de electrones involucrados α2 1 Densidad de mezcla ρ(Kg/m3) 0,08988 1,2 Cuadro 3.2: Cuadro par´ametros simulaci´on * Este valor por encima de 1 es consecuencia de tomar valores de densidad constantes 25 3. SIMULACI ´ ON NUM´ ERICA para las mezclas de gases, lo que obviamente no es la situaci´on real (ver cap´ıtulo 2) 26 Cap´ıtulo 4 Resultados Num´ericos En esta secci´on se discuten los resultados de la simulaci´on num´erica de la pila de combustible. Se obtienen la curva de polarizaci´on de la pila, curva de eficiencia, y otras curvas de importancia en el estudio del funcionamiento de una pila. En los resultados de algunas gr´aficas que se presentan a continuaci´on, se han separado los canales (´anodo y c´atodo) de forma artificial, ´unicamente mostrando el canal del ´anodo, GDLs, capas catal´ıticas, y canal del c´atodo, ya que en la membrana no se resuelve realmente ninguna ecuaci´on. Lo que se hace es un “mapping” de la cantidad de movimiento (corriente) saliente del ´anodo con el modelo asociado a conductividad baja, que implica movimiento unidireccional perpendicular a la membrana. A continuaci´on se presentan resultados para un caso espec´ıfico de corriente de 5,000A/m2 4.1. Velocidad En la Figura 4.1 se presentan los valores del m´odulo de la velocidad tanto para la regi´on an´odica como para la cat´odica. Para una mejor comprensi´on, se presenta esta cantidad primeramente utilizando un mapa de colores adecuado a las variaciones de su magnitud en el ´anodo (Figura 4.1 (a)), y posteriormente, adecuado a las variaciones de su magnitud en el c´atodo (Figura 4.1 (b)). 27 4. RESULTADOS NUM´ ERICOS (a) (b) Figura 4.1: Modulo de la Velocidad en el ´anodo(izquierdo) y c´atodo(derecho) 28 4. Resultados num´ericos En la figura 4.2 se observa el perfil de velocidad extra´ıdo para un corte transversal del ´andodo a media altura. Por las condiciones de contorno a la entrada y el consumo de H2correspondiente a la densidad de corriente demandada (5000A/m2), la magnitud de la velocidad que atraviesa la GDL es muy peque˜na comparada con la del canal en la parte central mostrada. De aqu´ı que el perfil de velocidad sea aproximadamente parab´olico en el canal, con una peque˜na perturbaci´on a la entrada de la capa difusora. Figura 4.2: Perfil de velocidad en corte transversal del ´andodo a media altura. En la figura 4.3 se observan los vectores de velocidad para las regiones del ´anodo y c´atodo. En la imagen aumentada se observa un detalle para mostrar con m´as claridad las direcciones de la velocidad. En la zona del ´anodo la GDL se ha coloreado de morado, la capa catal´ıtica de color amarillo, y se observa en estas zona la unidireccionalidad del flujo y que es normal a la superficie. Obs´ervese que esto est´a de acuerdo con el principio de m´ınima disipaci´on de Onsager Onsager [1931], debido al mayor car´acter disipativo de los medios porosos. En la regi´on del c´atodo se observa el mismo comportamiento que en el ´anodo. En este caso, la zona de color azul representa la GDL y la zona en rojo es la capa catal´ıtica. Se observa an´alogamente en estas zonas c´omo el flujo tambi´en es unidireccional y normal a la superficie en la frontera de la capa catal´ıtica. 29 4. RESULTADOS NUM´ ERICOS Figura 4.3: Vectores de velocidad ´anodo y c´atodo. En esta figura 4.4 se muestra el perfil de velocidad en salida capa catal´ıtica del ´anodo. Claramente se observa c´omo en la zona de entrada del canal la velocidad es mayor. Esto se debe b´asicamente a que la concentraci´on de H2es mayor en las zonas cercanas a la entrada y por tanto (Butler-Volmer) mayor la corriente generada que atraviesa la membrana. Por conservaci´on de masa tambi´en ha de se mayor el flujo de H2saliente y por tanto la velocidad. Figura 4.4: Perfil de velocidad salida capa catal´ıtica del ´anodo. 30 4. Resultados num´ericos (a) Barbir (b) Num´erica Figura 4.13: Comparativo densidad potencia frente a densidad de corriente Curva de potencial-eficiencia frente a la densidad de potencia En la figura 4.14 se presenta el potencial-eficiencia frente a la densidad de potencia. Se observa que al igual que la curva de potencia frente a la densidad de corriente, existe una potencia m´axima que la pila puede alcanzar, debido a que la eficiencia de la pila es directamente proporcional al potencial. Se pueden obtener mayores eficiencias con menores densidades de potencia, ya que el punto de operaci´on se puede seleccionar donde convenza o donde se necesite operar la pila. Tambi´en se puede observar que las gr´aficas de la figura 4.14 muestran un acuerdo cualitativo. (a) Barbir (b) Num´erica Figura 4.14: Comparativo de curva potencial-eficiencia frente a la densidad de potencia 37 4. RESULTADOS NUM´ ERICOS 38 Cap´ıtulo 5 Conclusiones En este trabajo se ha mostrado un nuevo tipo de modelizaci´on, basada en el principio de m´ınima disipaci´on de Onsager, que permite acoplar con fundamentos f´ısicos, los distintos fen´omenos de transporte que tienen lugar en una monocelda de una pila PEM. Se ha utilizado una geometr´ıa 2-D en una configuraci´on f´ısica relativamente sencilla que ya permite mostrar que esta aproximaci´on es correcta. Para ello se han implementado las necesarias subrutinas en un c´odigo num´erico sobre vol´umenes finitos de uso libre, OpenFOAM. Este acoplo distingue entre los canales m´as las capas difusoras (GDL), regiones donde circulan fluidos, y la membrana, regi´on donde circulan los protones. Las capas catal´ıticas marcan la frontera entre ambos tipos de regiones, ya que en este trabajo se consideran infinitamente delgadas (l´ıneas en el caso 2D). En los canales y las capas difusoras (GDL) se aplican las ecuaciones de Navier-Stokes, con la formulaci´on de Ochoa en medios porosos (GDL) que evita la introducci´on de condiciones de contorno entre el canal y la capa difusora. En las capas catal´ıticas se generan las reacciones electrol´ıticas, que se expresan como flujos entrantes o salientes de los reactantes. En la membrana el flujo de protones en las fronteras se adec´ua a los anteriores mediante la correspondiente ley de conservaci´on y las ecuaciones de Butler-Volmer. El modelo desarrollado para la membrana bajo la hip´otesis de Onsager implica para los protones un movimiento unidireccional normal a la membrana, lo que se refleja como una simple ley de Ohm. Los resultados num´ericos para una cierta intensidad de demanda t´ıpica se calculan en primer lugar. Las distintas magnitudes, velocidad y fracciones m´asicas de los reactantes (hidr´ogeno y ox´ıgeno), se muestran en las correspondientes secciones de la monocelda. Los resultados son congruentes con el comportamiento t´ıpico de una monocelda. Como validaci´on definitiva se calcula la curva de polarizaci´on y otras curvas caracter´ısticas de la monocelda que comparan adecuadamente con los resultados mos- 39 5. CONCLUSIONES trados por Barbir. Como trabajo futuro fundamental resta extender la formulaci´on para situaciones con membranas m´as conductoras, donde la aproximaci´on de corrientes prot´onicas perpendiculares a la membrana se relaje. Se hace notar, sin embargo, que a d´ıa de hoy, este continua siendo el caso para las membranas existentes para pilas PEM. En paralelo se pueden aplicar las aproximaciones mostradas en el presente trabajo a geometr´ıas 3D, con situaciones f´ısicas m´as complejas, que incluyan el efecto de las distintas densidades de las mezclas de gases o el transporte de calor en toda la monocelda o el transporte electr´onico en GDLs, incluyendo resistencias de contacto, etc. 40 Bibliograf´ıa F. Barbir, PEM Fuel Cells, Theory and practice, Elsevier Academic Press, 2005. 2,33, 34 J. A. Ochoa-Tapia, S. Whitaker, Momentum transfer at the boundary between a porous medium and a homeogenos fluid-i. theoretical development, International Journal of Heat and Mass Transfer 38 (6) (1995) 2635–2646. 13,14 Z. Fishman, A. Bazylak, Heterogeneous through-plane distributions of tortuosity, effective diffusivity, and permeability for pemfc gdls, Journal of The Electrochemical Society 158 (2) (2011) B247–B252. 14 K. Steinkamp, J. O. Schumacher, F. Goldsmith, M. Ohlberger, C. Ziegler, A nonisothermal pem fuel cell model including two water transport mechanismm in the membrane, Journal of Fuell Science and Technology 8 (1) (2008) 0110071–001100716. 17 L. Onsager, Reciprocal relations in irreversibles processes, Physical Review 37 and 97 (2 and 155) (1931) 406 and 1463. 19,29 W. Horne, K. Karamcheti, Extrema principles of entropy production and energy dissipation in fluid mechanics, in: Nasa technical memorandum, 1988, p. 100992. 19 J. A. Ochoa-Tapia, S. Whitaker, Momentum transfer at the boundary between a porous medium and a homeogenos fluid-ii. comparison with experiment, International Journal of Heat and Mass Transfer 38 (14) (1995) 2647–2655. H. Ju, C. Wang, Experimental validation of a pem fuel cell model by current distribution data, Journal of Electrochemical Society 151 (11) (2004) A1954–A1960. K. Lum, J. McGuirk, Three dimensional model of a complete polymer electrolyte membrane fuel cell-model formultaion, validation and parametric studies, Journal of Power Sources 143 (2005) 103–124. 41 BIBLIOGRAF´ IA F. J. Vald´es-Parada, B. Goyeau, J. A. Ochoa-Tapia, Jump momentum boundary condition at a fluid-porous dividing surface: Derivation of the closure problem, Chemical Engineering Science 67 (2007) 4025–4039. F. Barreras, A. Lozano, L. Vali˜no, C. Marin, A. Pascau, Flow distribution in a bipolar plate of a proton exchange membrane fuel cell: experiments and numerical simulation studies, Journal of Power sources 144 (2005) 54–55. F. Barreras, A. Lozano, L. Vali˜no, R. Mustata, C. Mar´ın, Fluid dynamics performance of different bipolar plates part i.velocity and pressure fields, Journal of Power sources 175 (2) (2008) 841–850. A. Lozano, L. Vali˜no, F. Barreras, R. Mustata, Fluid dynamics performance of different bipolar plates part ii. flow through the diffusion layer, Journal of Power sources 179 (2008) 711–722. L. Carrette, K. Friedrich, Fuels cells-fundamentals and applications, fuel cells 1 (1) (2001) 5–39. E. Hontan, M. Escudero, C. Bautista, L. Garca-Ybarra, Optimization of flow-field in polymer electrolyte membrane fuel cells using computational fluid dynamics techniques, Journal of Power Source 86 (1-2) (2000) 363–368. R. Perry, D. Green, J. Maloney, Chemical Engineers HandbookPEM Fuel Cells, Theory and practice, sixth edition Edition, Mc Graw-Hill Ed., 1984. F. White, Fluid Mechanics, Mc Graw Hill Ed., USA, 1979. J. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, 2nd Edition, Ed. Springer-Verlag, Berlin, 1999. 42