Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Study of fluid flow in the porous media of gas diffusion layers in proton exchange membrane fuel cells Shahraeeni, Mehdi 2013

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2013_spring_shahraeeni_mehdi.pdf [ 6.37MB ]
Metadata
JSON: 24-1.0073812.json
JSON-LD: 24-1.0073812-ld.json
RDF/XML (Pretty): 24-1.0073812-rdf.xml
RDF/JSON: 24-1.0073812-rdf.json
Turtle: 24-1.0073812-turtle.txt
N-Triples: 24-1.0073812-rdf-ntriples.txt
Original Record: 24-1.0073812-source.json
Full Text
24-1.0073812-fulltext.txt
Citation
24-1.0073812.ris

Full Text

Study of Fluid Flow in the Porous Media of Gas Diffusion Layers in Proton Exchange Membrane Fuel Cells by Mehdi Shahraeeni  B.A.Sc., University of Tehran, 2005 M.A.Sc., University of Tehran, 2007  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE COLLEGE OF GRADUATE STUDIES (Mechanical Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Okanagan) April 2013 c Mehdi Shahraeeni, 2013  Abstract A proton exchange membrane (PEM) fuel cell is an energy converting system generating electricity by oxidation of hydrogen and reduction of oxygen with water and heat as the only waste products. Despite the huge market potential of the fuel cell, its performance and cost must be improved significantly before constituting a viable market. One of the major problems of current fuel cells is water management: at energy demanding conditions where the cell is operating at high current densities, excessive water produced restricts the access of reactant gases and hence reduces the performance of the cell. To improve water management, it is necessary to study water transport mechanisms in the internal network of the cell, especially in the porous gas diffusion layer (GDL) through which transport of electrons, reactant gases, and water occurs. In this thesis, fluid flow through the GDL is studied experimentally and numerically using fluorescence microscopy and a pore network modeling approach, respectively. The images obtained from the microscope are analyzed to find patterns of flow inside the GDL samples with different hydrophobicity. Three different flow patterns are observed: initial invasion, progression, and pore-filling. The observations show that liquid water flows into the majority of available pores on the boundary of the hydrophilic (untreated) GDL and several branches segregate from the initial pathways. For the hydrophobic (treated) GDL, however, a handful of boundary pores are invaded and the original pathways extend toward the other side of the medium with minimum branching. In addition to flow visualization, the experimental setup facilitates the precise measurement of pressure and time of breakthrough which are used as boundary condition and the validation criterion for the numerical simulation, respectively. The numerical model, developed based on an invasion percolation algorithm, is used to study the effects of GDL hydrophobicity and thickness on the flow configuration and breakthrough time as well as to determine the flow rate and saturation in different GDL samples. The developed model can be used to optimize the GDL properties for designing porous medium with an effective transport characteristic.  ii  Preface I am the principal author of this dissertation. I developed the idea of the research, designed and built the experimental setup configuration, and defined and regulated the experimental procedure for collecting data. Graduate student, Siddiq Husain Tahseen, and undergrad student, Samuel Yew, contributed in measuring the expansion of setup components. I interpreted the experimental observations based on which the numerical scheme of this work is developed. I designed the pore-network algorithm, checked its correctness, implemented and programmed it, and generated the numerical results. My academic advisor, Dr. Mina Hoorfar, co-authored all of the chapters of this thesis, as well as the publications resulted from this work which are listed:  Journal Articles − M. Shahraeeni and M. Hoorfar, Experimental and numerical comparison of water transport in untreated and treated diffusion layers of proton exchange membrane (PEM) fuel cells, Journal of Power Sources, 238, 29-47, 2013. − M. Shahraeeni and M. Hoorfar, Effect of Gas Diffusion Layer Properties on the Time of Breakthrough, Journal of Power Sources, 196(14), 5918-5921, 2011. − M. Shahraeeni and M. Hoorfar, Pore-network modeling of liquid water flow in the Gas Diffusion Layer of proton exchange membrane (PEM) fuel cells, International Journal of Hydrogen Energy, In preparation.  Conference Proceedings − M. Shahraeeni and M. Hoorfar, Pore-network Modeling of Fluid Flow in Thin Porous Media: Application in PEMFC, Proceedings of ICNMM2012, 10th International Conference on Nanochannels, Microchannels, and Minichannels, July 2012, Puerto Rico, USA. iii  Preface − M. Shahraeeni and M. Hoorfar, Capillary Pressure-Saturation at the Time of Breakthrough for Thin Porous Media of GDLs, Proceedings of ICNMM2012, 10th International Conference on Nanochannels, Microchannels, and Minichannels, July 2012, Puerto Rico, USA. − S.H.Tahseen, K. Chen, M. Shahraeeni, S.C.M. Yew and M. Hoorfar, Measurement of Liquid Water Content inside the Gas Diffusion Layer, Proceedings of ESFuelCell2012, 10th Fuel Cell Science, Engineering and Technology Conference, July 2012, San Diego, USA. − M. Shahraeeni and M. Hoorfar, Simulation of Fluid Flow in the Gas Diffusion Layer of PEM fuel cell, Proceedings of WHEC2012, 19th World Hydrogen Energy Conference, June 2012, Toronto, Canada. − S.H.Tahseen, K. Chen, M. Shahraeeni, S.C.M. Yew and M. Hoorfar, Measurement of Liquid Water Content inside the Gas Diffusion Layer after the Time of Breakthrough, Proceedings of WHEC2012, 19th World Hydrogen Energy Conference, June 2012, Toronto, Canada. − M. Shahraeeni and M. Hoorfar, Numerical Investigation of Fluid Flow Inside the Porous Medium of GDL, Proceedings of ICNMM2010, 8th International Conference on Nanochannels, Microchannels, and Minichannels, August 2010, Montreal, Canada. − M. Shahraeeni and M. Hoorfar, Effect of PTFE Loading on Performance of the GDL in Water Removal, Proceedings of FUELCELL2010, 8th International Fuel Cell Science, Engineering and Technology Conference, June 2010, New York, USA. − M. Shahraeeni and M. Hoorfar, A Pore-network Model for Gas Diffusion Layers of PEM Fuel Cells, Proceedings of EFC2009, 3rd European Fuel Cell Technology & Applications Conference, December 2009, Rome, Italy. − B. R. Friess, M. Shahraeeni and M. Hoorfar, Water Management in PEM Fuel Cells, AIChE Conference 2008, Philadelphia, USA.  iv  Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xvi Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii Chapter 1: Introduction . . . . . . . . . . . . . . . . . . 1.1 Basic operation of PEM fuel cells . . . . . . . . . . . 1.1.1 The Gas Diffusion Layer . . . . . . . . . . . . 1.1.2 Flooding . . . . . . . . . . . . . . . . . . . . . 1.2 Literature review and theory . . . . . . . . . . . . . 1.2.1 Modeling of fluid flow in porous media . . . . 1.2.2 Network models of fluid flow in Gas Diffusion 1.2.3 Experimental techniques . . . . . . . . . . . . 1.3 Hypothesis . . . . . . . . . . . . . . . . . . . . . . . 1.4 Thesis organization . . . . . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . . . . . . . . . . . Layers . . . . . . . . . . . .  . . . . . . . . . .  1 3 5 6 8 9 11 16 25 27  Chapter 2: Experimental measurements 2.1 Setup . . . . . . . . . . . . . . . . . 2.1.1 Fluorescence microscopy . . . 2.1.2 Injection channel . . . . . . . 2.1.3 Procedure . . . . . . . . . . . 2.2 Results . . . . . . . . . . . . . . . . .  . . . . . .  . . . . . .  28 29 30 31 32 34  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . . .  . . . . . .  . . . . . .  v  TABLE OF CONTENTS 2.3  Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Flow characteristics . . . . . . . . . . . . . . . . . . . 2.3.2 Pressure and time at breakthrough . . . . . . . . . . .  48 48 52  Chapter 3: Numerical model . . . . . . . . . . . . . . . . . . . . 3.1 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Network generation . . . . . . . . . . . . . . . . . . . 3.1.2 Network treatment . . . . . . . . . . . . . . . . . . . . 3.1.3 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Comparison between the flow patterns obtained from numerical results and experimental observations . . . 3.3.2 Effect of hydrophobicity on the flow pattern . . . . . .  59 61 62 64 66 74 101  Chapter 4: Model implementation . . . . . . . . . . . 4.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Transient analysis of fluid flow . . . . . . . . 4.2.2 Analysis of breakthrough . . . . . . . . . . .  112 114 125 125 127  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  101 104  Chapter 5: Conclusions and future work . . . . . . . . . . . . . 135 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 Appendix A: Code structure . . . . . . . . . . . . . . . . . . . . . . 151  vi  List of Tables Table 3.1 The geometrical and hydrodynamic properties of the networks and fluids . . . . . . . . . . . . . . . . . . . . 67 Table 3.2 Normalized time of breakthrough obtained from the simulations (Averages of the replicates are provided in brackets). . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Table 4.1 Normalized time of breakthrough and breakthrough pressure (experiment). . . . . . . . . . . . . . . . . . . . . . 114 Table 4.2 Network properties and boundary condition. . . . . . . 114 Table 4.3 Normalized values of time of breakthrough obtained from numerical model; For each case the average value is provided (second row) with the value of experimental breakthrough pressure used as boundary condition (in brackets). . . . . . . . . . . . . . . . . . . . . . . . . . . 127  vii  List of Figures Figure 1.1 Figure 1.2 Figure 1.3  Figure Figure Figure Figure Figure  Figure  Figure  Figure  Figure  Figure  Figure  2.1 2.2 2.3 2.4 2.5  Different layers in a PEM fuel cell . . . . . . . . . . . Schematic of a PEM fuel cell . . . . . . . . . . . . . . Typical polarization curve of a PEM fuel cell (Schematic of polarization curve: Courtesy of UC-San Diego). . .  Experimental setup. . . . . . . . . . . . . . . . . . . . Details of fluorescence microscopy. . . . . . . . . . . . Details of injection channel. . . . . . . . . . . . . . . . The placement of the sample. . . . . . . . . . . . . . . Change in the volume of the setup versus capillary pressure: low flow rate (solid line), high flow rate (dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Flow configuration and capillary pressure variation versus time (sample thickness=110µm, loading=0%, 3 flow rate=1.11 × 10−11 ms ) . . . . . . . . . . . . . . . . 2.7 Flow configuration and capillary pressure variation versus time (sample thickness=110µm, loading=40%, 3 flow rate=1.11 × 10−11 ms ) . . . . . . . . . . . . . . . . 2.8 Flow configuration and capillary pressure variation versus time (sample thickness=110µm, loading=0%, 3 flow rate=1.11 × 10−10 ms ) . . . . . . . . . . . . . . . . 2.9 Flow configuration and capillary pressure variation versus time (sample thickness=110µm, loading=40%, 3 flow rate=1.11 × 10−10 ms ) . . . . . . . . . . . . . . . . 2.10 Flow configuration and capillary pressure variation versus time (sample thickness=190µm, loading=0%, 3 flow rate=1.11 × 10−11 ms ) . . . . . . . . . . . . . . . . 2.11 Flow configuration and capillary pressure variation versus time (sample thickness=190µm, loading=40%, 3 flow rate=1.11 × 10−11 ms ) . . . . . . . . . . . . . . . .  2 3 5 36 37 38 38  39  40  41  42  43  44  45  viii  LIST OF FIGURES Figure 2.12 Flow configuration and capillary pressure variation versus time (sample thickness=190µm, loading=0%, 3 flow rate=1.11 × 10−10 ms ) . . . . . . . . . . . . . . . . Figure 2.13 Flow configuration and capillary pressure variation versus time (sample thickness=190µm, loading=40%, 3 flow rate=1.11 × 10−10 ms ) . . . . . . . . . . . . . . . . Figure 2.14 Comparison of initial invasion patterns: Less surface area of the hydrophobic sample is invaded (left) compared to the hydrophilic area (right) (the encircled areas compare the invaded clusters). . . . . . . . . . . . Figure 2.15 Progression pattern for hydrophobic and hydrophilic samples (the encircled areas in Figures (b) and (d) show the new developed water frontier; note that for the hydrophobic sample (b) only three more clusters are invaded as water progresses, while for the hydrophilic sample (d) the number of invaded clusters is twice. . . Figure 2.16 Pore filling process for hydrophobic and hydrophilic GDLs; the encircled areas show the high saturated regions. For the hydrophobic sample four clusters are indicated as high saturated pores, while for the hydrophilic sample only one relatively-high saturated cluster is visible. . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.17 Measured pressure and time of breakthrough versus 3 the thickness of the sample (flow rate of 1.11×10−11 ms ) for different PTFE loadings. . . . . . . . . . . . . . . . Figure 2.18 Measured pressure and time of breakthrough versus 3 the thickness of the sample (flow rate of 1.11 × 10−10 ms ). Figure 2.19 Variations of the pressure and time over PTFE loadings for the thickest sample (370µm) (flow rate of 1.11× 3 10−11 ms ). . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.20 Variations of the pressure and time over PTFE loadings for the thickest sample (370µm) (flow rate of 1.11× 3 10−10 ms ). . . . . . . . . . . . . . . . . . . . . . . . . . Figure 3.1  Consecutive steps of pore-network modeling approach (Representation of the porous matrix: Courtesy of Institute of Energy Technology, Department of Mechanical and Process Engineering, Swiss Federal Institute of Technology Zurich ). . . . . . . . . . . . . . . . . . .  46  47  50  54  55  55 56  57  58  62  ix  LIST OF FIGURES Figure 3.2 Figure 3.3 Figure 3.4  Figure 3.5  Figure 3.6 Figure 3.7 Figure 3.8 Figure 3.9 Figure 3.10 Figure 3.11 Figure 3.12 Figure 3.13 Figure 3.14 Figure 3.15 Figure 3.16 Figure 3.17 Figure 3.18 Figure 3.19 Figure 3.20 Figure 3.21  Configuration of the unit cell and the adjacent cells. . A sample of pore-network model with index notations (15 × 15 × 10). . . . . . . . . . . . . . . . . . . . . . . Initial flow rate versus total exposed area for the numerical simulations; the distribution suggests that the initial flow rate is reduced by a factor of 100, as the network exposed area is reduced with the same factor. Treatment of networks with the fraction of hydrophobic agent: (a)f = 0.10, (b)f = 0.20, (c)f = 0.30 and (d)f = 0.40 . . . . . . . . . . . . . . . . . . . . . . . . Flowchart of the displacing algorithm. . . . . . . . . . Connecting throat between two adjacent cell. . . . . . Three replications of treated network (f = 0.4) . . . . Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.01. . . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.25. . . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.50. . . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.75. . . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 1.00. . . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 0.01. . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 0.25. . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 0.50. . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 0.75. . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 1.00. . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 0.01. . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 0.25. . . . . . . . . . . . . . . . Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 0.50. . . . . . . . . . . . . . . .  63 64  65  68 72 73 75 76 77 78 79 80 81 82 83 84 85 86 87 88  x  LIST OF FIGURES Figure 3.22 Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 0.75. . . . . . . . . . . . . . . . Figure 3.23 Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 1.00. . . . . . . . . . . . . . . . Figure 3.24 Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 0.01. . . . . . . . . . . . . . . . Figure 3.25 Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 0.25. . . . . . . . . . . . . . . . Figure 3.26 Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 0.50. . . . . . . . . . . . . . . . Figure 3.27 Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 0.75. . . . . . . . . . . . . . . . Figure 3.28 Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 1.00. . . . . . . . . . . . . . . . Figure 3.29 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.01. . . . . . . . . . . . . . . . Figure 3.30 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.25. . . . . . . . . . . . . . . . Figure 3.31 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.50. . . . . . . . . . . . . . . . Figure 3.32 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.75. . . . . . . . . . . . . . . . Figure 3.33 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 1.00. . . . . . . . . . . . . . . . Figure 3.34 Numerical results demonstrating the initial invasion pattern at z ∗ = 0.25: (a) f = 0.0 (b) f = 0.10 (c) f = 0.20 (d) f = 0.30 (e) f = 0.40. . . . . . . . . . . . Figure 3.35 Progression pattern for the treated sample at: (a) t∗ = 0.25 ; (b) t∗ = 0.50. . . . . . . . . . . . . . . . . . Figure 3.36 Progression pattern for the untreated sample at: (a) t∗ = 0.25 ; (b) t∗ = 0.50. . . . . . . . . . . . . . . . . . Figure 3.37 Numerical results demonstrating the pore-filling pattern at z ∗ = 0.25 for: (a) treated sample; and (b) untreated sample. . . . . . . . . . . . . . . . . . . . . . Figure 3.38 Effect of hydrophobic fraction (f ) on flow configuration at breakthrough. . . . . . . . . . . . . . . . . . . . Figure 3.39 Effect of hydrophobic fraction (f ) on saturation at the breakthrough. . . . . . . . . . . . . . . . . . . . . . . . Figure 3.40 Effect of hydrophobic fraction on the time of breakthrough. . . . . . . . . . . . . . . . . . . . . . . . . . .  89 90 91 92 93 94 95 96 97 98 99 100  102 106 107  108 109 110 111 xi  LIST OF FIGURES Figure 4.1  Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.01. . . . . . . . . . . . . . . . . Figure 4.2 Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.25. . . . . . . . . . . . . . . . . Figure 4.3 Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.50. . . . . . . . . . . . . . . . . Figure 4.4 Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.75. . . . . . . . . . . . . . . . . Figure 4.5 Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 1.00. . . . . . . . . . . . . . . . . Figure 4.6 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.01. . . . . . . . . . . . . . . . Figure 4.7 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.25. . . . . . . . . . . . . . . . Figure 4.8 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.50. . . . . . . . . . . . . . . . Figure 4.9 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.75. . . . . . . . . . . . . . . . Figure 4.10 Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 1.00. . . . . . . . . . . . . . . . Figure 4.11 Effect of hydrophobic fraction on flow configuration at breakthrough. . . . . . . . . . . . . . . . . . . . . . Figure 4.12 Numerical comparison of through-plane saturation profiles at different time steps for hydrophilic (solid line) and hydrophobic (dash line) samples. . . . . . . . . . . Figure 4.13 Analysis of error variance for the numerical results. . . Figure 4.14 Comparison of the numerical and experimental results for time of breakthrough (the shadow regions show the error interval for each case). . . . . . . . . . . . . . . Figure 4.15 Comparison of flow rates for hydrophobic and hydrophilic samples. . . . . . . . . . . . . . . . . . . . . . Figure 4.16 Effect of sample thickness on total saturation. . . . . .  115 116 117 118 119 120 121 122 123 124 126  130 131  132 133 134  xii  List of Symbols cr Ca dp , dt ∆t e F f I i i0 i, j, k IM AX JM AX KM AX k K L lthroat ls m M M mf iber mP T F E msample muntreated N n P , Pc P∗  Surface concentration of reactant Capillary number Pore diameter and throat diameter Simulation time step Charge of an electron Faraday constant Hydrophobic fraction Electric current Current density Exchange current density Pore indexes Number of pores in x−direction Number of pores in y−direction Number of pores in z−direction Boltzmann constant Permeability Tube length(Poiseuille equation) Throat length Length scale Mass of reactant (Faraday’s law) Molar mass of reactant Mobility ratio Fiber mass in the GDL sample Mass of PTFE in the GDL sample Mass of the GDL sample Mass of untreated GDL sample Total number of hydrophobic pores Number of displaced ions (Nernst equation) Fluid pressure and capillary pressure Normalized pressure  (mol/m2 ) (−) (µm) (s) (C) (C/mol) (−) (A) (A/cm2 ) (A/cm2 ) (−) (−) (−) (−) (m2 kg/s2 K) (m2 ) (m) (µm) (m) (kg) (kg/kmol) (−) (kg) (kg) (kg) (kg) (−) (−) (P a) (−)  xiii  List of Symbols Pcell Phydrodynamic Pref Q , Qnom Qthroat ˜ Q R Ri R1 , R2 rp , rt Re s Spore , Scell T t∗ Thydrophilic Thydrophobic tBT tcalibration V V Vcell Vexp Vf iber Vf luid Vpore Vtotal Vwater x ∗ ∗ x , y , z∗ z  Netowrk cell pressure Hydrodynamic pressure Reference (reservoir) pressure Flow rate, nominal flow rate Throat flow rate in the network Time-averaged flow rate Ideal gas constant Resistance of cell component Principal radii of curvature Pore radius and throat radius Reynold number Phase saturation Pore saturation Absolute temperature Normalized time Breakthrough time (hydrophilic sample) Breakthrough time (hydrophobic sample) Time of breakthrough Calibration time Velocity (definition of Ca number) Volume Volume of a generic cell in the network Expansion volume of setup components Volume of the fiber in the GDL sample Fluid volume within a particular pore (cell) Pore volume in the GDL sample Total volume of the GDL sample Total volume of water in the GDL sample Mole fraction of reactants (Nernst equation) Dimensionless cross sections in the network Number of transferred electrons (Faraday’s law)  (P a) (P a) (P a) (m3 /s) (−) (m3 /s) (J/molK) (Ω) (m) (µm) (−) (−) (−) (K) (−) (s) (s) (s) (s) (m/s) (m3 ) (m3 ) (m3 ) (m3 ) (m3 ) (m3 ) (m3 ) (m3 ) (−) (−) (−)  xiv  List of Symbols Greek symbols α δ ηact ηmass ηohm γ µ ω, %wt ρ θ  Charge transfer coefficient Thickness of PTFE layer covering a single pore Porosity Activation loss Mass transport loss Ohmic loss Surface/interfacial tension Viscosity Weight percentage of PTFE loading in the sample Density Contact angle  (−) (µm) (−) (V ) (V ) (V ) (mN/m) (N s/m2 ) (−) (kg/m3 ) (−)  xv  Acknowledgements I would like to thank my academic advisor, Dr. Mina Hoorfar, for her help and support during my work at UBC Okanagan. I would also like to thank the advisory committee Drs. Homayoun Najjaran and Claire Yu Yan, for their suggestions and advice throughout the course of my study. I would like to thank my friends, office mates and lab mates in Kelowna: Ali, Nima, Sina, Amin, Mojtaba, Milad, AliEb and Maryam. I would like to thank my brothers, Drs. Ebrahim Shahraeeni and Mohammad Sadegh Shahraeeni, and my sister-in-laws, Dr. Leila Sharifian and Leila Mohammadi for providing warm support and encouragement specially in the final stages of my PhD work. The financial supports from the British Columbia Innovation Council through the BCIC Graduate Scholarship, and BC Ministry of Advanced Education, Innovation and Technology through the Pacific Century Graduate Scholarship are gratefully acknowledged.  xvi  Dedication To my mother, Kobra...  xvii  Chapter 1  Introduction Global warming and the shortage of hydrocarbon fuels as well as the instabilities in their production are among main reasons that make researchers seek for new and reliable energy alternatives. Burning hydrocarbons directly to release chemical energy is still the main way for producing power in industry and transportation. If the heat produced during the burning process is not appropriately controlled and re-used, it would be a waste of exergy. To recover this waste energy, co-generation systems and combined-power cycles were developed. Based on the Carnot cycle principle, the efficiency of such systems depends on the temperature of the hot and cold reservoirs. The larger the difference between the temperatures of the hot and cold reservoirs, the higher the efficiency of the engine. Since there are limitations to decreasing the temperature of the cold reservoir, the temperature of the hot reservoir needs to be increased to have higher efficiency. This is equivalent to direct burning of hydrocarbons which results in production of nitrogen oxides (N Ox ) and carbon dioxide (CO2 ), causing air pollution and contributing to atmospheric green house gases (GHG). As a short–term solution, releasing the energy of hydrocarbons at a lower temperature is favorable though the efficiency of the system is reduced. The long–term solution is to find an environmentally friendly energy source which produces energy in low temperatures with zero emissions. Hydrogen is the best choice as the by-product of the oxidation of hydrogen is water. Proton Exchange Membrane (PEM) fuel cells are primarily designed to use hydrogen as fuel. The working conditions of a PEM fuel cell, as well as the wide range of power that it could provide, make it ideal as a future alternative energy converter. The essentials of a PEM fuel are very simple, with almost no moving parts. This can lead to high reliability and a long-lasting system as well as a low level of noise in operation. The latter feature is crucial in portable power applications and particularly for local power generation in combined heat and power schemes. A PEM fuel cell consists of a membrane electrode assembly (MEA) sandwiched between two flow channels, as presented schematically in Figure 1.1. The MEA contains a polymer electrolyte membrane (e.g. Nafion R ) embed1  Chapter 1. Introduction  Figure 1.1: Different layers in a PEM fuel cell ded between two porous gas diffusion electrodes (GDE). The GDE is composed of a platinum catalyst and a gas diffusion layer (GDL) constructed from macro-porous substrates (i.e., carbon fiber or carbon cloth impregnated with polytetrafluroethylene (PTFE)) coated with one or more micro-porous layers (i.e., amorphous mixture of carbon and PTFE) [FSH08]. The manufacturing technology of PEM fuel cells has come a long way in the last few years. However, their performance still must be significantly improved before they can contribute to a sustainable market. Recent experimental and numerical investigations identify water management as a critical factor in the design of robust and high-efficiency fuel cells. In essence, the ionic conductivity of the electrolyte is dependent on the hydration level of the membrane as water molecules transport hydrogen ions across the electrolyte. However, excessive water vapor condensation, due to lengthy operation or large output current, forms micro-droplets that cover the active sites on the catalyst layers, fill the pores of the GDL, and block access of reactant gases to the reaction site. Typically, this is the origin of a limiting current for PEM fuel cells [FSH08]. Thus, water management is one of the major problems for the fuel cell on the way to being commercialized. To improve water management, it is necessary to understand water transport mechanisms in the internal network of the cell, especially the porous medium of the GDL through which transport of electrons, reactant gases, and water occurs. The next section (1.1) introduces basic operation of a PEM fuel cell, its components and the role of each component in the operation of the cell, and the polarization curve as the main performance characterization method. A description of the material used for gas diffusion layers (GDLs), as the target component of the cell in this research, is presented. Then the electrochemistry of the flooding phenomenon, the sources of the voltage loss in 2  1.1. Basic operation of PEM fuel cells a cell, and the effect of flooding on the performance of the cell are thoroughly discussed. The section is summarized by linking the flooding to the fluid flow inside the GDL porous medium. Then, a thorough review of the current models and experimental techniques used to study fluid flow in fuel cell components is presented in Literature review and theory (1.2). Experimental visualization techniques focused on the formation and flow of liquid water in the PEM fuel cell and its components are thoroughly reviewed and discussed in Section 1.2.3. A brief introduction on the mechanism of water droplet formation and its role in limiting the current of the cell is restated in the Hypothesis section (1.3). This section also includes the motivation of this research as well as the experimental and numerical approaches adopted. The chapter is finally closed with a section describing the organization of the thesis.  1.1  Basic operation of PEM fuel cells  The basic operation of a PEM fuel cell is presented in Figure 1.2. Oxygen  Figure 1.2: Schematic of a PEM fuel cell and hydrogen are fed through the flow channels to the cathode and anode sides, respectively. The hydrogen fuel is delivered through the gas diffusion layer (GDL) to the catalyst layer. At the catalyst layer, the H2 molecules are converted to protons and electrons by means of the catalyst (usually Pt particles deposited on a porous substrate or membrane). The protons are 3  1.1. Basic operation of PEM fuel cells transferred to the cathode side through the membrane which is only proton conductive. The electrons flow through the external circuit and produce power. These electrons are conducted through the GDL to the surface of the catalyst at the cathode side. On this side oxygen, electrons and hydrogen ions react and produce water molecules. The polymer electrolyte membrane is a Nafion R -based material. This layer is primarily based on polyethylene, but the hydrogen molecules are substituted with fluorine to form polytetrafluroethylene (PTFE) or Teflon. The basic PTFE polymer is then “sulphonated” by adding a side chain that ends with sulphonic acid HSO3 . The HSO3 group is ionically added to the main chain results in the existence of an SO3− ion at the end of this side chain. This negative ion of SO3− attracts the H + ion and makes a strong mutual bond between the opposite charges of each molecule. Thus, the final polymer conducts protons from the anode side to the cathode side, depending on the hydration of the membrane [LDF03]. The catalyst layer is a combination of carbon powder and platinum particles. This mixture is deposited on both sides of the membrane to form the membrane electrolyte assembly (MEA). In an alternative method, the catalyst powder is deposited on the surface of the GDL to form a gas diffusion electrode (GDE). The gas diffusion layers (GDLs) are carbon fiber products (carbon paper or carbon cloth). In addition to providing mechanical support for the membrane, the GDLs conduct the reactant gas (oxygen at the cathode side and hydrogen at the anode side) to the reaction sites, where the oxidation and reduction take place. The GDLs must also be electrically conductive to transfer the electrons from the catalyst layer to the gas channels in the anode side and vice versa in the cathode side. Another characteristic of the GDLs is porosity, as the layers provide paths to deliver gases to the reaction sites on the catalyst layers as well as deliver the water produced at the cathode to the corresponding gas channel. Water is collected and evaporated from the interface of the GDL and the gas channel by means of the air which flows parallel to the interface. However, at high current densities, when the rate of reaction increases, water droplets are not removed from this interface at an appropriate rate. Thus, they partially block the surface of the catalyst layer, and if the cell keeps on working, a layer of water covers the whole active area of the catalyst layer. This results in a significant voltage drop. A typical polarization curve presented in Figure 1.3 represents this sudden voltage drop at high current densities. This problem is referred to as “flooding” which is one of the major issues limiting the performance and impeding the commercialization of industry-standard fuel cells. In this research, the role of the GDL in the flooding phenomenon is studied. The 4  1.1. Basic operation of PEM fuel cells  Figure 1.3: Typical polarization curve of a PEM fuel cell (Schematic of polarization curve: Courtesy of UC-San Diego). research is mainly focused on the flow of liquid water inside the GDL and the effect of different GDL properties on the flow.  1.1.1  The Gas Diffusion Layer  The Gas Diffusion Layer (GDL) is a porous diffusive medium made from carbon fiber. The layer is sandwiched between the gas channel and the catalyst layer. The main roles of the GDL are: 1) providing mechanical support for the membrane, 2) providing passages for the reactant gases to be delivered to the reaction sites and 3) conducting the electrons between the gas channel (bipolar plates) and the catalyst layer, and 4) removing water produced on the catalyst layer. GDLs are usually treated with a hydrophobic agent to improve their performance in delivering reactant gases to the catalyst layer. The hydrophobic agent prevents access of water droplets to portions of the pores in the medium. These portions eventually form hydrophobic passages through the medium. The reactant gases are meant to find their ways through the GDL to the reaction sites through these passages. 5  1.1. Basic operation of PEM fuel cells The interface of the GDL and the catalyst layer is usually treated with a micro porous layer (MPL). The MPL is a mixture of carbon black powder and PTFE particles deposited on the GDL. This layer, which has lower porosity compared to the GDL, reduces the electrical and heat resistivity of the electrodes by providing better contact between the GDL and the catalyst layer, although the porosity of the layer is lower compared to the GDL. Many studies have shown that the MPL significantly improves the performance of the cell at high current densities [PLP08].  1.1.2  Flooding  On the cathode side of a PEM fuel cell, the hydrogen ions (H + ) arrive through the membrane, combine with oxygen and form liquid water (see Figure 1.2). The electrochemical reaction chain occurring in a cell is as follow: 2H2 (g) → 4H + + 4e− (1.1) O2(g) + 4H + + 4e− → 2H2 O(l)  (1.2)  The open circuit potential corresponding to this reaction is 1.223V . Three main losses occur as the circuit is closed. The losses are mainly: 1) activation losses (ηact ), 2) ohmic losses (ηohm ), and 3) mass transport losses (ηmass ). Activation losses (ηact ) are estimated using the Tafel equation which relates the electrochemical reaction to the overpotential [LDF03]: ηact =  kT i ln( ) eα i0  (1.3)  In this equation k is the Boltzmann constant (1.3806×10−23 JK −1 ), T is the absolute temperature (◦ K), e is the electron charge (1.6022 × 10−19 C), and α is the charge transfer coefficient which represents the fraction of the additional energy required for the cathodic reduction reaction at the electrode. In the Tafel equation (1.3), i and i0 are the current density and exchange current density (A/m2 ), respectively. The exchange current density (i0 ) shows the rate of electron transfer between the electrode and the electrolyte. It is a function of catalytic activity of the electrode catalyst material and the concentration of the reactants at the surface of the corresponding electrode. Ohmic losses (ηohm ) are due to electrical resistance Ri (in ohms, Ω)of different components of the cell, while current I (A)is passing through the  6  1.1. Basic operation of PEM fuel cells components. In addition, resistance of the components is one of the major sources of heat generation in the cell [LDF03]: ηohm = I  Ri  (1.4)  i  As the current density of the cell increases, mass transport losses significantly reduce the open circuit voltage. Two main factors contributing to the reduction of the cell voltage are: 1) the reduced concentration of the mole fractions of the reactants (xH2 and xO2 ) on the catalyst layer, and 2) the modified exchange current density (i0 in equation 1.3) affects the kinetics of the reaction on the catalyst, which significantly reduces the cell open circuit potential. The first factor is estimated by the Nernst equation: ηmass,1 =  RT 1 ln( ) nF xH2 x0.5 O2  (1.5)  In this equation R is the ideal gas constant (8.31446J/molK), T is the absolute temperature of the reaction in K, n is the number of ions displaced in the electrochemical reaction, F is the Faraday constant (96, 485C/mol), and xi is the reactant’s mole fraction. The second factor is estimated by ηmass,2 =  cr,i kT ln( 0 ) eα cr  (1.6)  where cr,i0 is the surface concentration of reactant corresponding to the current density of i0 , and cr is the concentration at the catalyst surface (mol/m3 ). Based on these three losses (ηact , ηohm and ηmass ), the cell potential can be determined by Vcell = VOC − ηact − ηohm − ηmass  (1.7)  As the ohmic and mass transport losses are functions of the current density, they can be ignored at low current density and the largest affect on the cell voltage is the activation losses. As the current density increases, ohmic losses become significant. At high current densities, where the rate of water production increases, mass transport losses considerably reduce the cell potential. This is seen as a sudden drop in the voltage and power generated by the cell in a typical polarization curve (Figure 1.3). The polarization curve shows the maximum power obtainable from a typical fuel cell. The power generation monotonically increases with the current 7  1.2. Literature review and theory density. At the onset of a critical value for the current, the power drastically drops due to the increase in mass transport losses. Thus, although the activation and ohmic losses are important and need to be addressed, the major challenge in power gained from a fuel cell is the losses due to mass transport limiting the current output. The first immediate solution to enhance the limiting current is to increase the active area of the cell. This eventually leads to a significant increase in the cost as it requires a larger amount of noble materials such as platinum, as well as larger active area and heavier cell in the stack which will limit the application of the fuel cell in the automotive industry and portable electronic devices. Another solution is to redesign the fuel cell components to transport water more effectively. One of these components includes the gas diffusion layer which interacts with the catalyst layer (where the water is generated) and the flow channel (where the water is removed from the cell). Thus, understanding fluid flow in the porous media of the GDL is crucial to the design of a more effective layer and hence resolving the problem of flooding.  1.2  Literature review and theory  Heat and fluid flow, dispersion and other transport mechanisms in porous media play an important role in engineering sciences including biotechnology, construction and petroleum industry and recently fuel cell technology. For instance, understanding the transport phenomena in porous media has resulted in the development of techniques for enhancing oil recovery [Sha81], filtering and removing groundwater contaminants [BTM09], and implementing porous coating on the surface of replacement joints to promote a stable bone-implant interface [BSH+ 99]. For fuel cell technology, understanding the transport phenomena inside the porous media of the cell is crucial to enhance water management and the efficiency of the cell. The entire network of the fuel cell, except for the gas channels, is considered as porous media with significant structural complexity which can affect fluid flow and displacement within the media. The transport mechanisms occurring in the media can be complex by itself. For instance, heat and mass transfer, change of fluid phase (including evaporation and condensation), and interaction of existing forces in the domain (including viscous, capillary and buoyancy forces) must be considered in the study of flow behavior in porous media. However, a detailed analysis of such processes and the priority of their involvements in the final configuration help in finding the dominant processes in each applications. 8  1.2. Literature review and theory  1.2.1  Modeling of fluid flow in porous media  Two approaches exist for modeling transport phenomena in porous media: 1) discrete, and 2) continuum. Based on conservation of momentum, continuum models deploy a representative volume of the porous medium over which the corresponding equations are averaged and solved with appropriate boundary conditions. Implementation of the averaging technique would result in the extraction of parameters which represent the properties of porous media. Basically, the porosity and permeability of the medium are directly involved when one tries to solve the transport problem inside a porous medium. Once the proper equations are derived, common CFD1 techniques are applied to solve the equations on the domain. However, the complicated configuration of the boundary due to the structure of the medium requires additional techniques to track the boundary even in the steady-state condition. In addition, there would be an open discussion regarding the length scales required to define the representative volume and the time scale required to capture the details of phenomena involved. Similar to the problem of turbulence modeling, the more details expected from the model, the more costly the model will be in terms of calculation and time. Discrete models, on the other hand, split the complexity of the problem of fluid flow in porous media into two major sub-problems: representing the media with an appropriate model, and deploying simple physical rules for solving the flow in the model. As a very simple estimation, a network of pores connected to each other via throats could represent the porous media. The specifications of the network are determined based on the experimental data available for the media under the study. Such data include porosity, permeability, conductivity and tortuosity. Thus, the diameter of pores and throats, the length of throats, the coordination number, the thickness of the network and other parameters are synchronized with the data to develop an optimum network. Once the network is prepared, depending on the phenomenon under study, a constitutional law is employed to model the transport phenomenon. Providing a sufficiently accurate model developed for the medium, there would be no limitation for studying the phenomenon in the medium. In other words, any transport phenomenon such as momentum transfer, diffusion, electrical conduction, heat transfer, and dispersion can be studied. Moreover, any combination of these basic phenomena can also be investigated if the appropriate laws are already known. In essence, pore-network modeling is a platform for studying transport phenomena in1  Computational Fluid Dynamics  9  1.2. Literature review and theory side a porous medium. In continuum modeling, however, the corresponding equation for each phenomenon should be derived independently, as the calculation of the topology of the medium is already considered in deriving the equations. A comprehensive review of the discrete and continuum modeling approaches is provided in by Sahimi [Sah93] and the references therein. Discrete models can be categorized into two main groups: 1) probabilistic models, and 2) deterministic models. The models in the literature as percolation and invasion percolation models are considered as probabilistic models. In these models, the probability of a pore or throat being invaded is determined by the value of the capillary pressure. Thus, they are suitable for very low capillary numbers (Ca2 ). However, as the capillary number increases these models fail to capture the details of the process since viscous force would also be comparable to the capillary force. Percolation models are considered as binary models in which a pore or throat can only have two states: occupied or empty. Thus, the smallest time resolution they can capture is limited to the size of the smallest pore or throat being invaded. In addition, the algorithm seeks the network in each time step to find the pore/throat with higher probability to be invaded. Therefore, in each time step, the invading fluid can only fill one of the pores. This would increase the time of computation drastically without providing comparable details. In addition, percolation–based models are restricted to low Ca numbers as already discussed. In fact, any changes in the boundary condition, either flow rate or pressure, have no effects on the performance of these models; whereas, in practice, the configuration of fluid distribution inside a porous medium highly depends on the flow rate (Ca number) and the mobility ratio (M 3 ). As a concluding remark, probability-based models such as percolation and invasion percolation might be good models to be matched with experimental data, but their application as a tool for studying the phenomena underlying fluid flow inside porous media is highly disputable. The second class of discrete models are referred to as pore-network models. Similar to percolation models, pore-network models also adopt a network of pores and throats as representative of the porous media. However, the rules for displacing the fluids inside the porous media are deterministic. This would ensure that the phenomena taking place at the microscopic scale is actually following basic physical laws rather than statistical laws. Conse2 Capillary number is the ratio of viscous force to the surface/interfacial tension force in the porous medium, i.e. Ca = µV /γ where µ is the phase viscosity, V is the phase velocity, and γ is the surface tension. 3 Mobility ratio, M , is the ratio of viscosity for the displaced and displacing fluids in the porous medium, i.e. M = µ2 /µ1 .  10  1.2. Literature review and theory quently, the predictions of these models are more realistic and can be implemented and interpreted more accurately. While pore-network models take advantage of being discrete models, they also implement and use (in a local sense) the conservation laws from the continuum models. In the capillarycontrolled displacement process, as occurs in the gas diffusion layers of PEM fuel cells, viscous forces are not dominant. Having small Ca numbers, porenetwork models appear to be the proper approach towards the modeling of such displacements [BKS92], [GR93], [BK91], [Blu98], [Blu01].  1.2.2  Network models of fluid flow in Gas Diffusion Layers  Nam and Kaviany [NK03] were the first to implement a pore-network model to gas diffusion layers. They assumed that intersecting fibers surrounded the void volume of the pore body[NK03]. These intersecting fibers overlapped on each so that the solid matrix of the medium was modeled by stacks of these fibrous structures. By changing the position of the fibrous layers, the whole porous medium was reconstructed making the coordination number to be 8. The model enabled the study of formation and distribution of condensed water in the GDL and hence the reduction in medium’s effective diffusivity. In essence, the local water saturation and porosity described the local effective diffusivity of the GDL. The distribution of local effective diffusivity was used along with a model for two-phase flow based on capillary hydrodynamics, and the kinetics of water formation to determine water saturation in the GDL. They also studied the role of the fiber diameter, porosity, and capillary pressure on the cell performance using the network model. These results were then used to develop a two-layer diffusion medium for the enhancement of the cell performance. Using resistor network theory, Gostick et al. [GIFP07] implemented a pore-network model to compute water and gas relative permeabilities, and the effective diffusivity of gas as a function of water saturation. The model was primarily designed to estimate multiphase flow and properties of the GDL. Based on two commonly used GDL materials, the model was calibrated by adjusting its parameters to match the experimental data, especially the absolute permeability as well as the capillary pressure curves for drainage process. By implementing the boundary conditions of a working fuel cell and physical parameters on the network model, limiting current calculations were performed. This limiting current was calculated for different levels of water saturations in a section of the GDL. Half of this section was open to the gas flow channel and the other half was covered by the land. Similar to [NK03], Gostick et al. [GIFP07] calculated the effective diffusiv11  1.2. Literature review and theory ity as a function of porosity and water saturation (s). The relative effective diffusivity and relative permeability were correlated to the saturation with different suggested exponents in the form of sa . No experimental results were reported to verify the above correlation. Sinha and Wang [SW07] used a quasi-static approach to model the immiscible displacement of liquid water inside a GDL. They studied dynamics of water transport at the pore-level and obtained profiles of water saturation in the GDL while the fuel cell was operating in realistic conditions. They reported that the controlling force in transporting liquid water in the GDL is capillary, and water flows through the connected clusters with finger-like frontiers. They also reported the influence of GDL-channel interface covered by liquid water on the transport mechanisms inside the GDL. In their second paper, Sinha and Wang [SW08] considered a GDL with mixed-wettability properties. The effect of wettability distribution on the transport mechanism of liquid water was studied using a pore-network model. The GDL modeled in this study was assumed to be used in a fuel cell at realistic operating conditions. The study showed that for a mixed-wet GDL water transports through the network of hydrophilic pores avoiding the finger-like configuration which occurs in the hydrophobic GDL samples. However, the model did not consider any variable saturation for each pore. As a result, the invaded pore is filled immediately as the water finds a way to get through; whereas, current study shows that pore saturation occurs gradually in progressive time steps. Chapuis et al. [CPQ+ 08] combined pore-network simulations and visualizations obtained from transparent micromodels to study two-phase flow dominated by capillary effects. Using a Voronoi diagram, they developed a pore-network model by assigning particular sections of the diagram to pores, bonds (throats) and disks (which represents the fibrous structure). The paper suggested that the classical invasion percolation algorithm is capable of predicting and simulating the invasion process of liquid water occurring in a hydrophobic medium if the contact angle of hydrophilic portion is much smaller than 90◦ . For contact angles near 90◦ , the mechanisms for the local growth of the interface affect the invasion pattern. For a hydrophilic medium, the simulations suggested a completely different patterns for invasion process. Markicevic and Djilali [MD11] employed a pore-network model with trapping algorithm to study multiphase flow in porous media. They used a 2D regular network with varying diameter for the throat. They also assumed that only one phase can invade a particular throat. The throat radii were determined with respect to experimental measurements obtained for 12  1.2. Literature review and theory permeability (K) and capillary pressure (Pc ). From the model, they determined a length scale (ls ) representing the onset of phase invasion into the next fraction of the network. The length scale (ls ) seems to be a function of flow regimes: for the viscous dominant regime, ls is less than that of the capillary dominant regime. This suggests that the invasion process occurs more frequently in the viscous dominant regime leading to higher saturation and lower relative permeability. For each phase, two different saturations (mobile (si ) and immobile (sim,i )) were defined, and momentum was transferred between phases through the mobile clusters. Thus, phase saturations were defined with respect to the mobility of the phases. For invasion of one phase into the other, the pressure of the invading phase should be higher than the second phase and the capillary pressure of the connecting throat (pc ), which is defined using the Young-Laplace equation. The phase pressure was calculated using the conductance of the throats and the conservation of mass for every pore. The flow in each throat was also calculated using the Poiseuille law. Once the pressure distribution was determined the flow rate was estimated by integrating the flow rates of the inlet and outlet throats. The permeability for each phase was then determined using Darcy’s law. The local capillary pressures were averaged over the termination interface to determine the medium capillary pressure. Using this algorithm, two-phase flow in a gas diffusion layer was investigated. The capillary pressure and relative permeabilities were obtained as functions of the phase saturation, and general flow configuration was studied. Markicevic and Djilali [MD11] suggested that relative permeabilities are power law functions of medium’s saturation. A sensitivity analysis showed that the permeability of the invaded phase changes drastically with the saturation of the invaded phase compared to that of the invading phase. In addition, the numerical results suggested that the exponent is a stronger function of the cluster size than the medium heterogeneity. Thus, the capillary pressure marginally changes by varying the saturations while it significantly varies as the cluster size and the heterogeneity of the medium change. The cluster size also affects the phases’ percolation: phase channeling occurs only for the small clusters. The possible problem of this approach is the application of a 2D network for simulation of two-phase flow in the GDL. Although the assumptions made in the study reasonably predicted the configuration of flow and fluid fronts, they made the algorithm very complicated and difficult to comprehend. In addition, there is lack of physical evidence to validate the assumptions used in the algorithm (e.g. the length scale (ls ) being a function of the flow regime). Kuttanikkad et al. [KPP11] studied the effect of hydrophobicity on water 13  1.2. Literature review and theory flow in the GDL using a pore-network model. The study used a regular 3D network of cubical pores interconnected with ducts representing the throats. An average pore diameter of dp = 17.5µm on a 40 × 40 × 10 network was distributed. All the pores and throats were also considered as elements which can be either hydrophobic or hydrophilic. In order to incorporate the effect of hydrophobicity, a portion of the pores and throats (1 − f ) was considered as a hydrophobic portion by assigning a contact angle of 115◦ . For the remaining fraction (f ), the contact angle was assumed to be around 80◦ , representing the hydrophilic property of the GDL sample. Thus, each element (i.e., pore or throat) is attributed with a threshold pressure which is a function of the element geometry (radius of the throat and pore), interfacial properties of the element (contact angle (θ) and surface tension (σ)) and an adjusting pressure determined with respect to hydrophobicity of the medium to reflect the experimental measurement [FCS10]. The pore radii (rn ) used to evaluate the pore threshold pressure were calculated with respect to the initial pore size distribution and the status of the surrounding throats. For instance, the pores are larger if the surrounding throats are already invaded and smaller if the throats are not wetted yet. This would increase the chance of invasion of the pore if the surrounding throats are wet. The model was then implemented to different hydrophilic fractions (f ) and showed that the capillary pressure-saturation curves are identical for these cases when the liquid saturation is low (less than 50%). For higher saturation, the curves reproduce the results obtained experimentally [FCS10]. In order to estimate relative permeabilities, arbitrary pressure values were imposed as inlet and outlet boundary conditions of the network (Pinlet and Poutlet ). Then, the total flow rate (Q) over the entire network was calculated using the porenetwork model. Then, Darcy’s law was used to determine the permeability. The results show that at low liquid saturation, the relative permeabilities are almost identical for all the hydrophilic fractions (f ). The relative permeability of the gaseous phase slightly increases as the hydrophilic fraction increases while in contrary the relative permeability of liquid phase decreases with hydrophilic fraction. For GDL samples at low saturations (lower than the percolation threshold, fc ), they concluded that the macroscopic properties of the medium such as permeability, capillary pressure and saturation are independent of the hydrophobicity of the medium. El Hannach et al. [EHPP11] compared two different pore-network algorithms to study transport of liquid water in the cathode catalyst layer (CCL) of a PEM fuel cell. The study assumed that the platinum particles and carbon support in the form of spherical agglomerates represented the CCL. Water was assumed to be produced within these agglomerates. In 14  1.2. Literature review and theory the pore-network model, the void spaces in between the pores and throats were considered to be filled with the agglomerates. Water was injected into the network through the randomly distributed “active” agglomerates within the network. All the network components (pores and throats) were evenly open to water injection. The criteria for invasion followed the same procedure described in [KPP11], i.e. the pore diameters were enlarged if the adjacent pores/throats were filled with water as proposed in [Blu97],[Blu98]. The study shows that liquid water configuration follows the fingering regime and stable displacement (“fist-like” configuration) in the hydrophobic and hydrophilic media, respectively. Thus, for the hydrophobic medium the fingers grow independently making the total saturation higher compared to the hydrophilic medium. It also shows that water enters the gas diffusion layer through independent injection points, representing the breakthrough points of the CCL. In contrast to the study performed by the same group [KPP11] which concluded that the effect of GDL hydrophobicity is negligible at low saturation, El Hannach et al. [EHPP11] showed that wettability has significant effect on liquid water transport in the catalyst layer. Recently, Pauchet et al. [PPSK11] studied the effect of GDL hydrophobicity degradation on the performance of PEM fuel cells using a porenetwork model. They argued that the loss of hydrophobicity in the GDL can explain the performance loss during the life span of the PEM fuel cell. Similar to the algorithm implemented in [KPP11], they simulated the loss of hydrophobicity by increasing the fraction, f , of the hydrophilic pores. As the hydrophilic fraction (f ) value approaches the percolation threshold (fc ), value the effective gas diffusion coefficient drops significantly which induces a decrease in the cell potential. Kang et al. [KKNK11] developed a “similarity model experiment” to simulate the flow of a non-wetting fluid through a porous medium with the same dimensionless numbers as in a working PEM fuel cell. A physical domain of 30cm × 30cm × 7.2cm was modeled with uniform flux injection through a disk with diameter of 10cm. The invasion-percolation path finding procedure developed by Lee et al. [LNK10] was used to determine the saturation distribution of the non-wetting fluid at steady-state condition. The algorithm sought the transport pathways with the lowest values for entry pressure. The procedure starts from the pores at the inlet which were occupied by the non-wetting phase. The results of the visualization showed that an intermediate layer sandwiched between a fine and a coarse layer delays the process of capillary fingering in the coarse layer, which reduces the volume of the non-wetting fluid. The pore-network developed in this study measured the total volume of non-wetting fluid injected into 15  1.2. Literature review and theory the porous medium. The modeling results are in close agreement with the experimental results. Wu et al. [WLZW12a] developed a pore-network model based on an invasion-percolation algorithm to study transport of liquid water and reactant gas transport through the GDL. The average diameter and length of the throat in this model are considered to be 9µm and 25µm, respectively, within a 80 × 80 × 12 network size. The network was initially considered to be completely hydrophilic, and then a fraction (1 − f ) of the pores were rendered into hydrophobic pores by randomly assigning a constant contact angle of θ = 110◦ . The invasion process followed a sequential procedure, in which a number of the inlet throats are the injection points. In this procedure, liquid is injected into the medium from the first injection point till the breakthrough is achieved. The process is then repeated for the second injection point and so on. The study concludes that the liquid water is evenly distributed through the GDL thickness for low hydrophilic fractions (f < 0.2). It is also shown that the hydrophobicity does not have any influence on the limiting current of the cell when the number of injection points is low (which simulates the low current density condition). At higher current densities, however, the simulation shows that the limiting current is not changing monotonically with hydrophobicity, and hence the performance of the fuel cell is maximized at an optimal loading of PTFE. The same group investigated the transport and reaction process in the cathode catalyst layer (CCL) using a pore-network model [WLZW12b]. The catalyst layer was assumed to have primary pores (as the agglomerates) and secondary pores (as the void space within the catalyst layer). Thus, the model only considers the secondary pores to simulate the reaction and transport phenomena. The length and average diameter of the throats were assumed to be 100nm and 55.7µm, respectively. A 15 × 15 × 100 network was constructed corresponding to the 10µm-thick catalyst layer. The oxygen concentration and the reaction rate profiles along the CCL thickness were predicted using the model.  1.2.3  Experimental techniques  Validation of the models developed to predict fluid flow in porous media is essential in order to make the models reliable. Experimental results are used for this purpose. The experimental results include measurement of macroscopic properties of the medium such as porosity, permeability, diffusivity and effective diffusivity, capillary pressure, saturation and etc ([Whi67], [GFI+ 06], [NK03], [PMS12], [MBD07]). Once experimental mea16  1.2. Literature review and theory surements are available, models are modified to match the measurement results. A reliable model predicts the characteristics of flow as close as possible. In addition, being able to predict the flow behavior for different conditions is the essence of a reliable model. For the porous media involved in PEM fuel cell applications, direct visualization of fluid flow is an effective tool to study the behavior of such flows. Thanks to the very thin nature of the media and very slow flow of fluids involved, the phenomena can be directly observed. Considering the thickness of porous material used as the GDL (100 − 400µm) and its high porosity (70% < < 90%) the details of fluid flow through the medium can be captured with high spatial and temporal resolution. The captured images are then analyzed to develop an understanding of phenomena occuring through the medium. In addition, the data obtained from direct visualization can be used to validate models developed for fluid flow in porous media. Different techniques including magnetic resonance imaging (MRI), neutron imaging, X-ray microtomography, direct optical visualization, and fluorescence microscopy have been employed to study water formation and distribution in different components of the fuel cell. Magnetic Resonance Imaging (MRI) MRI technique is a non-destructive method widely used in medical applications. The technique uses a radio-frequency signal to excite particular nuclei in the presence of a strong static magnetic field. The excited nuclei absorb the signal and resonate accordingly. The measurement of emitted resonance determines the existence of nuclei at a particular location. The MRI technique has been used for liquid water detection and determination of diffusion coefficient in the membrane [ZTD+ 93], as well as the measurement of the distribution of water content in the membrane during the operation of a fuel cell [TTH04]. These measurements have shown that the water concentration gradient in the membrane and the overall water content decrease as the cell current increases. Also, MRI measurements indicated that the rapid drop in the cell voltage, which occurs within the first three minutes of the experiment, is associated with a decrease in water content on the anode side of the membrane [TTH04]. Effects of the membrane thickness and cell working current on the level of hydration and distribution of water in the membrane have also been studied using MRI techniques [TTH05]. It was shown that a thick membrane is better hydrated than a thin one and the total water content in the membrane decreases for higher current densities [TTH05]. The MRI technique has also 17  1.2. Literature review and theory been used to study the effect of diffusion in water transport between the catalyst and membrane [FLS+ 04]. The results reveal that radial gradient diffusion removes water from the catalysts into the surrounding membrane. Another use of MRI has been in the study of the effect of the gas flow field configuration on the distribution of water in the membrane and cathode flow field [FBW06]. The results have shown a uniform water distribution in the electrolyte due to counter-flow configuration. In essence, the maximum power output was obtained when liquid water is first visible in the MRI image of the cathode flow field, and subsequently the power decreases as the liquid water continues to accumulate [FBW06]. Neutron Imaging The neutron imaging technique is another method which has been used to visualize water inside the network of a fuel cell. In this technique, the material under study is bombarded with beams of neutrons. If water exists in the material the intensity of the neutron beam is impaired. This impairment is proportional to water content in the material. The technique can be employed while the fuel cell is working in real conditions, and unlike the MRI technique, neutron imaging is not sensitive to the material used in the fuel cell. However, the cost of this method has limited the use of this technique compared to other visualization techniques [Baz09]. Neutron imaging has been used to study different effects including water distribution along the thickness of the membrane [MGP96]; production, transport and removal of water throughout the cell components [SJAW04]; and visualization of two-phase flow in an operating PEM fuel cell [PHC+ 05]. Using neutron imaging, Pekula et al.[PHC+ 05] have found that at low power density the anode channels are blocked; while at high power density water is distributed more evenly. Trabold et al. [TOJ+ 06] used the neutron imaging technique to study the effect of the current density and the ratio of cathode stoichiometry on the volume of water. They reported a decrease in water content of the cell as the load of the cell increases indicating a higher gas velocity is required for removing water from the flow channels when the cell is working at the current density of 1.0 A/cm2 . In addition, the effect of the stoichiometric ratio of the cathode on the amount of water content is reported to be negligible. They concluded that the voltage loss at high current densities is the result of accumulation of liquid water in a cell’s porous material including the GDL and MEA as opposed to the impediments in water transport in the channel. Owejan et al. [OTJ+ 06] employed the neutron imaging technique to 18  1.2. Literature review and theory measure water transport through the GDL of a PEM fuel cell. They used an interdigitated flow field in order to conduct the reactant into the porous diffusion layers. The pressure differential from the inlet to the outlet was correlated to the level of saturation in the GDL. The study concluded that flooding primarily occurs in manifold at the outlet of the flow field. In addition, they reported a significant reduction in relative permeability as the water content in the GDL reaches the critical water mass. Thus, mass transport is the limiting factor when liquid water fills half of the void space in the cathode diffusion layer. The neutron imaging technique was used by Turhan et al. [THBM06] to study the effects of the reactant gas flow rates, operating cell pressures, and inlet relative humidity on the cell performance. The study showed that increasing the inlet gas flow rate significantly reduces the accumulation of liquid water in the cell. The results also suggested opposite trends for the effect of the over- and under-humidified inlet condition on the cell’s water content. For the over-humidified inlet condition, decreasing the pressure of the cell increases the amount of liquid water; while for the under-humidified inlet condition the trend is reversed. Zhang et al. [ZKS+ 06] studied the impacts of GDLs and the flow fields material on the formation of liquid water and its transport using the neutron imaging technique. For a serpentine flow field, they showed that the hydrophobicity of the material of flow field affects the presence of liquid droplets. They also studied the effect of the GDL material on accumulation of water on the cathode side of a quasi one-dimensional cell. It was shown that the cloth–based material generally held less water compared to the paper–based material. They reported that a significant content of water in the diffusion layer does not limit the performance of the cell. In contrast, Siegel et al. [SMS+ 08] reported significant voltage drop due to the accumulation of liquid water on the anode side of a PEM fuel cell. They used a cell with dead-end anode and deployed neutron imaging technique to observe water accumulation. For a range of temperature, cathode inlet RH, current density and air stoichiometric ratio, they found the rate of water accumulation. The results confirmed the importance of anode flooding even when the cathode is not plugged with liquid water. Bellows et al. [BLA+ 99] deployed neutron imaging in order to determine the gradient of water content inside the Nafion membrane. They addressed the experimental difficulties in neutron imaging by using thick membranes and careful alignment of the cell within the neutron beam. Semi-quantitative interpretation of their results suggests that the diffusivity of water in the membrane is higher than the values predicted by the models [ZTD+ 93], 19  1.2. Literature review and theory [SZG91]. Using the neutron imaging technique, Hickner et al. [HSC+ 06] reported a lag of 100 seconds between the water content of the PEM fuel cell and the current density. The study also reported less water content within the cell when the working temperature is higher. For the temperature of 60◦ C the study showed a peak in the water content for the current density of 650mA/cm2 . Further increase in the current density was associated with the lower water content. Owejan et al. [OTJ+ 07] used the neutron imaging technique to obtain two-dimensional distributions of water in a cell with an active area of 50cm2 . They reported that more water exists in the hydrophobic-coated flow field channels. However, the water slugs formed within the channel are smaller which improve the performance of the cell at high current densities. For some cases, they reported a significant difference in the cell performance when the accumulation of liquid water marginally increases. Thus, flooding of the GDL was recognized as the primary reason for the sudden voltage drop. Kim and Mench [KM09] used neutron imaging to investigate the phase-change-induced water transport mechanism in a PEM fuel cell. They recognized the gas phase as a key controlling parameter in this transport mechanism, i.e. if the gas phase exists in the catalyst layer or gas diffusion layer the phase-change-induced water transport is the dominant mechanism and the net water flux is observed from the hot to the cold region. X-ray microtomography The principal of X-ray microtomography, also referred to as micro-computed tomography (µCT ), is similar to that of the neutron imaging technique, but the static magnetic field is not necessary. The material under study is scanned with X-rays and, depending on the material, the X-ray beam is attenuated. An array of detectors receives the transmitted radiation from the material. Subsequently the material specifications are determined by analyzing the transmitted signals. Sinha et al. [SHW06] used this method to investigate the water saturation distribution through the GDL. They obtained the distributions along the thickness of the GDL as a function of time during the gas purging step and showed an exponential decrease in the drying rate with purge time. Manke et al. [MHG+ 07] used X-ray radiography to study the evolution and transport of water inside the GDL in an in-situ experiment. They reported the formation, growth and transport of liquid water and then correlated the operating conditions to the dynamics of droplet formation. Their results suggested that the pores of the GDL are filled continuously which demon20  1.2. Literature review and theory strates a continuous “capillary-tree-like” transport mechanism which is in agreement with the results reported by Pasaogullari and Wang [PW05]. This eruptive transport process observed in the experiment, also directed them toward the implementation of capillary tree patterns in their model at high current densities. Lee et al. [LLK+ 08] used X-ray radiography to determine quantitatively water distribution between the flow field and the GDL. They also employed an image processing technique to remove noise from the images obtained from the radiography. A very thin cross section of the cell including membrane sandwiched between two GDLs was investigated by Hartnig et al. [HMK+ 08]. The spatial and temporal resolutions in this study were 3mm and 5s, respectively, which is quite high resolution compared to other X-ray imaging employed to investigate water transport in fuel cells. The study claims that liquid water droplets cause two different diffusion barriers in the GDL when the cell is working at high current densities and surface properties of the GDL, such as hydrophobicity, determine the location of the barriers. Direct visualization techniques The direct methods to observe flooding and water transfer in a PEM fuel cell usually requires a cell with transparent components. For instance, the formation of water on the surface of the GDL as well as at the interface of the cathode gas channel has been observed using direct visualization [TPH03],[HMWH04],[YZLW04],[BKTO05],[SNY+ 05],[GW07],[KSM06]. A microscope is usually installed and a cell with transparent flow field is used. The level of insight that these methods can provide is limited due to the opaque nature of the cell components. However, the observation of water droplet formation on the GDL surface and flow field channels is correlated to the performance of the cell if the voltage and current of the cell are measured during the experiment [LMW06],[SPA07],[TOG+ 06],[KWKB08],[WSHL06], [LKR+ 09]. Tuber et al. [TPH03] studied cathode flooding in a small fuel cell using direct visualization. The voltage discharge performance of the cell was measured while the images of the cathode were being acquired during the operation of the cell. They reported cathode flow field blockage as a result of flooding, which impaired the performance of the cell significantly. The study reported a drop in the current density of the cell for standard and hydrophobic GDLs; while the hydrophilic GDLs showed a constant value over the first 40 minutes of operation. Although the paper argued that better 21  1.2. Literature review and theory performance of the fuel cell was achieved with treated GDLs, the experiment was performed for a limited time intervals. This prevents a concrete conclusion regarding the effects of hydrophobicity on the performance of GDLs in transferring water. Interestingly, the authors mentioned limited diffusivity due to more uniform distribution of water inside the hydrophilic GDL. An experimental model of an operating PEM fuel cell was developed by Borrelli et al. [BKTO05] to study the flow of water in the cathode channel. They qualitatively compared the similarities in behavior of water droplets between the model and the experimental observations reported by Tuber et al. [TPH03]. Hakenjos et al. [HMWH04] designed an experimental test cell to measure simultaneously the current, temperature and water distribution. They developed a segmented anode flow field for measuring the current. An optical window was located on the back plate of the cathode flow channel allowing direct visualization. They used three different rates of air flow and reported its influence on the temperature and water formation and distribution. For the smallest air flow rate, the study showed that a large area of flow field is flooded with water at a low current of 3.8A. The distribution of current showed that only active area in the inlet. No water is visible in the channel at the highest air flow rate. The distribution of current confirms existence of regions with high current densities near the cell center and the temperature distribution is almost uniform over the entire cell area. Yang et al. [YZLW04] used a cell with a transparent cathode and explored liquid water transport at high current densities. They reported nonuniform water distribution between the channels. The emergence of liquid water droplet from preferential openings on the GDLs were studied and discussed. They highlighted the contribution of the hydrophilic flow channel in removal of liquid droplets from the hydrophobic surface of the GDL. The results identified liquid film drainage from the gas channel as an important factor to avoid flooding in PEM fuel cells. Sugiura et al. [SNY+ 05] separated the gas water flow path by implementing a water absorption layer. They adopted the direct visualization technique to compare the gas flow characteristics of their design with the conventional cathode gas channels. The study showed a reduction in the flooding rate by a factor of 4 using the absorption layer. However, the polarization curve obtained for their design including the absorption layer showed a higher voltage drop at high current densities. Ge and Wang [GW07] adopted the direct visualization technique to investigate the formation of water on the anode side. They used both serpentine and parallel patterns as the flow channels. The study showed that 22  1.2. Literature review and theory the cathode GDL contained water droplets while the anode GDL was completely dry. They reported a significant effect of GDL hydrophobicity on water distribution on the anode side for the current densities smaller than 0.2A/cm2 . For the hydrophobic GDL, the water is condensed on the flow channel blocking the passage of reactant gases, whereas the hydrophilic GDL absorbs the liquid and prevents channel clogging. Although the argument is valid for transient conditions, higher saturation of GDLs is proved to reduce considerably the performance of the fuel cells. Kumbur et al. [KSM06] developed an experimental model to observe and study the removal of water droplets from the surface of the GDL using a direct visualization technique. This study was similar to the study conducted by Borrelli et al. [BKTO05]. They primarily studied the effect of the droplet shape and size and flow characteristics on the rate of water removal from the channels. The results showed that at the high flow rate regime (Re < 600), the effect of PTFE loading on the hysteresis of contact angle is more important; while at the low flow rate on the anode side, the surface hydrophobicity has negligible effect on the instability of the droplet and hence on droplet removal. Spernjak et al. [SPA07] examined the effectiveness of various GDL materials in water removal from the cathode flow field using direct visualization. They visually studied the effect of the microporous layer transport of liquid water in the anode channel. They showed the MPL of the cathode side is responsible for the formation of liquid water droplet on the anode side. They argued that the pressure gradient produced by the existence of the MPL on the cathode side diffuses the water droplet into the membrane and towards the anode side. They concluded that untreated GDLs cannot provide enough water to keep the membrane hydrated. In addition, the pores saturated with liquid water limit the transfer of reactant gases. They identified a narrow contact area between the sidewalls of the channel and the GDL as the primary area where water droplets emerge out of the GDL. In contrast, the GDLs treated with the hydrophobic agent convert liquid water into individual droplets and distribute the droplets over the entire interface of GDL and channel. This would enhance availability of the pores for transport of reactant gaseous. Using a transparent test cell, Theodorakakos et al. [TOG+ 06] studied droplet formation and detachment in the cathode channel. They confirmed that single-droplet detachment from the GDL surface is the main mechanism for water removal. They used the experimental results of direct visualization to refine and improve a numerical model developed for droplet detachment. Their study showed that an increase in the temperature helps the detach23  1.2. Literature review and theory ment of liquid water droplet primarily due to reduction in liquid surface tension. The same method was used by Weng et al. [WSHL06] to investigate the influence of the cathode gas concentration and humidification on the performance of a PEM fuel cell. They showed that higher stoichiometric values improve the performance of the cell under humidification. For the case of non-humidified conditions, membrane dehydration occurs and the cell response is not steady. They also concluded that for low stoichiometric values the change in humidification does not have a significant effect on the performance of the cell. Kimball et al. [KWKB08] studied the effect of different orientation of flow channels on the performance of the cell using direct visualization. They claimed that gravity is important and has effect on the water distribution along the channel and the regime of liquid water flow inside the channels. The results suggested that the invasion of liquid water into the GDL requires a minimum hydrostatic pressure. Following the initial invasion, water fills the largest available pores in the GDL and travels toward cathode gas channel. Water droplets formed on the surface of the GDL are connected to the liquid within the GDL pores. Once the shear force in the channel is enough, the droplet is detached and removed from the channel in the form of liquid slugs. Another type of visualization method is the fluorescence microscopy technique in which the pathways of water droplets inside the porous medium are determined by tracing the fluorescent dyes of the water solution. The solution is injected into the GDL sample. A UV light source excites the fluorescent dyes in the solution. The excited dyes start emitting light in a longer wavelength and the wavelengths are captured by a camera. Thus, the flow of liquid water inside the GDL can be captured while limited by the opaque nature of fiber structure. Litester et al. [LSD06] were the group first that applied the fluorescence microscopy to liquid water transport inside the GDL. They captured the images in consecutive time steps. The transient image intensity data was correlated to the height of the liquid water in the medium. The transport of liquid water through the GDL was determined using these high spatial resolution images. The observations of this study suggested that channeling and fingering mechanisms are dominating flow of water rather than the capillary tree mechanism mentioned in prior work and models [NK03], [PW05]. Bazylak et al. [BSLD07] employed fluorescence microscopy to study the effect of GDL compression on the formation of preferential pathways for water transport. They concluded the change in hydrophobicity due to compression of the GDL guided the liquid water to travel through the compressed 24  1.3. Hypothesis area. Alternatively, compression may cause damage to the porous structure of the medium which results in the formation of macroscopic cracks on the sample. These cracks may be chosen by water as the preferential pathways since their corresponding radius are much higher than the actual pore size of the medium. Thus, water would face less resistance and flow through the compressed region. Using scanning electron microscopy (SEM) the morphology of the samples was also studied after compression. The breakup of fibers was reported as well as deterioration of the hydrophobic coating. In their next work, Bazylak et al. [BSD08] studied dynamic water transport and emergence of droplets through the GDL. The study showed that droplets appear at breakthrough locations which change periodically. Thus, pathways of liquid water in the GDL are dynamically interconnected. The same group used fluorescence microscopy technique to study the interactions of the liquid water with a solid wall on a PTFE-treated GDL [BHDS08]. The effect of wetting properties of the plate on the stability of a droplet was investigated. The study proposed the use of hydrophobic land areas (at the GDL/land interface) as the channel hydrophobicity enhances the mitigation of the accumulated water droplet under the land area to the flow channel.  1.3  Hypothesis  In the Literature review and theory section, the network models of fluid flow in Gas Diffusion Layers were presented followed by the experimental techniques for visualization of fluid flow in porous components, as well as measurement of hydrodynamic properties of these components. In the last section of Experimental Techniques, direct visualization techniques were presented. These direct techniques are categorized as 1) in–situ experiments and 2) ex–situ visualizations. In the in–situ experiments, a fuel cell with transparent gas flow channels working at a range of current densities, is deployed to correlate the electrochemical responses of the cell (i.e. voltage variations and power output) to the hydrodynamic phenomena (e.g. qualitative flow patterns, evaporation/condensation rate, etc.) and hydrodynamic properties (e.g. saturation-capillary pressure curves), occurring in porous components (see [TPH03, HMWH04, YZLW04, BKTO05, SNY+ 05, GW07, KSM06, SPA07, WSHL06, KWKB08]). The advantage of the in–situ measurements is that they can primarily correlate impaired performance of the fuel cell (low level of voltage/power at high current densities) to the hydrodynamic properties of the porous media, i.e. these methods are successful in providing evidence of voltage drop as a direct result of flooding. How25  1.3. Hypothesis ever, these measurements fail to provide extensive details of the phenomena causing the weak performance of the cell on a microscopic level. The ex–situ visualization techniques, on the other hand, concentrate at the micro-scale investigation of fluid flow in porous media by isolating the porous layers, and by directly simulating the flooding condition in these layers. In this way, the complexity of the experiment is significantly reduced by removing the electrochemical components of the fuel cell operation (e.g. chemical reactions of the anode’s and cathode’s catalyst layers, proton exchange transport phenomenon of the electrolyte membrane, and evaporation/condensation occurring in the gas flow channels), while the effects of these electrochemical phenomena are reproduced in the experiment by direct injection of water to the GDLs (see [LSD06], [PW05], [BSLD07], [BHDS08], [BSD08]). This research uses an ex–situ visualization technique (fluorescence microscopy) which was originally used by Litster et al. [LSD06], and followed by Bazylak et al. ([BSLD07], [BHDS08], [BSD08]) to investigate water flow through the Gas Diffusion Layer of a PEM fuel cell, by isolating this porous layer (GDL) from the rest of the fuel cell and directly injecting fluorescent solution into it. The acquired images in the previous studies were processed by correlating the image intensity to the height of fluorescent liquid column in the GDL sample. In this study, however, the image intensity is correlated to the saturation of fluorescent solution in the sample. Thus, highly saturated and low saturated regions are distinguished using this method. By processing the experimental images, the patterns of fluid flow are recognized and presented in three categories: 1) initial invasion pattern, 2) progression pattern, and 3) pore-filling pattern. Based on these patterns, a pore-network algorithm is developed featuring a pressure-correction term, which correlates the pore saturation to the total pressure of the pore. Another major shortcoming of previous studies ([LSD06], [BSLD07], [BHDS08], [BSD08]) is that the flow rate, with which all the experiments were conducted, was significantly higher (usually 1000 times more; note that the flow rates of these experiments had never been directly reported but the order of magnitude can be estimated by comparing the breakthrough times provided in these studies with the breakthrough times measured in this research) than the actual flow rate corresponds to the real working condition of the fuel cell or the flooding condition. Since flow phenomena in porous media is a function of capillary number (Ca = µV /γ), a high flow rate alters the capillary number by several orders of magnitude which results in a totally different flow regime, and consequently different characteristics of flow. In this study, however, the flow rate associated with the flooding 26  1.4. Thesis organization is accurately computed using the hydro-electrochemical correlations, and all the experiments are conducted using this flow rate. It is worth mentioning that handling the experimental setup with such an ultra-slow flow rate (≈ 10−11 m3 /s) is an extremely difficult task and it requires meticulous methods and well-thought experimental procedures prior to the execution of measurements. In the light of this hypothesis, thesis organization is presented in the following section.  1.4  Thesis organization  The thesis is organized as follows: A thorough review of relevant literature is presented in the first chapter, Introduction. The statement of the research question and its goals are presented in the same chapter. The next chapter, Experimental measurements, describes the visualization technique used to study the formation and flow of liquid water in the GDL. The experimental setup, the process of sample preparation, image acquisition, calibration, measurements and data collection are presented in detail. Then, the results of the visualizations and measurements are presented, analyzed and discussed. The analysis of the patterns of liquid water flow through the GDL leads to the development of a model for this flow. In Chapter 3, Numerical model a pore–network model developed to study the flow of liquid water within the GDL is presented. The numerical scheme, including network generation and treatment, and the displacing algorithm are presented and discussed. The model is then qualitatively validated. In Chapter 4 , Model implementation, the pore–network model is applied to the experimental conditions and measurements obtained in Chapter 2 and validated against these results. In the last chapter, Conclusions, an overall analysis and integration of the research is presented. Conclusions regarding hypothesis of the dissertation are presented and discussed. The strengths and limitations of this research are also presented in the same chapter. The potential applications of the experimental method and the numerical scheme are introduced. The chapter ends with possible directions and aspects of future research.  27  Chapter 2  Experimental measurements Ex-situ experiments are designed to study a specific component of a whole system while the component is isolated. The real working conditions of the system are reproduced and set as the boundary conditions for the ex-situ measurements. Thus, the effects of the components properties and the boundary conditions on the phenomena occurring inside the component, and its response are studied. In this research, the gas diffusion layer (GDL) of a PEM fuel cell is chosen as the component under study, due to its critical role in flooding and water management. The GDL is isolated and flooding conditions are applied as the boundary conditions of the experiment. The interactions of the GDL with the cell are: 1) the GDL contact with the catalyst layer, and 2) the GDL contact with the gas flow channel. At the interface of the GDL and the catalyst layer, water is generated due to the chemical reactions. A portion of the produced water diffuses back into the electrolyte membrane (PEM), while the rest diffuses into the the GDL. In ex–situ measurements, this diffusion is simulated by directly injecting liquid water into the GDL sample. The interface of the gas channel and GDL is simply simulated by having the other surface of the GDL open to the air, which is a valid assumption for air breathing PEM fuel cells working in real conditions. In essence, the response of the GDL is studied with all the physical interactions of the GDL and other cell components are simulated. The pattern of liquid water moving through the GDL porous medium, the breakthrough time and the breakthrough pressure are the main quantities measured in this study. For the former, a fluorescence microscopy technique is adopted. In this technique, fluorescent dyes are dissolved in liquid water, and are traced during the injection of water into the the porous sample. During the injection process, sequence of images are captured while water is flowing through the medium. The images are used to estimate the breakthrough time. A pressure transducer is installed close to the sample registers the system pressure during the experiment and at the breakthrough time. The imaging results are used to understand flow phenomena happening inside the complex structure of the GDL. These results are useful for proposing the pore-network model presented in the numerical chapter (See Chapter 3). 28  2.1. Setup The pressure measurements are used as the boundary conditions in the simulation and the breakthrough time is the vital quantity for validation of the numerical results. The application of the fluorescence microscopy for ex–situ studying of fluid flow in GDLs was originally proposed by Litster et al. [LSD06] and followed by Bazylak et al. ([BSLD07], [BHDS08]). These studies solely used the fluorescence vertical illuminator to visualize water evolution through the GDL, without any side measurements, while in this study the same imaging setup is combined with 1) voltage sensors to estimate the flow rate prior to water injection, and 2) a pressure transducer to measure the variations of setup pressure during the injection and at the breakthrough. The measured data are extensively used as boundary condition and validation criteria, as is explained in detail in Chapter 3. This chapter is organized as follows: In the Setup section, an overview of the experimental setup and its components is presented, including the injection module, the imaging module and the registration module. The details of water injection procedure, pump and tubings and the calculation of flow rate are discussed. The process of acquiring the images during the experiment, measuring the time of breakthrough and the breakthrough pressures are also provided. The section also includes the calibration process used to calibrate the time of due to the expansion of the setup components. In the Results and discussion section the images of water invading the pores of the GDL are presented for both hydrophobic and hydrophilic samples. The images are throughly analyzed and the effects of hydrophobicity on the pattern of flow is discussed. Finally, the measurements for the breakthrough time and pressure are presented and discussed.  2.1  Setup  Three main modules of the experimental setup are presented in Figure 2.1. The injection module (including the syringe pump, tubing and the injection channel) introduces the fluorescent solution into the porous medium of the GDL. The image acquisition module consists of a fluorescence vertical illuminator and the associated software rendering and saving the captured images. The imaging section of the experimental setup (including the injection module and the image acquisitions module) is essentially the same as the original method proposed by Litster et al. [LSD06]. However, the registration module is unique to this particular setup. The components of 29  2.1. Setup the registration module are micro-needle sensors, the pressure transducer, the sensor terminal and a lab PC which records the data obtained from the sensor terminal. For each experiment, a fresh GDL sample is placed in the injection channel. The fluorescent solution is pumped into the GDL sample from the bottom surface. The image acquisition module captures the images of water flow through the GDL in incremental time steps. The registration module collects the data including the images and signals sent from the camera, micro-needles and the pressure transducer. The experimental setup was designed and developed from the scratch by the author in the UBC’s Advanced Thermofluidics Laboratory for the purpose of this research. Each module of the setup is explained in detail in the following sections.  2.1.1  Fluorescence microscopy  As explained above, a fluorescence microscopy technique is used to visualize water flow inside the GDL. In this technique, the specimen is illuminated with a UV light of short wavelength. The much weaker emitted fluorescence with longer wavelength is then separated from the excitation light. The emitted light reaches the detector (eye) in a properly configured microscope setup in a way that the obtained fluorescence structures are superimposed contrasting against the black background ([PM08]). Figure 2.2 shows a schematic of the fluorescence microscope. The UV light source (EXFO X-Cite120Q Mercury halide short arc), provides an ultraviolet light beam with a wavelength in the range of 400nm to 600nm. The light wave passes through an excitation filter inside the illuminator to filter out wavelengths other than 490nm. Reflected from the dichroic mirror, the short wavelength light baths the GDL sample and excites the fluorescence dyes in the solution. The emitted light passes through the dichroic mirror and a barrier filter (550nm) and is captured by the CCD camera (DFC 340 FX). The vertical illuminator is basically an APO Zoom microscope with a photo tube (HC L 2TU 1.25X). The magnification of the main objective lens 1.0X and the C-Mount adapter also provides 0.55X zoom. The field of view for each image is a window of 3.3mm × 3.3mm. The total magnification is 25X and numerical aperture (NA) is 0.049 making the depth of field of 0.392mm. This assures the entire GDL thickness to be in focus, as the liquid is flowing inside and images are captured. The lateral resolution is 6.8µm and the time interval between two consecutive images is 1s. Since the fluorescence dyes are significantly shiny compared to 30  2.1. Setup the black background of the image (the fibrous material of the GDL), the time of exposure must be optimized. For a long exposure time, the image will be bleached out; whereas for a short exposure time the low intensity spots (presenting locations with the lower solution content) might not be appropriately captured. The optimum exposure time for the experiments was set as 10.9ms.  2.1.2  Injection channel  The injection channel is fabricated from a polytetrafluroethylene (PTFE) rod (ASTM D 1710) which is hydrophobic. The hydrophobicity of the PTFE channel prevents any unwanted leakage in the gap between the GDL sample and the surface of the injection channel. This ensured that all the solution delivered by the pump is injected into the GDL sample. The detailed schematic of the injection channel is shown in Figure 2.3. The multiport injection channel features unique characteristics which make it distinguished from the injection modules used in previous studies ([LSD06], [BSLD07], [BHDS08], and [BSD08]). This includes installation of micro-needles serving as voltage sensors, a pressure transducer, and a very thin layer of ultra-hydrophobic material on its surface. The channel is composed of a vertical hole with a diameter of 3mm drilled through the PTFE rod, and four horizontal holes with a diameter of 300µm. Four micro–needles are inserted into the horizontal holes. The first needle applies voltage to the solution passing through the channel. The other needles are attached to the voltage sensors. Knowing the distance between the needles and the time sequence of the signals provided by the sensors, the flow rate can be estimated. The estimated flow rate is then compared to the flow rate provided by the syringe pump to detect any possible leakage in the injection module (tubing, syringe and the channel). The surface of the channel is treated with a very thin layer of Polydimethylsiloxane (PDMS). In addition to the extra hydrophobicity that PDMS provides, its flexibility prevents the GDL sample from breaking under the excessive pressure imposed by the PTFE cylinder and the channel surface. However, the PDMS layer needs to be replaced after each run of experiment because its surface properties decays when it is exposed to the fluorescent solution. The thickness of each PDMS layer is estimated through a calibration process. This thickness is important for the calculation of the actual time of breakthrough from the total time of breakthrough.  31  2.1. Setup  2.1.3  Procedure  As explained earlier, the injection module delivers the fluorescence solution to the GDL sample by means of the injection channel. The nominal flow rate for the injection is 1.11 × 10−11 m3 /s. This flow rate is calculated based on the liquid water produced on the catalyst layer of a PEM fuel cell, when the cell is working at the flooding condition. It is worth mentioning that none of the previous studies ([LSD06], [BSLD07], [BHDS08], and [BSD08]) used the exact flooding flow rate which determines the capillary number (Ca) of the flow, and consequently the regime of flow in porous media. For an electrolysis reaction, Faraday’s law is used to correlate the mass of species produced in the reaction to the electric charges transferred. This law states ([LDF03]): QM m= (2.1) F z − During electrolysis, the mass of byproduct (m) formed on the electrode is directly proportional to: 1) the amount of electricity conducted by means of the same electrode (Q), and 2) the molar mass of the byproduct (M ). − During electrolysis, the mass of byproduct (m) formed on the electrode is inversely proportional to the number of electrons transferred by means of each ion (z) multiplied by Faraday constant (F = 96485C/mol). The time derivative of Equation 2.1 gives dm iM = dt zF  (2.2)  Dividing Equation 2.2 by the density of the substance (ρ) will result in dV iM = dt zF ρ  (2.3)  where i presents the current. For the electrochemical reaction occurring in the cathode catalyst layer of a PEM fuel cell, two ions (H + ) are produced per water molecule, i.e. z = 2. Substituting the values of M = 18.01528g/mol and ρwater = 1000kg/m3 in the equation 2.3, the volumetric rate of water production is obtained as dV dt  m3 sec  = 9.335 × 10−11  m3 col  i(A)  (2.4)  32  2.1. Setup Substituting the limiting current density of i = 1.4A/cm2 (representing a flooding condition) and the reactive area corresponding to a disk diameter of d = 3.3mm (considered in the experiment) results in the water flow rate value of Q = dV /dt = 1.11 × 10−11 m3 /s. In order to investigate the effect of flow rate on the time of breakthrough and the breakthrough pressure, all the experiments are conducted with two flow rates. The nominal flow rate corresponds to the flooding condition, and an arbitrary flow rate which was chosen to be ten times more than the nominal flow rate (1.11 × 10−10 m3 /s). Previous studies for determination of GDL saturation used a higher flow rate to simulate the condition of flooding [SH11]. In this study, both the flooding flow rate and the high flow rate are used to investigate the effect of increased flow rate on the breakthrough time. A GDL sample with the diameter of 25.4mm is placed into the injection channel. A cap is placed on the GDL applying enough pressure on the sample to prevent any leakage. Then, the system is covered by a PTFE cylinder providing a dark environment (see Figure 2.4). Four Toray R (Toray Industries, Inc., Chuo-ku, Tokyo, Japan) GDL samples are used in this study (TGP-H-30 (110µm), TGP-H-60 (190µm), TGP-H-90 (280µm) and TGP-H-120 (370µm)). For each thickness, a sample with no hydrophobic treatment (hydrophilic sample) and with a 40% PTFE treatment (hydrophobic sample) are considered to study the effect of hydrophobicity on the measurements. Each experiment was repeated three times to ensure the reproducibility of the results. Calibration Two aspects of the calibration process is addressed in this section: 1) calibration due to the gap between the last micro-needle, and 2) calibration due to the expansion of the setup components. As explained earlier, the last micro-needle sends the signal indicating the presence of liquid water at the bottom surface of the GDL. Although the needle is installed close to the bottom surface of the GDL, but there is a distance between the needle and the GDL surface. In addition, the thin PDMS layer covering the surface of the injection channel provides an additional gap between the needle and the GDL. Thus, the time required for liquid water to pass across the gap between the needle and the GDL sample should be considered in estimating the breakthrough time. In order to measure the calibration time (tcalibration,gap ) for each new PDMS layer treated on the injection channel, a voltage sensor was installed at the exit port of the channel. The time 33  2.2. Results between the signal sent by this voltage sensor and the signal of the last needle is considered as the calibration time. Since the calibration time is a function of the flow rate, it is measured for each flow rate. The second portion of calibration time is the time required for the setup components to expand prior to the injection process starts (tcalibration,expansion ). The flow rate and the time required for breakthrough are two quantities used in order to estimate the water content in GDLs. The experimental setup components were perfectly sealed to avoid any leakage. However, expansion (even small) of the tubing can drastically affect the results, as the working flow rate corresponding to the limiting current density is very low. Thus, the expansion problem should systematically be addressed. In this work, a technique is developed to measure the expansion of the experimental components at particular capillary pressures (Vexp = Vexp (Psys )). Knowing the nominal flow rate (Qnom ) and the time of breakthrough (tBT ), the calibrated volume of the injected water is estimated. Figure 2.5 shows the overall expansion of the components in the setup as a function of pressure for low and high flow rates. The measured breakthrough time is the total time elapsed between the signal sent by the last needle and the time the droplet is observed on the top surface of the GDL. The actual breakthrough time is obtained by deducting the two calibration times from the measured breakthrough time: tactual = tmeasured − tcalibration,gap − tcalibration,expansion  2.2  (2.5)  Results  The images captured during the experimental procedure are presented in this section. Figures 2.6 to 2.13 present one replicate of the experiments 3 3 conducted at two different flow rates (1.11 × 10−11 ms and 1.11 × 10−10 ms ) for the GDLs with two different thicknesses (110µm and 190µm) and two different PTFE loadings (0% and 40%). For each sample, four consecutive time steps are presented namely: t∗ = 0.25, t∗ = 0.50, t∗ = 0.75 and t∗ = 1.004 . The intensity–based images are then converted to contours which accurately show the location of highly saturated areas (red color mapping) and low saturated areas (blue color mapping). The variation of normalized pressure 5 with time is also presented for The time corresponding to each image (t∗ ) is normalized using the breakthrough time (tBT ) i.e. t∗ = t/tBT . 5 The pressure variations are normalized using the maximum value of the breakthrough pressure of all three replicates obtained for each thickness of GDL. 4  34  2.2. Results each case. The pressure variations in all the figures presented here show the same trend. Pressure values start increasing from the start of the injection process (t∗ = 0), till the breakthrough. The maximum value of capillary pressure is obtained at the breakthrough, and the pressure suddenly decreases as water forms a passage from the injection channel to the atmosphere. The measurements presented for both thicknesses show that average value of normalized pressures are higher for treated samples (40%) compared to untreated samples (0%) (compare Figures 2.6(e) and 2.7(e); and also Figures 2.10(e) and 2.11(e)). Pressure slightly increased through the 3 process of injection for the case of low flow rate (1.11 × 10−11 ms ) and as the breakthrough approaches, the change in pressure is visible. While for the 3 case of high flow rate (1.11 × 10−10 ms ) the pressure variation is considerable from the onset of injection. This can be observed by comparing the Figures 2.6(e) and 2.8(e), and also 2.10(e) and 2.12(e). The effect of hydrophobic content on the flow patterns are extensively analyzed and discussed in the next section.  35  2.2. Results  (a) The schematic of the experimental setup: the injection module (syringe pump, tubing, and the injection channel), the registration module (voltage sensors, pressure transducers, sensor terminal and data acquisition system), and the imaging module (vertical illuminator)  (b) Actual image of experimental setup developed at the UBC’s Advanced Thermofluidics (ATFL) Lab  Figure 2.1: Experimental setup. 36  2.2. Results  (a) The schematic of fluorescence microscopy technique  (b) Vertical illuminator used in this study  Figure 2.2: Details of fluorescence microscopy.  37  2.2. Results  (a) The schematic of injection channel and sensors  (b) Injection channel  Figure 2.3: Details of injection channel.  (a) GDL sample is placed into the channel.  (b) Sample holder provides pressure on the GDL sample.  (c) The cap closed the injection module providing dark environment.  Figure 2.4: The placement of the sample. 38  2.2. Results  Figure 2.5: Change in the volume of the setup versus capillary pressure: low flow rate (solid line), high flow rate (dashed line).  39  40  (c) t∗ = 0.75  (e) Variation of normalized pressure versus time  (b) t∗ = 0.50  (d) t∗ = 1.00  Figure 2.6: Flow configuration and capillary pressure variation versus time (sample thickness=110µm, 3 loading=0%, flow rate=1.11 × 10−11 ms )  (a) t∗ = 0.25  2.2. Results  41  (c) t∗ = 0.75  (e) Variation of normalized pressure versus time  (b) t∗ = 0.50  (d) t∗ = 1.00  Figure 2.7: Flow configuration and capillary pressure variation versus time (sample thickness=110µm, 3 loading=40%, flow rate=1.11 × 10−11 ms )  (a) t∗ = 0.25  2.2. Results  42  (c) t∗ = 0.75  (e) Variation of normalized pressure versus time  (b) t∗ = 0.50  (d) t∗ = 1.00  Figure 2.8: Flow configuration and capillary pressure variation versus time (sample thickness=110µm, 3 loading=0%, flow rate=1.11 × 10−10 ms )  (a) t∗ = 0.25  2.2. Results  43  (c) t∗ = 0.75  (e) Variation of normalized pressure versus time  (b) t∗ = 0.50  (d) t∗ = 1.00  Figure 2.9: Flow configuration and capillary pressure variation versus time (sample thickness=110µm, 3 loading=40%, flow rate=1.11 × 10−10 ms )  (a) t∗ = 0.25  2.2. Results  44  (c) t∗ = 0.75  (e) Variation of normalized pressure versus time  (b) t∗ = 0.50  (d) t∗ = 1.00  Figure 2.10: Flow configuration and capillary pressure variation versus time (sample thickness=190µm, 3 loading=0%, flow rate=1.11 × 10−11 ms )  (a) t∗ = 0.25  2.2. Results  45  (c) t∗ = 0.75  (e) Variation of normalized pressure versus time  (b) t∗ = 0.50  (d) t∗ = 1.00  Figure 2.11: Flow configuration and capillary pressure variation versus time (sample thickness=190µm, 3 loading=40%, flow rate=1.11 × 10−11 ms )  (a) t∗ = 0.25  2.2. Results  46  (c) t∗ = 0.75  (e) Variation of normalized pressure versus time  (b) t∗ = 0.50  (d) t∗ = 1.00  Figure 2.12: Flow configuration and capillary pressure variation versus time (sample thickness=190µm, 3 loading=0%, flow rate=1.11 × 10−10 ms )  (a) t∗ = 0.25  2.2. Results  47  (c) t∗ = 0.75  (e) Variation of normalized pressure versus time  (b) t∗ = 0.50  (d) t∗ = 1.00  Figure 2.13: Flow configuration and capillary pressure variation versus time (sample thickness=190µm, 3 loading=40%, flow rate=1.11 × 10−10 ms )  (a) t∗ = 0.25  2.2. Results  2.3. Discussion  2.3 2.3.1  Discussion Flow characteristics  In this section, the patterns of water flow inside the GDL are analyzed from the images obtained. Although not all the active area and thickness of the GDL is visible in fluorescence microscopy due to the opaque nature of carbon fibers, the images provide basic insight into important phenomena occurring during water flow through the GDL. The application of fluorescence microscopy for studying water flow in thin porous media such as GDLs, is an emerging technique. Prior studies helped establishing the ground knowledge for this application, however the assumptions made in these studies evolved as the method matured ([LSD06], [BSLD07], [BHDS08], and [BSD08]). Qualitative analysis of images obtained from this technique helps to understand the hydrodynamics of flow which eventually leads to a more accurate model for liquid flow inside the porous medium. Based on this understanding, a numerical pore-network model is developed and presented in the next chapter. The critical concepts and criteria used in the model (e.g., pore saturation, invading criteria, hydrodynamic pressure term) are derived based on the observation and analysis of the flow configuration and the measurement of the flow characteristics discussed in this section. The presented images are for a fuel cell working at the flooding condi3 tion which is equivalent to the injection flow rate 1.11 × 10−11 ms . All the prior studies which employed fluorescence microscopy as the visualization technique used a high flow rate to inject fluorescent solution into the GDL ([LSD06], [BHDS08], and [BSD08]). However, in this study, the injection flow rate is calculated base on the rate of water production at the flooding condition. The sample with the smallest thickness (110µm) is considered as it provides more details on the flow patterns. For this thickness, images obtained for both treated (hydrophobic) and untreated (hydrophilic) samples are presented, compared and discussed. Invasion pattern Water invades the pores quite differently in a hydrophilic medium compared to a hydrophobic medium. In a hydrophilic medium, the majority of pores are available to be invaded by liquid water due to distribution of local contact angle for each pore. In a hydrophobic medium, on the other hand, the pores are less vulnerable for water invasion. In previous studies, only one GDL sample was used and the observations from fluorescence microscopy 48  2.3. Discussion focused on the evolution of liquid water in that particular sample without any attempt to study the effect of GDL’s hydrophobicity on the flow of water ([BSD08]). In this study, an untreated GDL sample without any PTFE loading is used as the base sample, along with a treated GDL sample with 40% PTFE loading. In this way, the effect of GDL’s hydrophobicity on the invasion pattern is systematically addressed and studied. For two-phase flow in porous media, the capillary pressure is defined as the pressure difference between the phases. In the case of the GDL, water is the invading (or displacing) phase and air is the invaded (or displaced) phase. The Young-Laplace equation correlates the capillary pressure to the hydrodynamic and geometrical properties of the porous medium. For a drainage process to occur inside the GDL, the liquid pressure should be high enough to overcome the pore capillary pressure such that 4γ cos θ (2.6) d where γ is the interfacial tension and θ is the contact angle. When a GDL is treated with a hydrophobic agent (PTFE), the local contact angle of pores changes. Thus, the liquid pressure should increase to overcome the resistance of the pore and make a droplet invade it. In addition, fewer numbers of pores are vulnerable to be invaded by water in a hydrophobic medium; whereas in a hydrophilic medium, the number of available pores for invasion is higher for a given liquid pressure. Thus, less surface area of the hydrophobic medium is expected to be invaded compared to the hydrophilic medium. Figure 2.14 compares the invaded area for the hydrophilic and hydrophobic samples. As it is encircled in the figure, several areas of the hydrophilic sample are invaded by water forming clusters (Figure 2.14(b)), while for the hydrophobic sample (Figure 2.14(a)), only two clusters are distinguished. Pliquid − Pgas >  Progression pattern The criteria for invading a particular pore are the availability of the pore (proximity to the water frontier) and the capillary pressure required for the pore to be invaded. After the initial invasion, liquid water flows through the thickness of the samples. Litster et al.[LSD06] briefly studied development of “additional paths” after the initial invasion of water. They observed only one dominant fluid path for a treated sample (with 10% PTFE loading). No comparison was made for the progression patterns between the treated and an untreated 49  2.3. Discussion  (a) Hydrophobic sample  (b) Hydrophilic sample  Figure 2.14: Comparison of initial invasion patterns: Less surface area of the hydrophobic sample is invaded (left) compared to the hydrophilic area (right) (the encircled areas compare the invaded clusters). GDL sample in their study. However, the pattern for further developments depends on the characteristics of the medium. For a hydrophobic medium, water prefers to move through a path that has already been developed. Thus, the excess water injected into the medium after the initial invasion continues to fill the pores which were already invaded. The higher contact angle of the dry pores makes them less vulnerable to future invasion. There would be no or very few new developments in the water frontier. Figure 2.15(a) and 2.15(b) compares two consecutive time steps for the hydrophobic sample. To facilitate the comparison between the images, the corresponding time for each image is normalized by the total time of breakthrough (t∗ = t/tBT ).The first image is chosen as the reference image at t∗ = 0.487 and the second image shows the new development of the invasion pattern in the sample at t∗ = 0.553. Only three new frontiers are visible in this figure (encircled). For the hydrophilic medium, on the other hand, there is no preferred pattern of invasion. The comparison between Figures 2.15(c) and 2.15(d) (captured at t∗ = 0.487 and t∗ = 0.554, respectively) reveals that the new development occurs for the hydrophilic sample and a relatively large area is invaded as flow is developing in the hydrophilic medium compared to the hydrophobic sample. As stated above, the only prior work which briefly discussed the effect of hydrophobicity on progression pattern is Litster et  50  2.3. Discussion al.[LSD06]. However, they only used one GDL sample with 10% PTFE loading, and they did not conduct any experiments on an untreated sample to enable the comparison for water flow between treated and untreated GDL samples.  Pore filling pattern For a single empty pore to be invaded by liquid water in a porous medium, a threshold pressure is required. The threshold pressure is determined with respect to the capillary properties of the medium. If the medium is hydrophilic the threshold pressure is negative (cos θ < 0) which means that the medium intakes the water without resistance and the system reaches to equilibrium once the medium is saturated. For a hydrophobic medium, the threshold pressure is positive. Thus, the liquid pressure increases, and once it reaches to the threshold pressure it invades the pore. As the saturation of the pore increases the chance for water to invade the adjacent pores increases accordingly. The pressure of the liquid phase also contributes to the invasion process. Thus, one can quantify the chance of invasion by introducing a hydrodynamic pressure as a pressure correction term given by Phydrodynamic = Spore × Pref (2.7) where Spore is the saturation of an individual pore and Pref is the reference pressure which can be replaced by the reservoir pressure for the numerical simulation (explained in the next chapter). For the hydrophobic and hydrophilic samples, the difference between pore saturation and pore–filling processes is illustrated in Figure 2.16. The images correspond to t∗ = 0.992, i.e., an instance right before the breakthrough. For the hydrophobic sample, at least four zones are determined with a high concentration of water (red color-mapping) indicating saturated pores. For the hydrophilic sample, on the other hand, there is only one highly saturated region though the numbers of invaded pores are higher (as discussed previously). The studies which used fluorescence microscopy technique to investigate water flow through GDLs correlated the image intensity obtained from the microscopy to the height of liquid column in the sample ([LSD06], [BSLD07], [BHDS08], and [BSD08]). In this study, however, the image intensity is correlated to the pore saturation. This idea is supported by experimental observations and the limitation of the method in terms of in-plan and throughplan resolutions. In summary, the main difference between the pore–filling process in the hydrophobic and hydrophilic media is that in the hydrophobic sample water continuously fills a single pore until the pore is fully saturated 51  2.3. Discussion or the phase pressure of water at other locations of the frontier is high enough for another invasion process. In the hydrophilic sample, multiple invasions at different locations of the frontier occur due to a low or negative threshold pressure.  2.3.2  Pressure and time at breakthrough  In this section, the pressure measurements and time required for breakthrough are presented and discussed. GDL samples with four different thicknesses are used. For each thickness, a sample with no hydrophobic treatment (hydrophilic sample) and with 40% PTFE treatment (hydrophobic sample) are considered to study the effect of hydrophobicity on the pressure and time of breakthrough. Experiments with each sample was conducted with 3 3 low and high flow rates (1.11 × 10−11 ms and 1.11 × 10−10 ms ). Figures 2.17 and 2.18 show the results and the corresponding errors for different samples. For the low flow rate and high PTFE loading samples, the capillary pressure is larger and the resistance to the flow increases as the thickness of the sample increases (Figure 2.17). This is in agreement with experimental results reported in previous studies ([BNB+ 05] and [FCSPS07]). However, the increase in the breakthrough pressure is not proportional to the thickness of the sample. For instance, for the two thin GDLs (110µm and 190µm), the difference in the breakthrough pressure is significant, i.e., a 22% increase from the 110µm sample to the 190µm sample; while the change in the pressure for the two thick samples (280µm and 370µm) is marginal (2%). The same trend is observed for the sample with no PTFE loading, i.e. the increase in the breakthrough pressure is 29% for two thin samples, and as the thickness increases the pressure increases only by 3%. Figure 2.17 also presents the corresponding times of breakthrough for the low flow rate. For the samples with no PTFE loading, less time is required for water to penetrate and reach the gas channel. The effect of the thickness on the time of breakthrough is insignificant for thicker samples, i.e. the time of breakthrough for the two thick samples (280µm and 370µm) differs marginally (only by 4.6%) compared to the difference between the time of breakthrough of the two thin samples (i.e. 13.9%). The breakthrough pressures for the high flow rate are shown in Figure 2.18. Compared to the low flow rate, the pressure values are generally larger in this case (similar to the findings reported in [BNB+ 05]). For the samples with no PTFE loadings, as the thickness increases the breakthrough pressure increases similar to the trend observed for the low flow rate cases. 52  2.3. Discussion For the samples treated with PTFE, the breakthrough pressure increases by the thickness. However, the breakthrough pressure for the thickest sample (370µm) is slightly lower (350P a less) than that obtained for the second thickest sample (280µm). The same trend is observed in the variation of the time of breakthrough, i.e. the maximum time of breakthrough has been obtained for the 280µm samples with and without PTFE loadings. For both cases of loadings, the time of breakthrough unexpectedly decreases for the thickest sample (370µm). The behavior of the breakthrough pressure and time changes as the flow rate increases by a factor of ten. For this case, both parameters do not increase monotonically as the thickness of the sample increases. For the hydrophilic sample, the breakthrough pressure and the time of breakthrough increase up to the thickness of 280µm, where the pressure is maximized and the time marginally increases afterwards. For the hydrophobic samples, both variables are larger compared to the hydrophilic samples since the resistance of the medium to the flow increases due to the presence of the hydrophobic agent. This has also been reported in previous studies ([BNB+ 05],[GIFP08] and [FCSPS07]). Similar to the hydrophilic samples, the maximum pressure was obtained for the 280µm sample. In this case, the value of the time of breakthrough is also maximum for the second thickest sample (280µm). For the thickest sample (370µm), the effect of different PTFE contents on the time and pressure at breakthrough were also studied. Figure 2.19 shows the variations of the time and pressure at the low flow rate. Both variables increase as the PTFE loading of the sample increases. The maximum values are obtained at the loading of 40wt%. For the high flow rate case, for which the pressures are larger than those obtained during the low flow injection rate, the pressure and time of breakthrough increase monotonically as the PTFE loading of the sample increases (see Figure 2.20). For both flow rates, the variation of the pressure follows the same trend as the time of breakthrough.  53  2.3. Discussion  (a) Hydrophobic sample at t∗ = 0.486  (b) Hydrophobic sample at t∗ = 0.553  (c) Hydrophilic sample at t∗ = 0.487  (d) Hydrophilic sample at t∗ = 0.554  Figure 2.15: Progression pattern for hydrophobic and hydrophilic samples (the encircled areas in Figures (b) and (d) show the new developed water frontier; note that for the hydrophobic sample (b) only three more clusters are invaded as water progresses, while for the hydrophilic sample (d) the number of invaded clusters is twice.  54  2.3. Discussion  (a) Hydrophobic sample at t∗ = 0.992  (b) Hydrophilic sample at t∗ = 0.992  Figure 2.16: Pore filling process for hydrophobic and hydrophilic GDLs; the encircled areas show the high saturated regions. For the hydrophobic sample four clusters are indicated as high saturated pores, while for the hydrophilic sample only one relatively-high saturated cluster is visible.  Figure 2.17: Measured pressure and time of breakthrough versus the thick3 ness of the sample (flow rate of 1.11×10−11 ms ) for different PTFE loadings. 55  2.3. Discussion  Figure 2.18: Measured pressure and time of breakthrough versus the thick3 ness of the sample (flow rate of 1.11 × 10−10 ms ).  56  2.3. Discussion  Figure 2.19: Variations of the pressure and time over PTFE loadings for the 3 thickest sample (370µm) (flow rate of 1.11 × 10−11 ms ).  57  2.3. Discussion  Figure 2.20: Variations of the pressure and time over PTFE loadings for the 3 thickest sample (370µm) (flow rate of 1.11 × 10−10 ms ).  58  Chapter 3  Numerical model A pore–network modeling approach is adopted in this research. This approach provides details of the transport phenomena occurring at the microscopic scale. The results of the model developed in this study can be interpreted and used to understand the phenomena involved in porous media. Moreover, the model can easily be modified to capture the existence of another phase in the porous media and the interaction between different phases. As discussed in Chapter 1, the pore–network approach has been used by researchers to model flow of liquid water through the porous layers of PEM fuel cells ([NK03], [GIFP07], [SW07], [CPQ+ 08], [MD11], [KPP11], [EHPP11] and [PPSK11]). Although all of these studies used pore–network modeling concepts, details of the algorithm employed are different depending on how the porous medium is represented, which transport mechanisms are assumed, and how the capillary pressures are treated within the algorithm. For instance, Nam and Kaviany [NK03] produced the pore–network model by stacking fibrous structures (fibrous structures represent the solid matrix of the porous medium (GDL)) and shifting the layers of fibrous structure on each layer. For the displacing algorithm, they used a much simpler discretization of species conservation equation [BSL06] (similar to conservation of mass in fluid dynamics). Gostick et al. [GIFP07] used a simpler structure for the pore–network representation by assuming a regular lattice of cubes interconnected with ducts of square cross-section as the throats. However, they simulated convection (using mass and momentum conservation equations) and diffusion transport processes (using Fick’s law for binary diffusion). Their algorithm used an interface-tracking subroutine, by scanning all of the throats in the network, to determine and label the potential open throats for the invasion process. This process of searching for the frontier can be computationally expensive if the size of simulated GDL sample is large. Sinha and Wang [SW07] constructed an irregular network of pores based on X–ray micro-tomography, and then used simple Poiseuille’s law for laminar flow to displace the liquid phase between the pores. Similar to the 59  Chapter 3. Numerical model pore–network model adopted in this research, Sinha and Wang recognized the throats for the drop in the pressure of liquid, and they assumed that all the volume of liquid in the network is accumulated within the pores and not the throats. Using a reverse approach, Chapuis et al. [CPQ+ 08] first constructed a 2D model for their numerical algorithm. The pore–network model in their study was formed using a randomly distributed non-overlapping disks which represents the fibrous medium. Then, the same specifications of the network were exported to a digitally controlled micro–milling machine to produce the artificial experimental porous medium. For the displacing algorithm, they used a “bond invasion percolation” algorithm, in which a threshold capillary pressure was pre-estimated for each “bond” (the “bond” was defined as the minimum distance between each two adjacent disks), and then the algorithm identified the lowest invasion threshold. Consequently, the invading phase (liquid) filled the “bond” and the procedure was repeated until breakthrough was achieved. Two additional constraints governed the invasion process in their study: 1) only the accessible bonds can be invaded by liquid, and 2) the “trapped” bonds cannot be invaded by the invading phase. Having the constraints to track the frontier of the invading cluster, and also the requirement to track the “trapped” cluster drastically increased the cost of computation in this method even though the model was presented in 2D. In addition, despite the fact that the invasion percolation algorithm proposed by Chapuis et al. successfully predicted the experimental observations, the applicability of the model to actual fluid flow in the Gas Diffusion Layer of the PEM fuel cell is disputed. In this research, a regular network of pores represents the porous structure of the GDL. Similar to the pore–network model proposed by Sinha and Wang [SW07], the model presented in this research uses a 3D network representation of the GDL. However, the network is regular and the its main characteristics such as pore diameter, throat diameter and the size of the network are determined with respect to the porosimetry data of the GDL sample used in experimental chapter (Chapter 2). The model assumes that the accumulation of water occurs in the pores and the throats do not withhold any liquid water [SW07]. However, any drop in fluid pressure is applied through the throats using Darcy’s law as the constitutive correlation between capillary pressure and the throat’s flow rate. The main driving mechanism is fluid’s pressure in the displacing algorithm, and a novel technique (a “pressure correction” term) is adopted by adjusting the pore pressure to employ the effect of pore saturation in the displacing mechanism. In this technique, the liquid saturation in a particular pore is considered while the fluid is in60  3.1. Numerical Scheme vading to that pore, i.e. if the invaded pore is a highly saturated pore then the “pressure correction” term resists the invasion process by increasing the total pressure of the pore, and if the invaded pore is a low–saturated pore, the “pressure correction” term facilitates the invasion process. The same mechanism exists for the invading pore: the ability of the pore for invasion increases as the pore saturation increases. By correlating this term to the pore–saturation, the algorithm will be capable of handling the cluster interface without the need for any additional interface–tracking subroutines. The chapter is organized as follows: Numerical Scheme provides a full description of the pore-network model. The structure and size of the network, the coordination number, and specifications and physical properties of the pores and throats are presented in this section. The generation of the network, and its treatment presenting the effect of the hydrophobic agent are described. The properties of the fluid and how they are implemented in the network are presented. The details of the displacing algorithm, the driving forces used to move the fluid in the network, the role of capillary pressure and pore saturation in fluid displacement, and the implementation of the initial and boundary conditions are also presented in this section. The next section, Verification, is a set of case studies applying the proposed algorithm to different networks with meaningful physical properties. Thus, an identical network is generated, and treated by the hydrophobic agent. Simulations are presented for samples with five hydrophobic fractions to qualitatively study the effect of the PTFE loading on the general configuration of fluid flow. In the Results section, the numerical scheme is modified and applied to the real experimental conditions which were discussed in the previous chapter (Chapter 2). The modifications made in the pore-network model in order to match the physical properties of the GDL samples are throughly discussed. The results of the simulations are then presented, discussed and compared against the experimental results obtained in previous chapter.  3.1  Numerical Scheme  The numerical method used in this study is “pore-network modeling” and it was originally developed to study oil and gas flow in reservoirs in petroleum industry [Blu01]. Figure 3.1 schematically describes the procedures involved in pore-network modeling. The idea is to represent the porous medium by a network of pores interconnected with throats. The network is modified to implement any hydrophobic properties in the medium and then 61  3.1. Numerical Scheme  Figure 3.1: Consecutive steps of pore-network modeling approach (Representation of the porous matrix: Courtesy of Institute of Energy Technology, Department of Mechanical and Process Engineering, Swiss Federal Institute of Technology Zurich ). flow equations are used to displace fluids in the medium. A typical porous medium has two sections: solid matrix which is the material of the medium, in the case of GDL, the carbon fiber is solid matrix. The second section is the void space in between the solid matrix which is called pore body. The GDL is represented by a regular network of pores. The x − y planes are parallel to the interface of the GDL and the catalyst layer, and z−direction is along the GDL thickness. Liquid water is injected into the GDL from the z = 0 interface. After the network is generated it is treated to become partially hydrophobic. This is done by assigning a high contact angle to a number of the pores. Then a displacing algorithm is used to move fluid between the pores and to simulate flow of water inside the porous medium.  3.1.1  Network generation  The porous medium is represented by a regular network of cells interconnected via throats. In a regular 3D network of pores, each cell is connected 62  3.1. Numerical Scheme to six adjacent cells (this number is called the coordination number). Figure 3.2 shows a typical cell with its adjacent cells and the corresponding index notations namely rear ([i-1,j,k]) and front ([i+1,j,k]), left ([i,j-1,k]) and right ([i,j+1,k]), and top ([i,j,k+1]) and bottom ([i,j,k-1]). The entire network is presented in Figure 3.3. The x − y planes correspond to the interface of the GDL-catalyst layer and the planes parallel to that. The z-direction is along the GDL thickness. The liquid water is injected through the GDL from the z = 0 plane.  Figure 3.2: Configuration of the unit cell and the adjacent cells. The cell is a representation of the physical void volume in the porous medium (pore body). The dimensions of the pores are determined using the porosimetry data reported in [MRFL03]. In this study, the pore diameter is randomly assigned in the range of 10µm < dpore < 20µm with a uniform distribution. The pore body withholds the fluid inside the porous medium, and the interface of the adjacent pores is the place where the pressure drop occurs. Similarly, in the simulation, the accumulation of fluid in the porous medium occurs inside the cells and all the pressure variations occur in the throats. Thus, the dimensions of the throats are determined based on the pressure drop in the medium and the thickness of the medium. The length of throat is lthroat = 11µm, and its diameter is randomly distributed in the range of 5µm < dthroat < 10µm. The number of cells in the z–direction (KM AX) is determined based on the thickness of the GDL sample. For the samples 63  3.1. Numerical Scheme  Figure 3.3: A sample of pore-network model with index notations (15 × 15 × 10). with the thickness of 110µm, KM AX = 10. For a disk of diameter of 3.3mm, the injection flow rate of 1.11 × 10−11 m3 /s is used to mimic the mass flow rate at the flooding condition. If the same area is used for simulation, a 150×150 network should be adopted, which is computationally expensive. Thus, the exposed area is reduced by a factor of 100, which results in a 15 × 15 network. To ensure the validity of the simulation, the total flow rates of liquid water injected into each sample were analyzed against the exposed area. Figure 3.4 shows the initial flow rates at the onset of simulations versus the exposed area for different networks. The figure shows the injection flow rate is 100 times less than the experimental condition as the exposed area in the numerical simulation is reduced by the same factor.  3.1.2  Network treatment  In order to implement the effect of hydrophobicity of the GDL sample in the simulation, contact angle values are assigned to the pores in the network. First, the network is considered as a hydrophilic medium as the 64  3.1. Numerical Scheme  Figure 3.4: Initial flow rate versus total exposed area for the numerical simulations; the distribution suggests that the initial flow rate is reduced by a factor of 100, as the network exposed area is reduced with the same factor. contact angle of the solid matrix (carbon fibers) is less than 90◦ ([WDN04]). Therefore, a random distribution of contact angle values in the range of 60◦ < θ < 90◦ is assigned to all of the pores to present the interfacial properties of the solid matrix. Then, based on the hydrophobic fraction of the GDL sample, the number of pores required to modified is calculated (hydrophobic pores). Subsequently, the contact angle of these pores are modified randomly between 90◦ < θ < 120◦ . To determine the number of hydrophobic pores, it is assumed that the mass of GDL sample is the summation of the fiber and the PTFE masses (msample = mf iber + mP T F E ). Thus, the PTFE weight fraction of a treated sample with respect to an untreated sample is calculated as %wt = ω =  msample − muntreated × 100 muntreated  (3.1)  Therefore, the mass of PTFE loading is mP T F E = ωmf iber . The volume of the GDL fiber is calculated based on the total volume of the sample and  65  3.1. Numerical Scheme sample porosity, , Vf iber = Vtotal − Vpore = Vtotal (1 − )  (3.2)  The total volume of the sample can also be estimated using the throat length (lthroat ) and the network size 3 Vtotal = (IM AX − 1)(JM AX − 1)(KM AX − 1)lthroat  (3.3)  Thus, the mass of the fiber and PTFE can be estimated as 3 mf iber = ρf iber (IM AX − 1)(JM AX − 1)(KM AX − 1)(1 − )lthroat (3.4) 3 mP T F E = ωρf iber (IM AX −1)(JM AX −1)(KM AX −1)(1− )lthroat (3.5)  For every hydrophobic pore in the network, a very thin layer of PTFE is considered to cover the pore. If the number of hydrophobic pore is N and the thickness of the PTFE layer is δ, the mass of PTFE can be expressed as mP T F E = N · ρP T F E · Vlayer  (3.6)  in which Vlayer = πd2 δ and d presents the diameter of the pore. The two above equations (3.5 and 3.6), result in N=  3 ωρf iber (IM AX − 1)(JM AX − 1)(KM AX − 1)(1 − )lthroat ρP T F E · πd2 δ  (3.7)  For ρf iber = 400kg/m3 , ρP T F E = 2200kg/m3 , lthroat = 11µm, δ = 0.1µm and average porosity of = 0.78, the total number of hydrophobic pores will be: N ≈ 0.75ω(IM AX − 1)(JM AX − 1)(KM AX − 1)  (3.8)  Figure 3.5 shows the treated networks used in this study. The blue and red spheres represent the hydrophilic and hydrophobic pores respectively. The properties of the networks and the fluids are summarized in Table 3.1.  3.1.3  Algorithm  The pore-network modeling approach is used to study the transient flow of fluid through the porous medium. As discussed in the introductory paragraph of this chapter, the application of the pore–network modeling to water flow in the GDL is very recent, and a few researchers have adopted this approach to simulate water flow and flooding in the porous layer of PEM fuel 66  3.1. Numerical Scheme  Table 3.1: The geometrical and hydrodynamic properties of the networks and fluids Variable Value IMAX JMAX KMAX Coordination number Pore diameter (dpore ) Throat diameter (dthroat ) Throat length (lthroat )  15 15 10, 17, 25, 33 6 10-20 µm 5-10 µm 11 µm  Water viscosity(µ) Water density(ρ) Water interfacial tension(γ) Air viscosity Air density Contact angle range (hydrophilic) Contact angle range (hydrophobic) Density of carbon fiber (ρf iber ) Density of PTFE (ρP T F E ) Thickness of PTFE on the pore surface (δ) Average sample porosity ( )  1.002×10−3 N s/m2 1000 kg/m3 0.072 N/m 1.84×10−5 N s/m2 1.190 kg/m3 60-90 ◦ 90-120 ◦ 400 kg/m3 2200 kg/m3 0.1µm 0.78  67  3.1. Numerical Scheme  Figure 3.5: Treatment of networks with the fraction of hydrophobic agent: (a)f = 0.10, (b)f = 0.20, (c)f = 0.30 and (d)f = 0.40 cells ([NK03], [GIFP07], [SW07], [CPQ+ 08], [MD11], [KPP11], [EHPP11], [PPSK11]). In all of these studies, the porous medium is represented by a network of pores interconnected by throats. The network can be regular (e.g. [GIFP07], [MD11] and [EHPP11]) or irregular (e.g. [SW07] and [CPQ+ 08]). The geometrical properties of the network (i.e. network size and pore size distribution) are usually derived from the porosimetry data available for the GDL sample under study ([GIFP07]), while in some cases an artificial pore medium is developed based on the numerically generated pore structure ([CPQ+ 08]). Then a displacing algorithm is used to force the liquid water move through the network until a termination criteria (such as breakthrough) is met. In this research, the simulation is focused on the flow of water inside the GDL from the start of injecting the fluid into the medium up to the breakthrough (identical to the experimental conditions explained in chapter  68  3.1. Numerical Scheme 2). The data collected in the experimental section is used for modifying the model and setting the boundary condition. In particular, the measured pressure values are employed as the boundary condition. The first row of the network (plane z = 0) is the catalyst layer. It is assumed that all the surface of the catalyst layer is active at the condition of flooding. Thus, all the cells in this row are filled with pressurized liquid water (the value of the pressure is set according to the experimental measurement). During the simulation, these cells remain filled with water and provide water to the entire network. The saturation of the cells6 remain constant during the simulation (Scell = 1). Fluid inside a cell invades the adjacent cells if the following two criteria are met ([SW07]): 1) the fluid pressure in the invading cell is high enough to overcome the capillary pressure of the connecting throat and, 2) the pressure difference between the cells is high enough to overcome the resistance due to viscous friction in the throat. The first criterion considers and deploys the surface properties of medium. The significant of capillary-driven flow inside the porous medium is observed in this term: Pliquid − Pgas > Pcapillary  (3.9)  in which Pcapillary is determined with respect to the diameter of the tube (throat) and surface properties of the fluid using Young-Laplace equation: ∆P = γ  1 1 + R1 R2  (3.10)  where ∆P is the pressure difference across the fluid interface, γ is the interfacial tension of the fluid, and R1 and R2 are the principal radii of curvature. For a circular tube with the diameter of d, the Young-Laplace equation (3.10) can be simplified to ∆P =  4γ cos θ d  (3.11)  in which θ is the contact angle between the fluid and the medium. For the second criterion, the friction resistance of the throat is estimated using the Poiseuille equation for laminar flows ∆P = 6  128µLQ πd4  Cell saturation is the ratio of fluid volume to cell volume: Scell =  (3.12) Vf luid Vcell  69  3.1. Numerical Scheme where ∆P is the pressure gradient at two ends of the tube (throat) with length L, µ is the dynamic viscosity of the moving fluid, Q is the flow rate of the fluid and d is the diameter of the tube (throat). As a novel technique in the displacing algorithm of this study, the pressure of each cell includes two terms as presented in equation 3.13: 1) the hydrodynamic term which is a function of local saturation (Scell ) and the reference pressure (the pressure at the boundary), and 2) the capillary term which depends on the fluid contact angle and the diameters of the throat connected to the cell. Pcell = Phydrodynamic + Pcapillary  (3.13)  The direct observation of pore–filling pattern presented in the experimental chapter (2) motivated the introduction of the hydrodynamic term as a pressure correction term (see Chapter 2, “Discussion” section under “Pore filling pattern” subsection). As discussed in that section (“Pore filling pattern”), the probability of water invasion from a particular pore to an adjacent pore increases monotonically as the saturation of the source pore increases. Thus, the pressure of liquid phase contributes to the invasion process. This “probability” can be quantified by introducing the hydrodynamic pressure. The hydrodynamic term (Phydrodynamic ) is zero or has a positive value for the cases that liquid is injected into the sample (Pref > 0). If no fluid exists in the cell (Scell = 0), the hydrodynamic term would be zero. The capillary term (Pcapillary ) is a function of contact angle (θ). For the hydrophilic cells, where cos θ > 0, the capillary term is positive. This implies that the capillary is the driving force in hydrophilic pores, contributing to the total pressure of the cell. For the hydrophobic cell, on the other hand, the capillary term is negative (cos θ < 0). Similar to all the porenetwork models handling the capillarity (e.g. [GIFP07], [SW07], [CPQ+ 08], [MD11], [KPP11], [EHPP11]), the capillary term is only active at the fluidfluid interface. In the displacing algorithm adopted in this research, the pore saturation determines if the capillary pressure is active for a pore and should be considered. Thus, if the pore is full (S = 1), the capillary term is deactivated, otherwise the capillary term is added or subtracted from the hydrodynamic pressure depending on the value of the contact angle. Either term (Phydrodynamic or Pcapillary ) can be the dominant term depending on the condition of the flow and the regime under study. For a high capillary number7 (Ca >> 1), the viscous force is dominant over the surface tension. Similarly in the simulation, the capillary term is negligible and the share of 7  Ca =  µV γ  70  3.1. Numerical Scheme the hydrodynamic term in the total pressure of the cell is higher. For a low capillary number (Ca << 1), on the other hand, the surface tension force dominates and affects the flow configuration. Figure 3.6 summarizes the displacing algorithm used in this study. In each time step, the cell pressure (Pcell ) is estimated for each throat connected to the cell using equation 3.13. In fact, the cell has six pressures for every six direction as it is shown in Figure 3.2. Thus, every throat has two pressures at the two ends which correspond to the cells attached to it (see Figure 3.7). Using equation 3.12, the flow rate in the throat is calculated by  71  Figure 3.6: Flowchart of the displacing algorithm.  3.1. Numerical Scheme  72  3.1. Numerical Scheme  Figure 3.7: Connecting throat between two adjacent cell.  Qthroat =  πd4 (Pcell,1 − Pcell,2 ) 128µlthroat  (3.14)  The calculated flow rates in each throat are used to update the volume of the fluid in each cell. Thus, 5  Qthroat,i × ∆t  Vf luid,t = Vf luid,t−1 +  (3.15)  i=0  where Vf luid,t and Vf luid,t−1 , are the volume of the fluid at time t and t − 1 respectively. Then the saturation of each cell is estimated as Scell =  Vf luid Vcell  (3.16)  Using the cell saturation (Scell ), the hydrodynamic pressure of each cell is updated based on Phydrodynamic = Scell × Pref (3.17) Then, the total pressure of the cell is updated using Equation 3.13. The updated cell pressures are used to calculate the flow rate in each throat in next time step. The algorithm stops when the condition of breakthrough is achieved. Breakthrough is defined as when the pressure of at least one of the cells in the last row (z = KM AX) is equal to the breakthrough pressure. The value of the breakthrough pressure is obtained from the experiment (see chapter 2).  73  3.2. Results  3.2  Results  Numerical schemes are usually verified based on two aspects: 1) verification of the code, for which incorrect implementation of conceptual models, the error in inputs and other sources of errors are studied, and 2) verification of the calculation, which involves error estimation for a single calculation and grid convergence. Consistency checks for examining basic physical relationships expected in the solution, and grid refinement are also included in the verification process [AIA98]. In the current pore-network modeling, Network generation (section 3.1.1) and Network treatment (section 3.1.2) are two main random processes when the network is created. The effect of randomness on network generation is studied in Chapter(4) where the actual conditions of the experiment are reproduced by generating a fresh network geometry for each case under study. The effect of randomness on network treatment is studied in this section. Herein, a unique reference network is generated and treated with a specific hydrophobic fraction (f ). For each fraction, the reference network is treated three times to produce three replications. Thus, the randomness effect in network treatment is studied through analysis of the errors in this section. Finally, the model is used to determine the total saturation of invading fluid (water), as the hydrophobic fraction increases. The algorithm is applied to five hydrophobic fractions (f ) to verify the numerical scheme. For all the cases studied here, an identical geometry of network (reference network ) is generated (i.e. a network with the same pore size and throat size distributions). The reference network is then treated with five fractions of hydrophobic agent (f = 0.0, 0.1, 0.2, 0.3 and 0.4) to study the effect of hydrophobicity on the flow configuration. Three replicate simulations are performed for each case. Finally, the model is used to determine the total saturation of invading fluid (water), as the hydrophobic fraction increases. The criteria for verification are: 1) the total saturation of invading fluid (water)as a function of the hydrophobic fraction, and 2) the time of breakthrough as a function of the hydrophobic fraction. A network of 15 × 15 × 10, corresponding to the thinnest GDL sample (TGP-H-30, 110 µm), is developed. Figure 3.5 presents an example of treated networks. As explained earlier, for each hydrophobic fraction, three random treatments of the network are developed. Figure 3.8 shows an example of different treatments for f = 0.4. Following the above procedure, the saturation contours of flow through the GDL for five different loadings are obtained. Figures 3.9 to 3.33 presents these contours for different loadings. For each case, five sections in the x− 74  3.2. Results  Figure 3.8: Three replications of treated network (f = 0.4) and y− directions (x∗ , y ∗ = 0.05, 0.25, 0.50, 0.75 and 1.00) and four cross sections in the z− direction (z ∗ = 0.25, 0.50, 0.75 and 1.00) are presented. The simulation results are also presented for five time instances presenting three periods: the initial invasion (t∗ = 0.01), progression of flow during the injection (t∗ = 0.25, 0.50 and 0.75) and at breakthrough (t∗ = 1.00). These figures are then analyzed in the next section to extract the details of flow in the GDL, and to compare them with the patterns obtained from the experimental observations.  75  (b) y−planes  (c) z−planes  Figure 3.9: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.01.  (a) x−planes  3.2. Results  76  (b) y−planes  (c) z−planes  Figure 3.10: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.25.  (a) x−planes  3.2. Results  77  (b) y−planes  (c) z−planes  Figure 3.11: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.50.  (a) x−planes  3.2. Results  78  (b) y−planes  (c) z−planes  Figure 3.12: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.75.  (a) x−planes  3.2. Results  79  (b) y−planes  (c) z−planes  Figure 3.13: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 1.00.  (a) x−planes  3.2. Results  80  (b) y−planes  (c) z−planes  Figure 3.14: Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 0.01.  (a) x−planes  3.2. Results  81  (b) y−planes  (c) z−planes  Figure 3.15: Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 0.25.  (a) x−planes  3.2. Results  82  (b) y−planes  (c) z−planes  Figure 3.16: Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 0.50.  (a) x−planes  3.2. Results  83  (b) y−planes  (c) z−planes  Figure 3.17: Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 0.75.  (a) x−planes  3.2. Results  84  (b) y−planes  (c) z−planes  Figure 3.18: Saturation contours (sample thickness=110µm, PTFE loading=10%) at t∗ = 1.00.  (a) x−planes  3.2. Results  85  (b) y−planes  (c) z−planes  Figure 3.19: Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 0.01.  (a) x−planes  3.2. Results  86  (b) y−planes  (c) z−planes  Figure 3.20: Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 0.25.  (a) x−planes  3.2. Results  87  (b) y−planes  (c) z−planes  Figure 3.21: Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 0.50.  (a) x−planes  3.2. Results  88  (b) y−planes  (c) z−planes  Figure 3.22: Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 0.75.  (a) x−planes  3.2. Results  89  (b) y−planes  (c) z−planes  Figure 3.23: Saturation contours (sample thickness=110µm, PTFE loading=20%) at t∗ = 1.00.  (a) x−planes  3.2. Results  90  (b) y−planes  (c) z−planes  Figure 3.24: Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 0.01.  (a) x−planes  3.2. Results  91  (b) y−planes  (c) z−planes  Figure 3.25: Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 0.25.  (a) x−planes  3.2. Results  92  (b) y−planes  (c) z−planes  Figure 3.26: Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 0.50.  (a) x−planes  3.2. Results  93  (b) y−planes  (c) z−planes  Figure 3.27: Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 0.75.  (a) x−planes  3.2. Results  94  (b) y−planes  (c) z−planes  Figure 3.28: Saturation contours (sample thickness=110µm, PTFE loading=30%) at t∗ = 1.00.  (a) x−planes  3.2. Results  95  (b) y−planes  (c) z−planes  Figure 3.29: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.01.  (a) x−planes  3.2. Results  96  (b) y−planes  (c) z−planes  Figure 3.30: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.25.  (a) x−planes  3.2. Results  97  (b) y−planes  (c) z−planes  Figure 3.31: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.50.  (a) x−planes  3.2. Results  98  (b) y−planes  (c) z−planes  Figure 3.32: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.75.  (a) x−planes  3.2. Results  99  (b) y−planes  (c) z−planes  Figure 3.33: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 1.00.  (a) x−planes  3.2. Results  100  3.3. Discussion  3.3 3.3.1  Discussion Comparison between the flow patterns obtained from numerical results and experimental observations  As explained in the Flow characteristics section of Chapter 2, three patterns are observed for the water flow through the GDL: initial invasion, progression and pore-filling patterns. The patterns differ as the hydrophobic fraction of the porous sample increases. In this section, water distribution obtained from the pore-network modeling for for an untreated (f = 0.0) and treated samples (f = 0.4) are qualitatively compared to the experimental observations for each pattern. Initial invasion pattern As observed from the experiment, fewer pores are vulnerable to be invaded by water in a hydrophobic medium compared to a hydrophilic medium where the number of available pores for invasion is higher. Thus, less surface area of the hydrophobic medium is expected to be invaded compared to the hydrophilic medium. Figure 3.34 shows the effect of the hydrophobic fraction (f ) on the invaded area during the initial stages of injection (t∗ = 0.25). The cross section is chosen to be close to the bottom layer of the sample where the chemical reaction occurs and water is produced (z ∗ = 0.25). For the sample with no hydrophobic agent (f = 0.0) all the pores are filled (totally or partially) with water (Figure 3.34(a)). As the hydrophobicity increases, the fraction of uninvaded pores increases. For the sample with f = 0.1, two uninvaded areas are distinguished (Figure 3.34(b)). The same trend is observed as the hydrophobic fraction increases where, for the highest hydrophobic fraction (f = 0.4), a comparable area of the sample remains untouched (as encircled in Figure 3.34(e)). This is in good agreement with the initial invasion pattern discussed in Chapter 2.  101  (e)  (d)  (c)  Figure 3.34: Numerical results demonstrating the initial invasion pattern at z ∗ = 0.25: (a) f = 0.0 (b) f = 0.10 (c) f = 0.20 (d) f = 0.30 (e) f = 0.40.  (b)  (a)  3.3. Discussion  102  3.3. Discussion Progression pattern As explained earlier, after the initial invasion, liquid water flows through the thickness of the samples and the characteristics of the medium determines the pattern for further developments. For the hydrophobic medium, excess water injected into the medium fills the pores which are already invaded. In essence, the higher contact angle of the dry pores makes them less vulnerable to future invasion and a very few new branching occurs on the water frontier. Figure 3.35 shows the numerical results obtained for the hydrophobic sample at x∗ = 0.50,y ∗ = 0.50 and z ∗ = 0.50 for two consecutive time steps t∗ = 0.25 (i.e., reference time) and t∗ = 0.5. As the figure suggests, only a few new frontiers are visible (encircled) for the hydrophobic medium. For the hydrophilic medium, on the other hand, there is no preferred pattern of invasion and a large area is invaded as flow is developing (Figures 3.36(a) and 3.36(b)). Pore-filling pattern As observed experimentally, water continuously fills a single pore in the hydrophobic medium until the pore is fully saturated or the phase pressure of water at other locations of the frontier is high enough for another invasion process. In a hydrophilic medium, on the other hand, multiple invasions at different locations of the frontier occur due to a low or negative threshold pressure, and it is not necessary to have highly saturated pores in order for the invasion process to occur. Figures 3.37(a) and 3.37(b) show the numerical results for the untreated and treated samples, respectively. The red color mapping indicates the highly saturated areas as the color spectrum gets close to blue, the saturation is reduced. For the untreated sample (Figure 3.37(a)), almost all of the pores in the cross section (z ∗ = 0.25) are either completely or partially filled with water. Nevertheless, the portion of fully saturated clusters is few. These clusters are encircled in Figure 3.37(a). For the treated sample, on the other hand, fewer pores are invaded as is indicated by the green-blue areas in Figure 3.37(b). In this case, however, water manages to completely fill the invaded pores. Thus, the majority of the invaded pores are highly saturated as is indicated in the figure.  103  3.3. Discussion  3.3.2  Effect of hydrophobicity on the flow pattern  The pore network model presented is used to study the effect of hydrophobicity on the flow pattern and sample saturation. Figure 3.38 shows the saturation of each cell at the breakthrough. For each hydrophobic fraction, the flow is presented in all three directions. Also, to give a better image of the flow, three slices of the sample are provided in the through-plane directions (at x∗ , y ∗ = 0.05, 0.50 and 0.95). For the in-plane direction, a plane close to the top surface of the GDL (z ∗ = 0.75) is presented. The figure suggests that the portion of untouched spots inside the sample is increased as the hydrophobic fraction (f ) is increased. In other words, the area with zero or low saturation (blue-green color mappings) are larger for the samples with larger f . Also, water is uniformly distributed for the sample with no hydrophobic loading (f = 0.0) as it is shown in the first column of the figure, which presents the change in the pattern of water distribution in the planes of x∗ =constant. For the same sample, saturation is constant through the thickness of the GDL, except in the area close to the catalyst layer (z ∗ = 0), where water is produced and obviously local saturation is higher. At the area close to the gas channels (z ∗ = 1), the configuration is not uniform. However, the constant saturation through the thickness of the GDL, especially for the area far from the boundaries (0 < z ∗ < 1), suggests a uniform distribution of water. Also, for the sample with no treatment (f = 0.0), there is no preferred path for water to flow through and to find its way out. As the hydrophobic fraction increases, the pattern of water distribution changes significantly. Although the capillary number (Ca) and mobility ratio (M = µ2 /µ1 , which is the ratio of viscosity for the displaced and displacing fluids) are almost constant for all the cases studied here, the configuration of water distribution shifts toward viscous-fingering as the hydrophobic fraction increases. Comparison between the flow pattern simulated for the sample with f = 0.2 and the sample with no treatment (f = 0.0) verifies this fact. For instance, the pattern of water distribution in the mid-plane (x∗ = 0.50) of the sample with treatment (f = 0.2) exhibits four distinguishable peaks at the water frontier; while only two peaks are observed for the sample with no treatment (f = 0.0). As the hydrophobic fraction further increases, non-uniformity in the flow pattern develops. As a result of this non-uniformity, a larger area of the cross section close to the gas channel (z ∗ = 1) is untouched by liquid water as it is shown in the last column of the figure which presents the distribution of water in a plane parallel to the catalyst layer and close to the gas channel (z ∗ = 0.75)). For the sample with no treatment (f = 0.0), almost all the surface of this plane is covered 104  3.3. Discussion  Table 3.2: Normalized time of breakthrough obtained from the simulations (Averages of the replicates are provided in brackets). f=0.0 0.3755  f=0.1  0.3461 0.3869 (0.3695)  0.3868  0.3690 0.4110 (0.3889)  f=0.2 0.5081  f=0.3  0.5551 0.4869 (0.5167)  0.8018  0.7452 0.7311 (0.7594)  f=0.4 0.9714  0.8974 1.0000 (0.9563)  with water. The saturation of water is higher compared to the other cases with higher hydrophobic fractions. As the hydrophobic fraction increases, a few blue areas, which indicate the existence of uninvaded cells and pores, on this plane (z ∗ = 0.75) starts to develop. These areas provide passages for the reactant gases to pass through and reach the reaction sites on the catalyst layer. The model developed here is used to determine saturation in samples with different hydrophobicity. Figure 3.39 shows the saturation curves for difference hydrophobic treatments. For the sample with the highest hydrophobic fraction, the majority of the cross sectional area has very low local saturation. The overall saturation increases as the hydrophobic content decreases. The pattern presented in this figure is in agreement with those reported in [PLP08] and [PLP12]. Finally, the time of breakthrough is used to verify the numerical scheme. As explained in the experimental chapter (Chapter 2), the time of breakthrough is directly proportional to the hydrophobicity of the GDL sample, i.e. as the amount of PTFE loading increases, the time required for water to breakthrough the sample and form a sample-spanning water cluster increases. Table 3.2 summarizes the time of breakthrough obtained from the numerical scheme. Figure 3.40 presents the effect of the hydrophobic fraction on the time of breakthrough obtained numerically. As confirmed by the experiments, the time of breakthrough increases as the hydrophobic fraction increases.  105  3.3. Discussion  (a)  (b)  Figure 3.35: Progression pattern for the treated sample at: (a) t∗ = 0.25 ; (b) t∗ = 0.50. 106  3.3. Discussion  (a)  (b)  Figure 3.36: Progression pattern for the untreated sample at: (a) t∗ = 0.25 ; (b) t∗ = 0.50. 107  3.3. Discussion  (a)  (b)  Figure 3.37: Numerical results demonstrating the pore-filling pattern at z ∗ = 0.25 for: (a) treated sample; and (b) untreated sample.  108  3.3. Discussion  (a) x-planes  (b) y-planes  (c) z-planes  Figure 3.38: Effect of hydrophobic fraction (f ) on flow configuration at breakthrough. 109  3.3. Discussion  Figure 3.39: Effect of hydrophobic fraction (f ) on saturation at the breakthrough.  110  3.3. Discussion  Figure 3.40: Effect of hydrophobic fraction on the time of breakthrough.  111  Chapter 4  Model implementation In the experimental chapter (Chapter 2), a setup was designed and employed to study the flow of liquid water through GDL samples. The effect of sample properties, such as thickness and PTFE loading, on the time of breakthrough and breakthrough pressure were studied. Four thicknesses of GDL samples were used as the experimental matrix. For each thickness, a sample with no treatment of hydrophobic agent (hydrophilic sample 8 ) , and 40% PTFE loading (hydrophobic sample 8 ) were considered. The obtained results suggested that the flow of water inside the GDL is not a linear phenomenon in terms of time of breakthrough, i.e. when the thickness of the sample is doubled, the time of breakthrough is not necassarily doubled. The breakthrough pressure measured was also higher for a hydrophobic sample compared to that of a hydrophilic sample. The same was observed for the time of breakthrough. In essence, the comparison made between the time of breakthrough for the samples, suggests that for the hydrophobic sample a longer time is required for water to penetrate into the sample, travel through it and reach the other side of the sample. As explained in Chapter 2, this trend suggests that the amount of water inside the GDL (Vwater ) and hence saturation (S) at the time of breakthrough is higher for the hydrophobic sample compared to the hydrophilic sample if the flow rate into the GDLs (with different hydrophobicity) is assumed to be constant (Vwater = Q×tBT ). However, the data reported in the past [LNK10],[PP09],[GFI+ 06], [GIFP08] suggests the opposite (i.e., the saturation of hydrophilic samples is larger compared to that of the hydrophobic one). Thus, the assumption of constant flow rate for injecting water into the sample may overestimate the saturation. As a result, correct and precise measurement of the flow rate is required. Accurate measurements of the flow rate, however are very difficult (if not impossible) considering the small value of the flow rate (1.11 × 10−11 m3 /s) used in this study. The numerical scheme, on the other hand, can be used as a tool to 8  Since only two PTFE loadings are studied in this chapter , the term hydrophobic sample is used for the treated sample, and hydrophilic sample is used for the sample with no treatment.  112  Chapter 4. Model implementation replace the challenge associated with varying injection flow rates. If the numerical model is well-established and tested against available experimental data (e.g. flow patterns and breakthrough time), it can be used to evaluate the flow rate. Consequently, a reliable value for the injection flow rate can be obtained which leads to a correct estimation for the saturation. Chapter 3 presented the qualitative analysis of water flow through the GDL. Only one thickness of GDL sample was studied, with five PTFE loadings, while only two boundary pressures were set the same as the experimental measurements (f = 0.0 and f = 0.4). The boundary conditions for three other PTFE loadings were extrapolated using two measured values. In this chapter, the data obtained from the experiments are used to refine the numerical scheme. The physical properties of the GDL samples are used to properly set the network characteristics. All four GDL thicknesses used in experiments are simulated and the measured pressures are used as the boundary conditions. The numerical model is then validated against the experimental data (time of breakthrough) and the saturation values are obtained numerically. It is also shown that when the pressure is used as the boundary condition, the flow rate is not constant. For the hydrophobic sample, the flow rate is considerably lower than the sample with no treatment. Four networks with different thicknesses are used to model four samples under study. For each thickness, untreated and treated networks are generated to reproduce the experimental conditions for samples treated with 0% and 40% PTFE loadings, respectively. Three replicates for each numerical simulation are performed. For each replicate, a network is randomly generated and treated independently. Thus, the effect of network generation on the performance of the numerical scheme is also studied. The properties of the networks were summarized in Table 3.1. As mentioned earlier, the pressure is used as the boundary condition. The values of the pressure implemented are those obtained from the pressure transducer installed close to the bottom layer of the GDL sample (see Chapter 2 for the details of pressure measurements). Table 4.1 summarizes the normalized time of breakthrough and the pressure measurements used as the boundary conditions for the simulation. Table 4.2 summarizes the size and boundary condition used for each simulation.  113  4.1. Results  Table 4.1: Normalized time of breakthrough and breakthrough pressure (experiment). Sample thickness (µm) 110 190 280 370  Time (Pressure) hydrophilic hydrophobic 0.5138 0.5856 0.7316 0.7657  (5318.6P a) (6495.2P a) (8502.5P a) (8687.5P a)  0.6052 (7041.7P a) 0.7529 (9124.3P a) 0.9253 (9834.5P a) 1.0000 (10121.0P a)  Table 4.2: Network properties and boundary condition.  4.1  Sample thickness (µm)  Network size  110 190 280 370  15 × 15 × 10 15 × 15 × 17 15 × 15 × 25 15 × 15 × 33  Boundary condition (P a) hydrophilic hydrophobic 5300 6500 8500 8700  7000 9100 9800 10100  Results  The results of the numerical simulation for the thinnest sample (110µm) are presented in the following figures (Figures 4.1 to 4.10). The first replicate of each simulation is presented here. For each case, five cross sections in x− and y− directions (x∗ , y ∗ =0.05, 0.25, 0.50, 0.75 and 1.00) and four cross sections in z− direction (z ∗ =0.25, 0.50, 0.75 and 1.00) are presented. In addition, time instances are chosen: the initial configuration (t∗ = 0.01), three instances correspond to the process of flow during invasion (t∗ =0.25, 0.50 and 0.75) and the final configuration corresponds to the breakthrough (t∗ = 1.00). Temporal presentation of water flow in the samples shows the evolution of flow until breakthrough is achieved. In addition, the effect of PTFE loading is investigated in the next section using the figures presented in this section.  114  (b) y−planes  (c) z−planes  Figure 4.1: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.01.  (a) x−planes  4.1. Results  115  (b) y−planes  (c) z−planes  Figure 4.2: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.25.  (a) x−planes  4.1. Results  116  (b) y−planes  (c) z−planes  Figure 4.3: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.50.  (a) x−planes  4.1. Results  117  (b) y−planes  (c) z−planes  Figure 4.4: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 0.75.  (a) x−planes  4.1. Results  118  (b) y−planes  (c) z−planes  Figure 4.5: Saturation contours (sample thickness=110µm, PTFE loading=0%) at t∗ = 1.00.  (a) x−planes  4.1. Results  119  (b) y−planes  (c) z−planes  Figure 4.6: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.01.  (a) x−planes  4.1. Results  120  (b) y−planes  (c) z−planes  Figure 4.7: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.25.  (a) x−planes  4.1. Results  121  (b) y−planes  (c) z−planes  Figure 4.8: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.50.  (a) x−planes  4.1. Results  122  (b) y−planes  (c) z−planes  Figure 4.9: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 0.75.  (a) x−planes  4.1. Results  123  (b) y−planes  (c) z−planes  Figure 4.10: Saturation contours (sample thickness=110µm, PTFE loading=40%) at t∗ = 1.00.  (a) x−planes  4.1. Results  124  4.2. Discussion  4.2 4.2.1  Discussion Transient analysis of fluid flow  Fluid flow in a porous medium is a complex phenomenon due to the structural complexity associated with the medium. It is hard to assume a steady-state flow inside a medium since the boundary of the flow domain changes depending on its hydrodynamics. Even if the macroscopic properties of flow remain constant with time, there would always be changes at microscopic scale. Due to this complexity, it is required to define appropriate time scale and length scale of flow prior to study. In the case of fluid flow inside the GDL, breakthrough time can be used as a time scale. Thus, the phenomenon of fluid flow through the medium can be divided into two time intervals: before breakthrough and after breakthrough. The most challenging part is the study of flow before breakthrough, since the flow is transient and the hydrodynamics of flow is changing with time. Due to very low rates of flow through the GDL, once a sample-spanning cluster of liquid water is formed inside the GDL, the cluster can be seen as a permanent passage for delivering liquid water from the catalyst layer to the gas channel. By adopting the developed pore-network model, the transient flow of fluid through the GDL is investigated in this section. The evolution of liquid water inside the medium and the effect of hydrophobic content are studied. Figure 4.11 compares the evolution of liquid water through the GDL, for the hydrophilic and hydrophobic samples (0% and 40% PTFE loadings, respectively) at mid-plane y ∗ = 0.5. The corresponding time for each slice is normalized by the total time of breakthrough of the sample, i.e. t∗ = t/tBT . For the hydrophilic sample (4.11(a)), water flows through the sample quite uniformly. The frontier of invading fluid is not distinguishable reflecting a stable displacement mechanism. In addition, the local saturation is reducing from the catalyst layer (Scell = 1.0) up to the gas channel (Scell ≈ 0.0). As the breakthrough approaches, four possible spots are visible for water to find its way out (see Figure 4.11(a) for t∗ = 0.75). The conditions are quite different for the hydrophobic sample (Figure 4.11(b)). The imbibition process in the hydrophobic sample is not as uniform as the hydrophilic sample. Although the capillary numbers (Ca) are the same for both samples, the configuration of liquid water inside the hydrophobic sample is close to the viscous-fingering regime compared to the stable displacement for the hydrophilic sample. Thus, only one water frontier is distinguished which is indicated in Figure 4.11(b) at t∗ = 0.75. The viscous-fingering regime 125  4.2. Discussion  (a) Hydrophilic sample at y ∗ = 0.50  (b) Hydrophobic sample at y ∗ = 0.50  Figure 4.11: Effect of hydrophobic fraction on flow configuration at breakthrough. of liquid water configuration through the hydrophobic sample leaves the majority area of the sample open to the reactant gases. This is confirmed by existence of low saturated regions (blue areas) in the same time steps (t∗ ) for the hydrophobic and hydrophilic samples. The saturation profiles shown in Figure 4.12 compare total water content in each sample. Note that this is the total saturation along the z−axis. The water content was integrated at each x − y plane, and divided by the total pore volume of the same plane: Sk =  i i  j j  Vliquid,k Vpore,k  (4.1)  The difference between the saturation profile of the hydrophilic and hydrophobic samples is negligible at first (see the profiles at t∗ = 0.0 and t∗ = 0.25). As water flows into the samples, it covers the majority of the hydrophilic sample while the saturation of hydrophobic sample increases slowly 126  4.2. Discussion due to the viscous displacement regime. This leads to a higher saturation for the hydrophilic sample compared to the hydrophobic sample at the time of breakthrough (t∗ = 1.00).  4.2.2  Analysis of breakthrough  The stable displacement regime increases the chance of breakthrough for the hydrophilic sample compared to the hydrophobic sample. Thus, breakthrough is achieved earlier for the case of hydrophilic sample (as confirmed by the experimental time of breakthrough presented in Chapter 2). Table 4.3 summarized the normalized time of breakthrough obtained from the numerical modeling for three replicates. Figure 4.13 presents these values Table 4.3: Normalized values of time of breakthrough obtained from numerical model; For each case the average value is provided (second row) with the value of experimental breakthrough pressure used as boundary condition (in brackets). Thickness (µm)  Time (Pressure) hydrophilic  hydrophobic  110  0.5033 0.4824 0.4853 0.4904 (5300P a)  0.5979 0.5975 0.5165 0.5706 (7000P a)  190  0.5574 0.6064 0.5454 0.5697 (6500P a)  0.7408 0.7359 0.7340 0.7369 (9100P a)  280  0.8249 0.7063 0.7359 0.7557 (8500P a)  0.8656 0.9380 0.8773 0.8936 (9800P a)  370  0.7038 0.7836 0.7780 0.7551 (8700P a)  0.9909 0.9471 1.0000 0.9793(10100P a)  against the thickness of the GDL. As it is shown the time of breakthrough is larger for the hydrophobic samples compared to the hydrophilic one for all the thicknesses. The increase in the breakthrough time is not proportional to the increase in the sample thickness. In addition, the difference between the values for two thick samples (280µm and 370µm) is marginal for the hydrophilic samples. Figure 4.14 compares the results of the numerical model with experimental results. Except for the sample with the thickness of 280µm, the trend predicted by the numerical scheme is in good agreement with the experimental measurements. Considering the experimental and numerical errors, the difference is small. 127  4.2. Discussion The next step is to determine the flow rate and saturation in the hydrophilic and hydrophobic samples after the verification of the model based on the time of breakthrough. In essence, the saturation of the GDL sample at any time (t) can be estimated by integrating instantaneous flow rate Q(t): t  S(t) ∝  Q(t)dt  (4.2)  0  One can define the average flow rates for each sample as ˜ hydrophobic = Q ˜ hydrophilic = Q  Thydrophobic 0  Qhydrophobic (t)dt  Thydrophobic Thydrophilic 0  Qhydrophilic (t)dt  Thydrophilic  (4.3) (4.4)  where Thydrophobic and Thydrophilic are the time of breakthrough for the hydrophobic and hydrophilic samples, respectively. As a result, the saturation of each sample at breakthrough will be ˜ hydrophobic × Thydrophobic Shydrophobic ∝ Q ˜ hydrophilic × Thydrophilic Shydrophilic ∝ Q  (4.5) (4.6)  Numerical results presented in Figure 4.14, as well as experimental evidence ([LNK10],[PP09],[GFI+ 06], [GIFP08]) suggest that saturation in the hydrophilic sample must be more than that of the hydrophobic sample (Shydrophilic > Shydrophobic ). Considering the fact that Thydrophilic < Thydrophobic (shown both experimentally and numerically), the equations 4.5 and 4.6 suggest that the average injection flow rate into the hydrophobic sample should be lower than the injection flow rate of the hydrophilic sample ˜ hydrophobic < Q ˜ hydrophilic Q  (4.7)  Figures 4.15 presents the comparison of instant flow rates for the hydrophobic and hydrophilic samples. Each figure corresponds to a particular sample thickness. In each thickness, the initial flow rates (Q0 ) for the hydrophobic and hydrophilic samples are the same. As water flows through the sample, the instantaneous flow rate (Q(t)) continuously decreases for both samples. However, the flow rate is lower for the hydrophobic sample although the pressure boundary condition is higher. This confirms the findings presented in Equation 4.7. The flow rates obtained numerically can be used to calculate the saturation in each sample at breakthrough. In essence, the total water saturation 128  4.2. Discussion for each sample can be obtained by integrating the water content within each cell: i j k Vliquid S= (4.8) i j k Vpore The saturation versus thickness is plotted in Figure 4.16. For each thickness, the hydrophilic sample contains more water. As the sample thickness increases, the effect of the hydrophobic content on saturation becomes insignificant, i.e. for the thickest sample the hydrophobicity does not have any effect on the total saturation. Thus, treating thick GDL samples with a hydrophobic agent has little effect in terms of saturation and water management. For hydrophilic samples, on the other hand, the difference in saturation values are marginal. This suggests that a thicker GDL sample may be employed for typical fuel cell operation at the cost of little increase in the level of saturation at breakthrough.  129  4.2. Discussion  (a) at t∗ = 0.0 and t∗ = 0.25  (b) at t∗ = 0.50  (c) at t∗ = 0.75  (d) at t∗ = 1.00  Figure 4.12: Numerical comparison of through-plane saturation profiles at different time steps for hydrophilic (solid line) and hydrophobic (dash line) samples.  130  4.2. Discussion  Figure 4.13: Analysis of error variance for the numerical results.  131  4.2. Discussion  Figure 4.14: Comparison of the numerical and experimental results for time of breakthrough (the shadow regions show the error interval for each case).  132  4.2. Discussion  (a) Thickness of 110µm  (b) Thickness of 190µm  (c) Thickness of 280µm  (d) Thickness of 370µm  Figure 4.15: Comparison of flow rates for hydrophobic and hydrophilic samples.  133  4.2. Discussion  Figure 4.16: Effect of sample thickness on total saturation.  134  Chapter 5  Conclusions and future work Overall analysis and integration of the research and conclusions of the dissertation In this research, a comprehensive study was performed to understand the complex phenomena of liquid water flow through the gas diffusion layers (GDLs) of fuel cells. The study includes both experimental measurements and numerical simulation to establish an accurate model for flow of liquid water. A visualization technique (fluorescence microscopy) was adopted and the evolution of liquid water flow inside the porous GDL was thoroughly investigated. The main characteristics of the flow were extracted and derived from the visualization study. The characteristics were then used to propose a model for water flow in the medium at the pore scale. A pore–network approach was adopted as for the numerical simulation of fluid flow. The developed model was then applied to several experimental cases, and the results of the simulations and the experimental measurements were compared. Using the developed model, the effects of GDL properties on the water flow and performance of a fuel cell were studied. In the following section, the specific conclusions and contributions of the research are presented.  Conclusions regarding hypotheses of the dissertation The analysis of liquid water visualization in the porous medium of a gas diffusion layer (GDL) shows that treatment of the medium with the hydrophobic agent (PTFE) modifies the flow pattern at the microscopic level. The main differences in the flow pattern between the untreated and treated media are: − Invasion pattern refers to the initial development of liquid flow in the porous GDL. Although the surface of treated and untreated samples are equally exposed to the same injection flow rate, the visualization results reveal that liquid water flows into the majority of the available pores on the boundary of the untreated GDL while for the treated GDL, only a handful of boundary pores are invaded. Re135  Chapter 5. Conclusions and future work duction of pore sizes and higher contact angle of the pores due to the treatment are the main reasons for the difference between the invasion patterns. − Progression pattern refers to evolution of liquid water patterns over time. While water flows through the GDL, it develops certain pathways toward the gas flow field. For the untreated GDL, several branches segregate from the initial pathways; while for the treated GDL, the original pathways developed are extended toward the other side of the medium with minimum possible branches. Thus, the pattern of water flow in the untreated GDL is shifted toward stable displacement. For the treated GDL, viscous fingering best describes the flow pattern and final configuration of liquid water in the medium. − Pore–filling pattern refers to the mechanism by which the liquid water fills the void space in the medium. It was shown that once a pathway is developed through the GDL, it starts feeding the pores connected to it. In a treated GDL at each time step, only a limited number of accessible pores are invaded by water and the process continues until water completely fills these pores. In contrast, for an untreated GDL, several pores are simultaneously being invaded by developed liquid water pathways. Moreover, for the untreated GDL the invading pore does not necessarily need to be full prior to the invasion of the adjacent pores. In essence, water invades the adjacent pores once its pressure overcomes the pore threshold pressure. The experimental setup developed in this research enabled the measurement of the time and pressure at the breakthrough. The results show that: − the breakthrough pressure follows the same pattern as the time of breakthrough. − the capillary pressure does not necessarily increase as the thickness of the porous medium increases. − different trends for the change of the capillary pressure as a function of the hydrophobic content of the GDL exist. For low flow rate, the maximum pressure is obtained for the PTFE loading of 40wt%; while for the higher flow rate, the capillary pressure monotonically increases as the hydrophobicity increases. − the assumption of constant flow rate for injecting water into the sample may overestimate the saturation. 136  Chapter 5. Conclusions and future work The analysis of liquid water patterns through the GDL, which were obtained from the visualization technique, was used to verify the pore–network algorithm developed in this research. In essence, the flow characteristics determined in the experimental visualization were adopted at microscopic pore–level to propose an accurate model for displacement and flow of liquid water in the pore network. The results of the simulations confirm that: − for the treated GDL, the local pressure of liquid increases to overcome the capillary pressure threshold. This, in turn, leads to a higher local saturation (pore saturation) for the treated GDL; while the overall medium saturation remains lower compared to the untreated GDL. − higher overall saturation for untreated medium means higher probability for breakthrough. Thus, the time required for breakthrough is shorter for the untreated GDL. − water flow rates into the treated sample are smaller than those into the untreated one; while the pressure boundary condition is higher. − the simulations suggest that treatment of a thin GDL with hydrophobic agent improves the cell performance by lowering the overall saturation. As the sample thickness increases, however, the effect of treatment becomes insignificant. Eventually, saturation of the treated and untreated GDLs become the same for the thickest sample modeled in this research.  Strengths and limitations of the dissertation research The method of fluorescence microscopy used in this research is categorized as a direct visualization method; thus, the results are directly used to develop the concepts regarding the liquid water flow through the porous medium. This removes any additional data processing procedure to interpret the experimental data for further analysis. The numerical scheme proposed in this study is very efficient. It provides data at the pore level. Basically in a pore–network approach, the complexity of the problem involved in fluid flow in porous media is broken down into two main sub-problems: reconstruction of the porous medium with a network of pores connected to each other via throats, and implementation of simple physical rules for displacing the fluid. The main strength of the research is direct implementation of the experimental findings and concepts into the process of model development. This leads to an accurate model for the water flow through the porous medium. 137  Chapter 5. Conclusions and future work PEM fuel cells work at typical temperatures in the range of 70 − 80◦ C. However, all the experiments of this research were conducted at the room temperature (20◦ C). The variation of temperature definitely has a great impact as it is directly contributing in evaporation/condensation occurring through the GDL. The numerical model developed in this research, also, does not include any condensation/evaporation occur in a fuel cell due to higher level of working temperature.  Potential applications of the research findings The experimental findings can be extensively used in future model developments in the field of fluid flow in micro-porous structures. The data collected for the breakthrough pressure and time for various GDL samples can be used as the boundary condition and validation criteria of any fluid flow model. The pore–network model developed in this research is a tool for studying and designing future GDLs through which flow of liquid water is controlled and managed to prevent flooding. The numerical scheme developed in this thesis is not limited to water flow through the GDL, but the concepts can be extensively employed and used to develop models for fluid flow in thin porous media.  Analysis of possible future research The challenging part of ex–situ study of flooding phenomena, and water flow through the GDL, is handling the very low flow rate. The determination of the flow rate is necessary to accurately estimate saturation of the medium in each instance. However, direct measurement of the exact flow rate is not a trivial task. Thus, flow rate estimation and the methods for its measurement have great practical application in this area. The working temperature of a PEM fuel cell is higher than the room temperature used in this research. The variation of temperature definitely has a great impact as it is directly contributing in evaporation/condensation occurring through the GDL. Thus, studying the effects of the temperature variation on water flow through GDL is recommended. Although the analysis of the image obtained in this research from fluorescence microscopy images differs considerably from the original method proposed by Litster et al. ([LSD06]), the image analysis can still be improved by considering different possible scenarios for fluid flow through the porous medium and finding the exact path of water flow.  138  Chapter 5. Conclusions and future work Following the enhancement in the image analysis technique, the assumptions of the model concepts developed for the fluid flow through the porous media can be significantly improved. This improvement will be enabled by studying different possible scenarios of fluid flow in the medium obtained by the image analysis described above. The proposed numerical scheme only considered water as the dynamic phase with variable hydrodynamic properties throughout the pore network. In other words, the gas phase is static with constant pressure for the whole domain. The model proposed in this research can be modified to include variable hydrodynamic properties for the gas phase. Also, the effect of temperature gradient on the gas phase properties and consequently on the fluid flow can be studied. Moreover, the interaction of both phases due to the temperature variation (evaporation/condensation) followed by the inclusion of second phase (gas) can be investigated and studied. The possible future research can be summarized as follows: − direct measurement of the flow rate − image enhancement for improvement of the model − including the hydrodynamic properties of the gas phase, the effect of temperature variation In addition, the experimental setup can be modified to study saturation after the time of breakthrough. The MPL can also be included as another pore-network structure attached to the GDL pore-network to study its effect on the GDL overall saturation.  139  Bibliography [AIA98] AIAA. Guide for the verification and validation of computational fluid dynamics simulations. AIAA G-077-1998, 1998. → pages 74 [Baz09] A. Bazylak. Liquid water visualization in pem fuel cells: A review. International Journal of Hydrogen Energy, 34(9):3845 – 3857, 2009. → pages 18 [BHDS08] A. Bazylak, J. Heinrich, N. Djilali, and D. Sinton. Liquid water transport between graphite paper and a solid surface. Journal of Power Sources, 185(2):1147–1153, 2008. → pages 25, 26, 29, 31, 32, 48, 51 [BK91] M. Blunt and P. King. Relative permeabilities from two-and three-dimensional pore-scale network modelling. Transport in Porous Media, 6(4):407–433, 1991. → pages 11 [BKS92] M. Blunt, M. J. King, and H. Scher. Simulation and theory of two-phase flow in porous media. Physical Review A, 46(12):7680–7699, 1992. → pages 11 [BKTO05] J. Borrelli, S.G. Kandlikar, T. Trabold, and J. Owejan. Water transport visualization and two-phase pressure drop measurements in a simulated pemfc cathode minichannel. In Proceedings of 3rd International Conference on Microchannels and Minichannels, Toronto, Canada, 2005. → pages 21, 22, 23, 25 [BLA+ 99] RJ Bellows, MY Lin, M. Arif, AK Thompson, and D. Jacobson. Neutron imaging technique for in situ measurement of water transport gradients within nafion in polymer electrolyte fuel cells. Journal of the Electrochemical Society, 146:1099, 1999. → pages 19  140  Bibliography [Blu97] M. Blunt. Effects of heterogeneity and wetting on relative permeability using pore level modeling. SPE Journal, 2(1):70– 87, 1997. → pages 15 [Blu98] M.J. Blunt. Physically-based network modeling of multiphase flow in intermediate-wet porous media. Journal of Petroleum Science and Engineering, 20(3):117–125, 1998. → pages 11, 15 [Blu01] M.J. Blunt. Flow in porous media–pore-network models and multiphase flow. Current opinion in colloid & interface science, 6(3):197–207, 2001. → pages 11, 61 [BNB+ 05] Jay Benziger, James Nehlsen, David Blackwell, Tom Brennan, and Johannah Itescu. Water flow in the gas diffusion layer of pem fuel cells. Journal of Membrane Science, 261(1):98–106, 2005. → pages 52, 53 [BSD08] A. Bazylak, D. Sinton, and N. Djilali. Dynamic water transport and droplet emergence in pemfc gas diffusion layers. Journal of Power Sources, 176(1):240–246, 2008. → pages 25, 26, 31, 32, 48, 49, 51 [BSH+ 99] JD Bobyn, GJ Stackpool, SA Hacking, M. Tanzer, and JJ Krygier. Characteristics of bone ingrowth and interface mechanics of a new porous tantalum biomaterial. Journal of Bone & Joint Surgery, British Volume, 81(5):907–914, 1999. → pages 8 [BSL06] R Byron Bird, Warren E Stewart, and Edwin N Lightfoot. Transport phenomena. Wiley, 2006. → pages 59 [BSLD07] A. Bazylak, D. Sinton, Z.S. Liu, and N. Djilali. Effect of compression on liquid water transport and microstructure of pemfc gas diffusion layers. Journal of Power Sources, 163(2):784–792, 2007. → pages 24, 26, 29, 31, 32, 48, 51 [BTM09] D. Blessent, R. Therrien, and K. MacQuarrie. Coupling geological and numerical models to simulate groundwater flow and contaminant transport in fractured media. Computers & Geosciences, 35(9):1897–1906, 2009. → pages 8 [CPQ+ 08] O. Chapuis, M. Prat, M. Quintard, E. Chane-Kane, O. Guillot, and N. Mayer. Two-phase flow and evaporation in model 141  Bibliography fibrous media Application to the gas diffusion layer of PEM fuel cells. Journal of Power Sources, 178(1):258–268, 2008. → pages 12, 59, 60, 68, 70 [EHPP11] M. El Hannach, J. Pauchet, and M. Prat. Pore network modeling: Application to multiphase transport inside the cathode catalyst layer of proton exchange membrane fuel cell. Electrochimica Acta, 2011. → pages 14, 15, 59, 68, 70 [FBW06] K.W. Feindel, S.H. Bergens, and R.E. Wasylishen. Insights into the distribution of water in a self-humidifying h2/o2 proton-exchange membrane fuel cell using 1h nmr microscopy. Journal of the American Chemical Society, 128(43):14192– 14199, 2006. → pages 18 [FCS10] J.D. Fairweather, P. Cheung, and D.T. Schwartz. The effects of wetproofing on the capillary properties of proton exchange membrane fuel cell gas diffusion layers. Journal of Power Sources, 195(3):787–793, 2010. → pages 14 [FCSPS07] Joseph D Fairweather, Perry Cheung, Jean St-Pierre, and Daniel T Schwartz. A microfluidic approach for measuring capillary pressure in pemfc gas diffusion layers. Electrochemistry communications, 9(9):2340–2345, 2007. → pages 52, 53 [FLS+ 04] K.W. Feindel, L.P.A. LaRocque, D. Starke, S.H. Bergens, and R.E. Wasylishen. In situ observations of water production and distribution in an operating h2/o2 pem fuel cell assembly using 1h nmr microscopy. Journal of the American Chemical Society, 126(37):11436–11437, 2004. → pages 18 [FSH08] B. R. Friess, M. Shahraeeni, and M. Hoorfar. Water management in pem fuel cells. In Proceedings of 2008 AIChE Conference, Philadelphia, USA, 2008. → pages 2 [GFI+ 06] J.T. Gostick, M.W. Fowler, M.A. Ioannidis, M.D. Pritzker, Y.M. Volfkovich, and A. Sakars. Capillary pressure and hydrophilic porosity in gas diffusion layers for polymer electrolyte fuel cells. Journal of Power Sources, 156(2):375–387, 2006. → pages 16, 112, 128 [GIFP07] J. T. Gostick, M. A. Ioannidis, M. W. Fowler, and M. D. Pritzker. Pore network modeling of fibrous gas diffusion layers 142  Bibliography for polymer electrolyte membrane fuel cells. Journal of Power Sources, 173(1):277–290, 2007. → pages 11, 59, 68, 70 [GIFP08] J.T. Gostick, M.A. Ioannidis, M.W. Fowler, and M.D. Pritzker. Direct measurement of the capillary pressure characteristics of water-air-gas diffusion layer systems for pem fuel cells. Electrochemistry Communications, 10(10):1520–1523, 2008. → pages 53, 112, 128 [GR93] P. A. Goode and T. S. Ramakrishnan. Momentum transfer across fluid-fluid interfaces in porous media: a network model. AIChE Journal, 39(7), 1993. → pages 11 [GW07] S. Ge and C.Y. Wang. Liquid water formation and transport in the pefc anode. Journal of the Electrochemical Society, 154:B998, 2007. → pages 21, 22, 25 [HMK+ 08] C. Hartnig, I. Manke, R. Kuhn, N. Kardjilov, J. Banhart, and W. Lehnert. Cross-sectional insight in the water evolution and transport in polymer electrolyte fuel cells. Applied Physics Letters, 92(13):134106–134106, 2008. → pages 21 [HMWH04] A. Hakenjos, H. Muenter, U. Wittstadt, and C. Hebling. A pem fuel cell for combined measurement of current and temperature distribution, and flow field flooding. Journal of Power Sources, 131(1):213–216, 2004. → pages 21, 22, 25 [HSC+ 06] MA Hickner, NP Siegel, KS Chen, DN McBrayer, DS Hussey, DL Jacobson, and M. Arif. Real-time imaging of liquid water in an operating proton exchange membrane fuel cell. Journal of the Electrochemical Society, 153:A902, 2006. → pages 20 [KKNK11] J.H. Kang, K.N. Kim, J.H. Nam, and C.J. Kim. Visualization experiments and pore network simulations for invasionpercolation transport of a non-wetting fluid through multiple porous layers. International Journal of Hydrogen Energy, 2011. → pages 15 [KM09] S. Kim and MM Mench. Investigation of temperature-driven water transport in polymer electrolyte fuel cell: phase-changeinduced flow. Journal of the Electrochemical Society, 156:B353, 2009. → pages 20 143  Bibliography [KPP11] S.P. Kuttanikkad, M. Prat, and J. Pauchet. Pore-network simulations of two-phase flow in a thin porous layer of mixed wettability: Application to water transport in gas diffusion layers of proton exchange membrane fuel cells. Journal of Power Sources, 196(3):1145–1155, 2011. → pages 13, 15, 59, 68, 70 [KSM06] EC Kumbur, KV Sharp, and MM Mench. Liquid droplet behavior and instability in a polymer electrolyte fuel cell flow channel. Journal of Power Sources, 161(1):333–345, 2006. → pages 21, 23, 25 [KWKB08] E. Kimball, T. Whitaker, Y.G. Kevrekidis, and J.B. Benziger. Drops, slugs, and flooding in polymer electrolyte membrane fuel cells. AIChE Journal, 54(5):1313–1332, 2008. → pages 21, 24, 25 [LDF03] J. Larminie, A. Dicks, and Knovel (Firm). Fuel cell systems explained. 2003. → pages 4, 6, 7, 32 [LKR+ 09] Z. Lu, SG Kandlikar, C. Rath, M. Grimm, W. Domigan, AD White, M. Hardbarger, JP Owejan, and TA Trabold. Water management studies in pem fuel cells, part ii: Ex situ investigation of flow maldistribution, pressure drop and two-phase flow pattern in gas channels. international journal of hydrogen energy, 34(8):3445–3456, 2009. → pages 21 [LLK+ 08] S.J. Lee, N.Y. Lim, S. Kim, G.G. Park, and C.S. Kim. X-ray imaging of water distribution in a polymer electrolyte fuel cell. Journal of Power Sources, 185(2):867–870, 2008. → pages 21 [LMW06] Z. Liu, Z. Mao, and C. Wang. A two dimensional partial flooding model for pemfc. Journal of power sources, 158(2):1229– 1239, 2006. → pages 21 [LNK10] K.J. Lee, J.H. Nam, and C.J. Kim. Steady saturation distribution in hydrophobic gas-diffusion layers of polymer electrolyte membrane fuel cells: A pore-network study. Journal of Power Sources, 195(1):130–141, 2010. → pages 15, 112, 128 [LSD06] S. Litster, D. Sinton, and N. Djilali. Ex situ visualization of liquid water transport in PEM fuel cell gas diffusion layers. Journal of Power Sources, 154(1):95–105, 2006. → pages 24, 26, 29, 31, 32, 48, 49, 51, 138 144  Bibliography [MBD07] B. Markicevic, A. Bazylak, and N. Djilali. Determination of transport parameters for multiphase flow in porous gas diffusion electrodes using a capillary network model. Journal of Power Sources, 171(2):706–717, 2007. → pages 16 [MD11] B. Markicevic and N. Djilali. Analysis of liquid water transport in fuel cell gas diffusion media using two-mobile phase pore network simulations. Journal of Power Sources, 196(5):2725– 2734, 2011. → pages 12, 13, 59, 68, 70 [MGP96] R. Mosdale, G. Gebel, and M. Pineri. Water profile determination in a running proton exchange membrane fuel cell using small-angle neutron scattering. Journal of membrane science, 118(2):269–277, 1996. → pages 18 [MHG+ 07] I. Manke, C. Hartnig, M. Grunerbel, W. Lehnert, N. Kardjilov, A. Haibel, A. Hilger, J. Banhart, and H. Riesemeier. Investigation of water evolution and transport in fuel cells with high resolution synchrotron x-ray radiography. Applied Physics Letters, 90(17):174105–174105, 2007. → pages 20 [MRFL03] M. Mathias, J. Roth, J. Fleming, and W. Lehnert. Diffusion media materials and characterization. Handbook of fuel cellsFundamentals, technology and applications, 3:517–537, 2003. → pages 63 [NK03] J. H. Nam and M. Kaviany. Effective diffusivity and watersaturation distribution in single-and two-layer PEMFC diffusion medium. International Journal of Heat and Mass Transfer, 46(24):4595–4611, 2003. → pages 11, 16, 24, 59, 68 [OTJ+ 06] JP Owejan, TA Trabold, DL Jacobson, DR Baker, DS Hussey, and M. Arif. In situ investigation of water transport in an operating pem fuel cell using neutron radiography: part 2–transient water accumulation in an interdigitated cathode flow field. International journal of heat and mass transfer, 49(25):4721– 4731, 2006. → pages 18 [OTJ+ 07] JP Owejan, TA Trabold, DL Jacobson, M. Arif, and SG Kandlikar. Effects of flow field and diffusion layer properties on water accumulation in a pem fuel cell. International Journal of Hydrogen Energy, 32(17):4489–4502, 2007. → pages 20 145  Bibliography [PHC+ 05] N. Pekula, K. Heller, PA Chuang, A. Turhan, MM Mench, ¨ u. Study of water distribution and JS Brenizer, and K. Unl¨ transport in a polymer electrolyte fuel cell using neutron imaging. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 542(1):134–141, 2005. → pages 18 [PLP08] S. Park, J. W. Lee, and B. N. Popov. Effect of PTFE content in microporous layer on water management in PEM fuel cells. Journal of Power Sources, 177(2):457–463, 2008. → pages 6, 105 [PLP12] S. Park, J.W. Lee, and B.N. Popov. A review of gas diffusion layer in pem fuel cells: Materials and designs. International Journal of Hydrogen Energy, 2012. → pages 105 [PM08] J. B. Pawley and B. R. Masters. Handbook of biological confocal microscopy. Journal of Biomedical Optics, 13:029902, 2008. → pages 30 [PMS12] L.M. Pant, S.K. Mitra, and M. Secanell. Absolute permeability and knudsen diffusivity measurements in pemfc gas diffusion layers and micro porous layers. Journal of Power Sources, 2012. → pages 16 [PP09] S. Park and B.N. Popov. Effect of hydrophobicity and pore geometry in cathode gdl on pem fuel cell performance. Electrochimica Acta, 54(12):3473–3479, 2009. → pages 112, 128 [PPSK11] J. Pauchet, M. Prat, P. Schott, and S.P. Kuttanikkad. Performance loss of proton exchange membrane fuel cell due to hydrophobicity loss in gas diffusion layer: Analysis by multiscale approach combining pore network and performance modelling. International Journal of Hydrogen Energy, 2011. → pages 15, 59, 68 [PW05] U. Pasaogullari and C.Y. Wang. Two-phase modeling and flooding prediction of polymer electrolyte fuel cells. Journal of the Electrochemical Society, 152:A380, 2005. → pages 21, 24, 26 [Sah93] M. Sahimi. Flow phenomena in rocks: from continuum models to fractals, percolation, cellular automata, and simulated 146  Bibliography annealing. Reviews of modern physics, 65(4):1393, 1993. → pages 10 [SH11] M. Shahraeeni and M. Hoorfar. Effect of gas diffusion layer properties on the time of breakthrough. Journal of Power Sources, 196(1):5918–5921, 2011. → pages 33 [Sha81] D.O. Shah. Surface phenomena in enhanced oil recovery. Plenum Pub Corp, 1981. → pages 8 [SHW06] P.K. Sinha, P. Halleck, and C.Y. Wang. Quantification of liquid water saturation in a pem fuel cell diffusion medium using x-ray microtomography. Electrochemical and solid-state letters, 9:A344, 2006. → pages 20 [SJAW04] R. Satija, DL Jacobson, M. Arif, and SA Werner. In situ neutron imaging technique for evaluation of water management systems in operating pem fuel cells. Journal of Power Sources, 129(2):238–245, 2004. → pages 18 [SMS+ 08] J.B. Siegel, D.A. McKay, A.G. Stefanopoulou, D.S. Hussey, and D.L. Jacobson. Measurement of liquid water accumulation in a pemfc with dead-ended anode. Journal of The Electrochemical Society, 155:B1168, 2008. → pages 19 [SNY+ 05] K. Sugiura, M. Nakata, T. Yodo, Y. Nishiguchi, M. Yamauchi, and Y. Itoh. Evaluation of a cathode gas channel with a water absorption layer/waste channel in a pefc by using visualization technique. Journal of power sources, 145(2):526–533, 2005. → pages 21, 22, 25 [SPA07] D. Spernjak, A.K. Prasad, and S.G. Advani. Experimental investigation of liquid water formation and transport in a transparent single-serpentine pem fuel cell. Journal of Power Sources, 170(2):334–344, 2007. → pages 21, 23, 25 [SW07] P. K. Sinha and C. Y. Wang. Pore-network modeling of liquid water transport in gas diffusion layer of a polymer electrolyte fuel cell. Electrochimica Acta, 52(28):7936–7945, 2007. → pages 12, 59, 60, 68, 69, 70 [SW08] P.K. Sinha and C.Y. Wang. Liquid water transport in a mixedwet gas diffusion layer of a polymer electrolyte fuel cell. Chemical engineering science, 63(4):1081–1091, 2008. → pages 12 147  Bibliography [SZG91] T.E. Springer, TA Zawodzinski, and S. Gottesfeld. Polymer electrolyte fuel cell model. J. electrochem. Soc, 138(8):2334– 2342, 1991. → pages 20 [THBM06] A. Turhan, K. Heller, JS Brenizer, and MM Mench. Quantification of liquid water accumulation and distribution in a polymer electrolyte fuel cell using neutron imaging. Journal of Power Sources, 160(2):1195–1203, 2006. → pages 19 [TOG+ 06] A. Theodorakakos, T. Ous, M. Gavaises, JM Nouri, N. Nikolopoulos, and H. Yanagihara. Dynamics of water droplets detached from porous surfaces of relevance to pem fuel cells. Journal of colloid and interface science, 300(2):673– 687, 2006. → pages 21, 23 [TOJ+ 06] TA Trabold, JP Owejan, DL Jacobson, M. Arif, and PR 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 journal of heat and mass transfer, 49(25):4712–4720, 2006. → pages 18 [TPH03] K. T¨ uber, D. P´ocza, and C. Hebling. Visualization of water buildup in the cathode of a transparent pem fuel cell. Journal of Power Sources, 124(2):403–414, 2003. → pages 21, 22, 25 [TTH04] S. Tsushima, K. Teranishi, and S. Hirai. Magnetic resonance imaging of the water distribution within a polymer electrolyte membrane in fuel cells. Electrochemical and solid-state letters, 7:A269, 2004. → pages 17 [TTH05] K. Teranishi, S. Tsushima, and S. Hirai. Study of the effect of membrane thickness on the performance of polymer electrolyte fuel cells by water distribution in a membrane. Electrochemical and Solid-State Letters, 8:A281, 2005. → pages 17 [WDN04] A.Z. Weber, R.M. Darling, and J. Newman. Modeling twophase behavior in pefcs. Journal of the Electrochemical Society, 151:A1715, 2004. → pages 65 [Whi67] Stephen Whitaker. Diffusion and dispersion in porous media. AIChE Journal, 13(3), 1967. → pages 16 148  Bibliography [WLZW12a] R. Wu, Q. Liao, X. Zhu, and H. Wang. Impacts of the mixed wettability on liquid water and reactant gas transport through the gas diffusion layer of proton exchange membrane fuel cells. International Journal of Heat and Mass Transfer, 2012. → pages 16 [WLZW12b] R. Wu, Q. Liao, X. Zhu, and H. Wang. Pore network modeling of cathode catalyst layer of proton exchange membrane fuel cell. International Journal of Hydrogen Energy, 2012. → pages 16 [WSHL06] F.B. Weng, A. Su, C.Y. Hsu, and C.Y. Lee. Study of water-flooding behaviour in cathode channel of a transparent proton-exchange membrane fuel cell. Journal of power sources, 157(2):674–680, 2006. → pages 21, 24, 25 [YZLW04] XG Yang, FY Zhang, AL Lubawy, and CY Wang. Visualization of liquid water transport in a pefc. Electrochemical and Solid-State Letters, 7:A408, 2004. → pages 21, 22, 25 [ZKS+ 06] J. Zhang, D. Kramer, R. Shimoi, Y. Ono, E. Lehmann, A. Wokaun, K. Shinohara, and G.G. Scherer. In situ diagnostic of two-phase flow phenomena in polymer electrolyte fuel cells by neutron imaging:: Part b. material variations. Electrochimica Acta, 51(13):2715–2727, 2006. → pages 19 [ZTD+ 93] Jr. Zawodzinski, A. Thomas, Derouin, Charles, Radzinski, Susan, Sherman, J. Ruth, Smith, T. Van, Springer, E. Thomas, Gottesfeld, and Shimshon. Water uptake by and transport through nafion (r) 117 membranes. Journal of the Electrochemical Society, 140:1041–1047, 1993. → pages 17, 19  149  Appendices  150  Appendix A  Code structure The numerical algorithm of this research is developed using C++ programming language, with an object-oriented approach. In this approach, the network data are stored in a 3D array of pores. Each pore is considered as an object with certain attributes. These attributes are defined in “classes”. Three main classes are defined: 1) class fluidC , 2) class throatC, and 3) class pore. The fluidC represents all the fluids and their properties within the simulation, including density, viscosity, contact angle, surface tension, phase pressure and temperature. The definition and attributes of class fluidC with two samples of public functions are presented in the following: class fluidC { private : double d e n s i t y ; double v i s c o s i t y ; double c o n t a c t a n g l e ; double t e n s i o n ; double p r e s s u r e ; double t e m p e r a t u r e ; double mass ; double volume ; public : //SET void s e t d e n s i t y ( double den ) { d e n s i t y=den ; }; ... //GET double g e t d e n s i t y ( ) { return d e n s i t y ; }; };  In the throatC class, the properties of a generic throat is defined including its diameter, length and the flow rates associated with each phase: c l a s s throatC { private : double d i a m e t e r ;  151  Appendix A. Code structure double l e n g t h ; double q [ n u m b e r o f p h a s e s ] ; public : //SET void s e t d i a m a t e r ( double d i a ) { d i a m e t e r=d i a ; }; void s e t l e n g t h ( double l e n ) { l e n g t h=l e n ; }; void s e t q ( double qq , i n t i ) { q [ i ]=qq ; }; //GET double g e t d i a m a t e r ( ) { return d i a m e t e r ; }; double g e t l e n g t h ( ) { return l e n g t h ; }; double g e t q ( i n t i ) { return q [ i ] ; }; };  The third generic class is the class pore, which only has the diameter as the private variable. This class includes an array of “fluid” with its dimension equals to the number of phases being simulated. In the “throat” array six elements associated with six throats connected to the pore are defined: class pore { private : double d i a m e t e r ; public : fluidC f l u i d [ number of phases ] ; throatC t h r o a t [ 6 ] ; void s e t d i a m e t e r ( double d i a ) { d i a m e t e r=d i a ; }; double g e t d i a m e t e r ( ) { return d i a m e t e r ; }; // w r i t e a l l t h e e x e c u t i v e f u n c t i o n s i n s i d e t h e p o r e c l a s s . double p o r e v o l u m e ( ) { return ( PI ∗ ( pow ( g e t d i a m e t e r ( ) , 3 . 0 ) ) / 6 . 0 ) ; }; ...  152  Appendix A. Code structure //UPDATE VOLUME OF EACH PHASE IN EACH PORE FOR DT void update volume ( void ) { int i , j ; double tmp1 ; f o r ( i =0; i <n u m b e r o f p h a s e s ; i ++){ tmp1 = 0 . 0 ; f o r ( j =0; j <6; j ++){ tmp1=tmp1+t h r o a t [ j ] . g e t q ( i ) ; }; f l u i d [ i ] . set volume ( f l u i d [ i ] . get volume ( ) }; }; };  Note that all the executive functions are embedded into the definition of the classes. In this way, the communication and the flow of information between various subroutines involved in the code will be strictly under controlled, which effectively makes the code readable and the debugging process efficient. In the declaration section, the main network and a temporary network are defined along with the index counters and the output files. The subroutines performing auxiliary calculations such as exporting the network data to an output file (WriteNet()), Generating the network(NetGen()), calculating the flow rates within the network’s throats (calc q()), copying a network into a new network (copyNet()), assigning the hydrodynamic properties to the generated geometry of the network (assign hydrodynamics(), and treatment of the network for random distribution of contact angle to the pores (treatment()), reading and loading a network from an existing file (readNet()) and saving the individual properties of the network to the files (writeNetDia(), writeNetPressure(), writePoreSaturation() and writeNetSaturation()) are also defined in this section. p o r e Mypore ; p o r e Net [IMAX ] [ JMAX ] [KMAX] , Net tmp [IMAX ] [ JMAX ] [KMAX] ; i n t INDEX I [IMAX ] ; i n t INDEX J [JMAX ] ; i n t INDEX K [KMAX] ; fstream f n e t ; o f s t r e a m p0 , s , ps ; void w r i t e N e t ( char f i l e n a m e [ ] ) { f n e t . open ( f i l e n a m e , i o s : : out | i o s : : b i n a r y ) ; f n e t . w r i t e ( reinterpret cast <char ∗>(&Net ) , s i z e o f ( p o r e ) ∗IMAX∗JMAX∗KMAX) ; f net . close ();  153  Appendix A. Code structure cout<<” \nNETWORK IS STORED IN : ”<<f i l e n a m e ; // This g e n e r a t e s t h e n e t w o r k and s t o r e s i t i n memory NET; // Note : I t doesn ’ t w r i t e i t i n t o t h e f i l e ! // Note : This p r o d u c e s a random n e t w o r k ! // Note : This o n l y g e n e r a t e s t h e g e o m e t r i c a l n e t w o r k . // Note : The hydrodynamic p r o p e r t i e s a r e n o t a s s i n g e d h e r e i n ! void NetGen ( ) { . . . } int i , j , k , l ; double tmp1 , tmp2 ; double rndd ; cout<<” \nGENERATING NEW NETWORK. . . ” ; s r a n d ( ( unsigned ) time ( 0 ) ) ; f o r ( i =0; i <IMAX; i ++) f o r ( j =0; j <JMAX; j ++) f o r ( k=0;k<KMAX; k++){ //1<PORE DIAMETER<2 // 0.5 <THROAT DIAMETER<1 ... cout<<” \nSYNCRONIZING NEW NETWORK. . . ” ; //SYNCRONIZING //ORIGIN i =0; j =0; k =0; Net [ i ] [ j ] [ k ] . t h r o a t [ 0 ] . s e t d i a m a t e r ( 0 . 0 ) ; Net [ i ] [ j ] [ k ] . t h r o a t [ 1 ] . s e t d i a m a t e r ( 0 . 0 ) ; Net [ i ] [ j ] [ k ] . t h r o a t [ 5 ] . s e t d i a m a t e r ( 0 . 0 ) ; ... //CORE f o r ( i =1; i <IMAX; i ++) f o r ( j =1; j <JMAX; j ++) f o r ( k=1;k<KMAX; k++){ Net [ i ] [ j ] [ k ] . t h r o a t [ 0 ] . s e t d i a m a t e r ( Net [ i ] [ j ] [ k − 1 ] . t h r o a t [ 2 ] . g e t d i a m a t e r ( ) ) ; Net [ i ] [ j ] [ k ] . t h r o a t [ 1 ] . s e t d i a m a t e r ( Net [ i ] [ j − 1 ] [ k ] . t h r o a t [ 3 ] . g e t d i a m a t e r ( ) ) ; Net [ i ] [ j ] [ k ] . t h r o a t [ 5 ] . s e t d i a m a t e r ( Net [ i − 1 ] [ j ] [ k ] . t h r o a t [ 4 ] . g e t d i a m a t e r ( ) ) ; }; }; }; // This c a l c u l a t e s t h e f l o w r a t e s w i t h i n t h e n e t w o r k . // Note : This s u b r o u t i n e d o e s n o t have i n p u t argument . // I t u s e s t h e G l o b a l Net t o perform c a l c u l a t i o n s ! void c a l c q ( ) { i n t i , j , k , l ,m; double tmp0 , tmp1 , tmp2 , tmp3 , tmp4 , tmp5 ;  154  Appendix A. Code structure f o r ( i =1; i <IMAX−1; i ++) f o r ( j =1; j <JMAX−1; j ++) f o r ( k=1;k<KMAX−1;k++) f o r ( l =0; l <n u m b e r o f p h a s e s ; l ++){ tmp0=(Net [ i ] [ j ] [ k ] . f l u i d [ l ] . get pressure () − Net [ i ] [ j ] [ k − 1 ] . f l u i d [ l ] . g e t p r e s s u r e ( ) ) ∗ ( PI ) ∗ ( pow ( Net [ i ] [ j ] [ k ] . throat [ 0 ] . get diamater () , 4.0) ) / ( 128.0 ∗ Net [ i ] [ j ] [ k ] . f l u i d [ l ] . get viscosity () ∗ Net [ i ] [ j ] [ k ] . t h r o a t [ 0 ] . get length () ); Net [ i ] [ j ] [ k ] . t h r o a t [ 0 ] . s e t q ( tmp0 , l ) ; ... }; void copyNet ( p o r e NetSource [IMAX ] [ JMAX ] [KMAX] , p o r e NetTarget [IMAX ] [ JMAX ] [KMAX] ) // This a s s i g n s t h e hydrodynamic p r o p e r t i e s ( i . e . d e n s i t y , // v i s c o s i t y , c o n t a c t a n g l e , and s u r f a c e t e n s i o n t o t h e // n e t w o r k . ) //NOTE: a s s i g n i n g o f c o n t a c t a n g l e i s a random p r o c e d u r e . //The n e t w o r k i s c h a n g i n g a f t e r w a r d s ! void a s s i g n h y d r o d y n a m i c s ( ) { . . . } ; //TREATMENT t r e a t s t h e n e t w o r k t o d i s t r i b u t e // c o n t a c t a n g l e s t o t h e p o r e s . . . void t r e a t m e n t ( ) { . . . } ; // This r e a d s t h e n e t w o r k from t h e f i l e : f i l e n a m e void readNet ( char f i l e n a m e [ ] ) { . . . } ; //−−−−−−−−−−−−WRITING THE Network Diameter void w r i t e N e t D i a ( ) { . . . } ; void w r i t e N e t P r e s s u r e ( ) { . . . } ; void w r i t e P o r e S a t u r a t i o n ( ) { . . . } ; void w r i t e N e t S a t u r a t i o n ( ) { . . . } ;  Following declaration, the main body of the code starts in which the subroutines and the functions defined in the declaration and within the classes are used in the order explained in the Chapter 3, under the “Algorithm” section, to displace liquid water within the network until the condition of breakthrough is achieved.  155  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0073812/manifest

Comment

Related Items