UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Numerical simulations of gas-liquid two-phase flow in Polymer Electrolyte Membrane fuel cells Ding, Yulong 2012

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
[if-you-see-this-DO-NOT-CLICK]
ubc_2012_fall_ding_yulong.pdf [ 8.96MB ]
Metadata
JSON: 1.0059060.json
JSON-LD: 1.0059060+ld.json
RDF/XML (Pretty): 1.0059060.xml
RDF/JSON: 1.0059060+rdf.json
Turtle: 1.0059060+rdf-turtle.txt
N-Triples: 1.0059060+rdf-ntriples.txt
Original Record: 1.0059060 +original-record.json
Full Text
1.0059060.txt
Citation
1.0059060.ris

Full Text

Numerical Simulations of Gas-Liquid Two-Phase Flow in Polymer Electrolyte Membrane Fuel Cells by Yulong Ding B.A.Sc., Tsinghua University, 2005 M.A.Sc., Tsinghua University, 2008  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Chemical and Biological Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) July 2012 © Yulong Ding, 2012  Abstract Water management in PEM fuel cells has received extensive attention for its key role in fuel cell operation. Several water management issues have been identified that needed further investigation, i.e., droplet behaviour on the GDL surface, twophase flow patterns in gas flow channels, impact of two-phase flow on PEM fuel cell performance, impact of flow mal-distribution on PEM fuel cell performance, and mitigation of flow maldistribution. In this work, those issues were investigated based on simulations using computational fluid dynamics (CFD) method. Using the Volume of Fluid (VOF) two-phase flow model, droplet behaviour and two-phase flow patterns in mini-channels were identified consistently in both simulations and experimental visualizations. The microstructure of the GDL was found to play a significant role in the formation of local two-phase flow patterns, and the wettability of both GDL and channel wall materials greatly impacted on the two-phase flow patterns. A novel 1+3D two-phase flow and reaction model was developed to study the impact of two-phase flow on PEM fuel cell performances. The existence of twophase flow, especially the slug flow, in gas flow channels was found to be detrimental to the fuel cell performance and stability. Uneven liquid flow distribution into two parallel gas channels significantly reduces the fuel cell output voltage because of the induced severe non-uniform gas distribution, which should be avoided in the operation due to its negative effect on the fuel cell performance and durability. Finally, several maldistribution mitigation methods were tested in the simulation. It was found that utilizing narrow communication channels or adding gas inlet resistances could effectively reduce the gas flow maldistribution.  ii  Preface Section 1.3 of Chapter 1 is based on a review paper written by our research group, which was published: R. Anderson, L. Zhang, Y. Ding, M. Blanco, X. Bi, D.P. Wilkinson, A critical review of two-phase flow in gas flow channels of proton exchange membrane fuel cells, J. Power Sources 195 (2010) 4531-4553. The paper was under the supervision of Dr. Xiaotao Bi and Dr. David Wilkinson. I was responsible for writing the „CFD simulations of two-phase flow in PEM fuel cell flow channels‟ section and the simulation part of „Summary and outlook‟ section. Dr. Lifeng Zhang wrote the section „In situ experimental two-phase flow studies in PEM fuel cells‟, and PhD candidate Ryan Anderson wrote „Introduction‟, „Experimental visualization techniques‟, and „Ex situ experimental two-phase flow studies relevant to PEM fuel cells‟. Dr. Mauricio Blanco wrote „Water mitigation strategies‟. Chapter 3 is based on a simulation work that I conducted, analyzed and wrote, which was published under the supervision of Dr. Xiaotao Bi and Dr. David Wilkinson: Y. Ding, H.T. Bi, D.P. Wilkinson. Three-dimensional numerical simulation of water droplet emerging from a gas diffusion layer surface in microchannels. J. Power Sources 2010; 195: 7278-7288. Part of simulation results were also presented in a conference: Ding, Y., H.T. Bi, and D.P. Wilkinson. 3-D Numerical Simulation of Water Droplet Emerging from a Non-uniform Gas Diffusion Layer Surface. 20th International Symposium on Transport Phenomena (ISTP-20). 2009. Victoria, Canada. Chapter 4 is based on a simulation paper which I conducted under the supervision of Dr. Xiaotao Bi and Dr. David Wilkinson: Y. Ding, H.T. Bi, D.P. Wilkinson. Three-dimensional numerical simulation of water droplet emerging from a gas diffusion layer surface in micro-channels. J. Power Sources 2010; 195: 7278-7288. Part of the simulation results was also presented in a conference: Ding, Y., X.T. Bi, and D.P. Wilkinson, 3-D Numerical Simulation of Gas-Liquid iii  Flow in a Minichannel With a Non-Uniform GDL Surface. 8th International Conference on Nanochannels, Microchannels, and Minichannels, 2010, Montreal, Canada. Chapter 6 is based on the work that will be submitted for publication: Y. Ding, R. Anderson, L. Zhang, H.T. Bi, D.P. Wilkinson. Simulations of Two-phase Flow Distribution in Communicating Parallel Channels for a PEM Fuel Cell. This paper was supervised by Dr. Xiaotao Bi and Dr. David Wilkinson. I was responsible for conducting all the simulations and writing most of the paper. PhD candidate Ryan Anderson wrote the „introduction‟ and „Implications to PEM fuel cells‟. Dr. Lifeng Zhang provided feedback and suggestions on the simulation approach and data analysis. Chapter 7 is based on the simulation work I conducted under the supervision of Dr. Xiaotao Bi and Dr. David Wilkinson, and it has been submitted for publication: Y. Ding, H.T. Bi, D.P. Wilkinson. 3D Simulations of the Impact of Two-phase Flow on PEM Fuel Cell Performance. Chapter 9 is based on the experimental work that I collaborated with PhD candidate Ryan Anderson and summer student Gary Wang. Gary and I conducted the experiment. I was in charge of designing the channel used in section 9.1, establishing the experimental procedure, and analyzing experimental data. Ryan designed the two parallel channels used in sections 9.2 and 9.3 and provided suggestions for the experimental procedure.  iv  Table of contents Abstract ................................................................................................................. ii Preface ................................................................................................................ iii Table of contents .................................................................................................. v List of tables .......................................................................................................... x List of figures ....................................................................................................... xi Lists of symbols and abbreviations ................................................................... xviii Acknowledgements ............................................................................................ xxi Dedication ..........................................................................................................xxii 1. Introduction ....................................................................................................... 1 1.1. Background of fuel cell technology ............................................................. 1 1.2. Proton exchange membrane (PEM) fuel cells ............................................ 2 1.2.1. PEM fuel cell components .................................................................... 2 1.2.2. PEM fuel cell performance ................................................................... 6 1.3. Water management in PEM fuel cells ....................................................... 13 1.3.1. Experimental visualization techniques ................................................ 15 1.3.2. Gas-liquid two-phase flow models ...................................................... 17 1.3.3. Two-phase flow studies in PEM fuel cells........................................... 21 1.4. Research objectives and thesis layout...................................................... 36 2. Development of two-phase flow model of PEM fuel cells ................................ 40 2.1. 1+3D model .............................................................................................. 40 v  2.2. 1D MEA model ......................................................................................... 42 2.3. 3D VOF method ........................................................................................ 45 2.4. Water treatment ........................................................................................ 47 2.5. CFD algorithm and solution procedure ..................................................... 49 2.5.1. Precision improvement ....................................................................... 50 2.5.2. Mesh dependency .............................................................................. 51 3. Droplet behaviours in gas flow micro-channel................................................. 58 3.1. Computational domain and boundary conditions ...................................... 59 3.2. Results and discussion ............................................................................. 60 3.2.1. Droplet evolution ................................................................................ 60 3.2.2. Effects of water inlet structures .......................................................... 63 3.2.3. Effects of GDL surface contact angles ............................................... 68 3.2.4. Effects of channel wall surface contact angles ................................... 72 3.2.5. Effects of liquid flow rates ................................................................... 76 3.3. Conclusions .............................................................................................. 79 4. Two-phase flow patterns in gas flow channel ................................................. 81 4.1. Computational domain and boundary conditions ...................................... 82 4.2. Results and discussion ............................................................................. 84 4.2.1. Two-phase flow patterns in gas flow channels ................................... 84 4.2.2. Effects of liquid flow rates ................................................................... 86 4.2.3. Effects of GDL surface wettability ....................................................... 91 vi  4.2.4. Effects of channel wall surface wettability .......................................... 95 4.3. Conclusions .............................................................................................. 99 5. Impacts of two-phase flow maldistribution on flow behaviours in parallel gas flow channels .................................................................................................... 101 5.1. Computational domain and boundary conditions .................................... 102 5.2. Results and discussion ........................................................................... 103 5.2.1. Effects of severe flow maldistibution ................................................ 103 5.2.2. Effects of the degree of flow maldistribution ..................................... 107 5.2.3. Effects of GDL surface hydrophobicity on flow maldistribution ......... 112 5.2.4. Effects of channel wall hydrophobicity on flow maldistribution ......... 116 5.2.5. Effects of liquid inlet position on flow maldistribution ........................ 119 5.3. Conclusions ............................................................................................ 121 6. Communicating channels for flow mal-distribution mitigation ........................ 123 6.1. Computational domain and boundary conditions .................................... 124 6.2. Results and discussion ........................................................................... 126 6.2.1. Effects of communication channels on flow maldistribution .............. 126 6.2.2. Effects of communication channels on liquid slug removal............... 132 6.3. Implications of communication channels to PEM fuel cells ..................... 137 6.4. Conclusions ............................................................................................ 139 7. 3D simulations of the impact of two-phase flow on PEM fuel cell performance .......................................................................................................................... 141 7.1. Computational domain and parameters .................................................. 142 vii  7.2. Results and discussion ........................................................................... 144 7.2.1. Effects of stoichiometric flow ratios for single phase flow ................. 144 7.2.2. Effects of two-phase flow on the PEMFC performance .................... 144 7.2.3. Effects of severe flooding on cell performance ................................. 150 7.2.4. Effects of gas stoichiometric flow ratio on the PEMFC performance 152 7.2.5. Effects of GDL surface wettability on fuel cell performance.............. 157 7.2.6. Effects of channel wall wettability on fuel cell performance .............. 161 7.3. Conclusions ............................................................................................ 165 8. Effects of two-phase flow maldistribution on PEM fuel cell performance ...... 167 8.1. Computational domain and parameters .................................................. 167 8.2. Results and discussion ........................................................................... 168 8.2.1. Effects of flow maldistribution severity on fuel cell performance....... 168 8.2.2. Effects of liquid generation rate on flow maldistribution .................... 173 8.2.3. Effects of GDL surface wettability on flow maldistribution ................ 175 8.2.4. Effects of channel wall wettability on flow maldistribution ................. 177 8.2.5. Effects of cathode side gas stoichiometric flow ratio on flow maldistribution ............................................................................................ 179 8.2.6. Effects of inlet resistances on flow maldistribution ........................... 182 8.3. Conclusions ............................................................................................ 187 9. Experimental observation and evaluation ..................................................... 189 9.1. Droplet development in a transparent mini-channel................................ 189  viii  9.2. Two-phase flow maldistribution in parallel channels ............................... 193 9.3. Effects of communication channels ........................................................ 197 9.4. Conclusions ............................................................................................ 199 10. Conclusions and recommendations ............................................................ 201 10.1. Conclusions .......................................................................................... 201 10.1.1. Droplet behaviour ........................................................................... 201 10.1.2. Two-phase flow pattern development ............................................. 202 10.1.3. Impact of two-phase flow on PEM fuel cell performance ................ 202 10.1.4. Flow maldistribution and its impact on the cell performance .......... 202 10.1.5. Impact of material wettability .......................................................... 204 10.1.6. Maldistribution mitigation methods ................................................. 205 10.2. Suggestions for ideal gas flow channel design ..................................... 206 10.3. Recommendations for future work ........................................................ 206 10.3.1. Modelling development................................................................... 207 10.3.2. Experimental validation .................................................................. 207 10.3.3. Further simulation/experiment work ............................................... 208 References ....................................................................................................... 210  ix  List of tables Table 2.1: Assumptions of 1+3D model. ............................................................. 48 Table 7.1: Constants and operating parameters. .............................................. 143  x  List of figures Figure 1.1: Schematic of a PEM fuel cell .............................................................. 2 Figure 1.2: Flow field designs: (A) straight and parallel; (B) interdigitated; (C) serpentine. (Reprinted from Li et al. [2] with permission from Elsevier) ................ 3 Figure 1.3: SEM micrographs of: (a) carbon paper and (b) carbon cloth. (Reprinted from Wang et al. [4] with permission from Elsevier) ............................ 4 Figure 1.4: Schematic of catalyst layer. (Reprinted from Djilali and Sui [5] with permission from Taylor & Francis) ........................................................................ 5 Figure 1.5: Typical polarization curve of PEM fuel cells. (Reprint from Larminie et al. [10] with permission from John Wiley and Sons).............................................. 7 Figure 1.6: Effect of GDL hydrophobicity on droplet shape. ................................ 23 Figure 1.7: Schematic of contact angle hysteresis (A - advancing contact angle, R - receding contact angle). ................................................................................ 24 Figure 1.8: Droplet detachment height versus air velocity. Both hysteresis and non-hysteresis results are presented for comparison with experimental results. (Reprinted from Chen et al. [91] with permission from Elsevier) ......................... 27 Figure 1.9: Typical flow patterns in PEM fuel cell flow channels. (Reprinted from Hussaini et al. [98] with permission from Elsevier) .............................................. 29 Figure 1.10: Flow patterns in PEM fuel cells. ...................................................... 30 Figure 1.11: Ex-situ flow patterns in terms of superficial gas and liquid velocities. (Reprinted from Lu et al. [99] with permission from Elsevier) .............................. 32 Figure 1.12: Flow maldistribution in both ex-situ and in-situ PEM fuel cell experiment. (Reprinted from Kandlikar et al. [113] with permission from Elsevier) ............................................................................................................................ 33 Figure 2.1: Illustration of 1+3D cathode side of PEM fuel cell. ............................ 41 Figure 2.2: 1+3D model coupling method. .......................................................... 42 Figure 2.3: Three dimensional computational domain. ....................................... 50 Figure 2.4: Effects of solver precision on droplet behaviour. .............................. 51 Figure 2.5: Top view of meshes for base case, improved case X1, and refined cases X2 and X4 ................................................................................................. 52 Figure 2.6: Time evolution of the water coverage ratio on the two side walls ..... 53 xi  Figure 2.7: Top view of symmetrical meshes. ..................................................... 54 Figure 2.8: Time evolution of water coverage ratio on the two side walls. .......... 55 Figure 2.9: Side views of non-uniform meshes ................................................... 56 Figure 2.10: Comparison of FLUENT simulation results with literature results. (Modified from Chen et al. [91] with permission from Elsevier) ........................... 57 Figure 3.1: Three dimensional computational domain. ....................................... 59 Figure 3.2: Three stages of the emerging water droplet into a channel: (a) emergence and merging, (b) accumulation, and (c) detachment. ....................... 61 Figure 3.3: Time evolution of water distribution and pressure drop..................... 62 Figure 3.4: Different water injection inlet structures representing the GDL: (a) Uniform, (b) 1-pore, (c) 4-pores, (d) 16-pores (base case), and (e) 64-pores. .... 64 Figure 3.5: Effects of water inlet configurations on the two-phase flow patterns: (a) Uniform, (b) 1-pore, (c) 4-pore, and (d) 64-pore. ................................................ 65 Figure 3.6: Effects of water inlet configurations on the time-averaged water distribution. ......................................................................................................... 66 Figure 3.7: Effects of water inlet configurations on the time-averaged pressure drop..................................................................................................................... 67 Figure 3.8: Effects of GDL surface wettabilities on the two-phase flow patterns:(a) θ=0o, (b) θ=30o, (c) θ=45o, (d) θ=60o,(e) θ=90o, (f) θ=120o, and (g) θ=180o. ...... 70 Figure 3.9: Effects of GDL surface wettabilities on time-averaged water distribution. ......................................................................................................... 71 Figure 3.10: Effects of GDL surface wettabilities on time-averaged pressure drop. ............................................................................................................................ 72 Figure 3.11: Effects of channel wall wettabilities on two-phase flow patterns:(a) θ=0o, (b) θ=30o, (c) θ=60o, (d) θ=90o,(e) θ=120o, (f) θ=140o, and (g) θ=180o. ... 74 Figure 3.12: Effects of channel wall wettabilities on time-averaged water distribution. ......................................................................................................... 75 Figure 3.13: Effects of channel wall wettabilities on time-averaged pressure drop. ............................................................................................................................ 76 Figure 3.14: Effects of liquid flow rates on two-phase flow patterns: (a) 1/10th, (b) 1/100th. ................................................................................................................ 77 Figure 3.15: Effects of liquid flow rates on time-averaged water distribution. ..... 78 Figure 3.16: Effects of liquid flow rates on time-averaged pressure drop. .......... 78 xii  Figure 4.1: Three-dimensional computational domain: (a) overview (b) simplified GDL microstructure (c) Local view of GDL microstructure with meshes. ............ 83 Figure 4.2: Two-phase flow patterns in the channel............................................ 85 Figure 4.3: Time-averaged water coverage ratio on different surfaces along the channel. .............................................................................................................. 86 Figure 4.4: Effects of liquid flow rate on the two-phase flow patterns in the channel: (a) x10 case, (b) x100 case, and (c) x1000 case. ................................. 88 Figure 4.5: Time-averaged water coverage ratio on different surfaces along the channel: (a) x10 case, (b) x100 case, and (c) x1000 case. ................................. 89 Figure 4.6: Effects liquid flow rate on the pressure drop and water volume fraction in the channel. ........................................................................................ 91 Figure 4.7: Effects of contact angle of GDL surface on the two-phase flow patterns in the channel :(a) θ=45o, (b) θ=60o,(c) θ=90o, (d) θ=120o, (e) θ=140o. 93 Figure 4.8: Effects of GDL surface contact angle on the water distribution in the channel. .............................................................................................................. 94 Figure 4.9: Effects of GDL surface contact angle on the pressure drop and water volume fraction in the channel. ........................................................................... 95 Figure 4.10: Effects of contact angle of channel side and top wall surfaces on the two-phase flow patterns in the channel: (a) θ=45o, (b) θ=60o, (c) θ=90o, (d) θ=120o, (e) θ=140o. ............................................................................................. 97 Figure 4.11: Effects of channel wall contact angle on the water distribution in the channel. .............................................................................................................. 98 Figure 4.12: Effects of channel wall contact angle on the pressure drop and water volume fraction in the channel. ........................................................................... 99 Figure 5.1: Three-dimensional computation domain of two-parallel channel. ... 102 Figure 5.2: Effects of maldistribution on two-phase flow patterns in parallel channels: (a) w/o maldistribution, (b) with maldistribution ................................. 104 Figure 5.3: Effects of maldistribution on averaged gas velocity in each branch. .......................................................................................................................... 105 Figure 5.4: Effects of maldistribution on total pressure drop and liquid volume fractions in each branch. ................................................................................... 106 Figure 5.5: Effects of liquid maldistribution on wall liquid coverage in each branch. .......................................................................................................................... 107  xiii  Figure 5.6: Effects of the degree of flow maldistribution on two-phase flow patterns. ............................................................................................................ 108 Figure 5.7: Effects of the degree of flow maldistribution on (a) gas velocity (b) gas flow rate ratio. ................................................................................................... 109 Figure 5.8: Effects of the degree of flow maldistribution on GDL surface liquid coverage ratio. .................................................................................................. 110 Figure 5.9: Effects of the degree of flow maldistribution on pressure drop. ...... 111 Figure 5.10: Effects of the degree of flow maldistribution on liquid volume fraction. .......................................................................................................................... 112 Figure 5.11: Effects of GDL surface hydrophobicity on the two-phase flow patterns :(a) θ=45o, (b) θ=60o, (c) θ=90o, (d) θ=120o. ....................................... 113 Figure 5.12: Effects of GDL surface hydrophobicity on (a) pressure drop (b) liquid volume fraction.................................................................................................. 114 Figure 5.13: Effects of GDL surface hydrophobicity on GDL liquid coverage ratio. .......................................................................................................................... 115 Figure 5.14: Effects of GDL surface hydrophobicity on flow maldistribution. .... 115 Figure 5.15: Effects of channel wall hydrophobicity on the two-phase flow patterns :(a) θ=60o, (b) θ=90o, (c) θ=120o, (d) θ=140o. ..................................... 116 Figure 5.16: Effects of channel wall hydrophobicity on (a) pressure drop (b) liquid volume fraction.................................................................................................. 117 Figure 5.17: Effects of channel wall hydrophobicity on GDL liquid coverage ratio. .......................................................................................................................... 118 Figure 5.18: Effects of channel wall hydrophobicity on flow maldistribution. ..... 119 Figure 5.19: Comparison of water inlet positions of two cases. ........................ 120 Figure 5.20: Simulated two-phase flow pattern for the mid-inlet case. .............. 120 Figure 5.21: Effects of liquid inlet position on the average gas velocity in each branch. .............................................................................................................. 121 Figure 6.1: Three-dimensional computational domain (a) illustration of parallel communicating channels (b) base case ............................................................ 125 Figure 6.2: Snapshots two-phase flow patterns in each case. .......................... 127 Figure 6.3: Time-averaged water volume fraction in each branch. ................... 128 Figure 6.4: Time-averaged total GDL water coverage ratio. ............................. 129 xiv  Figure 6.5: Time-averaged superficial velocity along the left branch. ............... 131 Figure 6.6: Time-averaged total pressure drop. ................................................ 132 Figure 6.7: Snapshots of initial and final water distribution in channels with communications. ............................................................................................... 134 Figure 6.8: Temporal variation of water volume fraction in each parallel channel branch. .............................................................................................................. 135 Figure 6.9: Temporal variation of total GDL water coverage ratio. ................... 136 Figure 6.10: Final state of superficial gas velocity along the channel length. .... 137 Figure 7.1: Illustration of 1+3D cathode side of a PEM fuel cell. ....................... 142 Figure 7.2: Effects of gas stoichiometric flow ratios for single phase flow. ........ 144 Figure 7.3: Two-phase flow pattern in the cathode side channel. ..................... 145 Figure 7.4: Effects of liquid water generation rates on two-phase flow patterns. (a) 10 times, (b) 100 times, (c) 1000 times. ............................................................ 147 Figure 7.5: Effects of liquid water generation rates on the cell voltage and pressure drop. ................................................................................................... 148 Figure 7.6: Effects of liquid water generation rates on liquid wall coverage and volume fraction.................................................................................................. 149 Figure 7.7: Fuel cell polarization curves under severe flooding conditions. ...... 150 Figure 7.8: Effects of operating current densities on overall pressure drop. ..... 151 Figure 7.9: Effects of operating current densities on water wall coverage and liquid holdup. ..................................................................................................... 152 Figure 7.10: Two-phase flow patterns at different gas stoichiometric flow ratios. (a) stoich.=2 (b) stoich.=3 (c) stoich.=4 (d) stoich.=5 (e) stoich.=6. .................. 154 Figure 7.11: Effects of gas stoichiometric flow ratio on the fuel cell voltage output and pressure drop. ............................................................................................ 155 Figure 7.12: Effects of gas stoichiometric flow ratio on water wall coverage ratio and volume fraction........................................................................................... 156 Figure 7.13: Effects of gas stoichiometric flow ratio on polarization curves. ..... 157 Figure 7.14: Effects of GDL surface contact angle on the two-phase flow patterns in the channel at 800 mA/cm2: (a) θ=45o, (b) θ=60o, (c) θ=120o, (d) θ=140o..... 158 Figure 7.15: Effects of GDL surface contact angle on water wall coverage and holdup at 800 mA/cm2. ...................................................................................... 159 xv  Figure 7.16: Effects of GDL surface contact angle on the fuel cell pressure drop at 800 mA/cm2. ................................................................................................. 160 Figure 7.17: Effects of GDL surface contact angle on the fuel cell voltage output at three different current densities..................................................................... 161 Figure 7.18: Effects of channel wall contact angle on the two-phase flow patterns in the channel at 800 mA/cm2: (a) θ=45o, (b) θ=60o, (c) θ=120o, (d) θ=140o..... 162 Figure 7.19: Effects of channel wall contact angle on water wall coverage and holdup at 800 mA/cm2. ...................................................................................... 163 Figure 7.20: Effects of channel wall contact angle on the fuel cell pressure drop at 800 mA/cm2. ................................................................................................. 164 Figure 7.21: Effects of channel wall contact angle on the fuel cell voltage output at three different current densities..................................................................... 165 Figure 8.1: Three-dimensional computation domain of two parallel channels. .. 168 Figure 8.2: Effects of flow maldistribution severity on the two-phase flow patterns in the channel under three different current densities: (a) 100 mA/cm2, (b) 400 mA/cm2, (c) 800 mA/cm2. .................................................................................. 169 Figure 8.3: Effects of flow maldistribution severity on the fuel cell polarization curve. ................................................................................................................ 170 Figure 8.4: Effects of flow maldistribution severity on the total pressure drop. . 171 Figure 8.5: Effects of flow maldistribution severity on the GDL water coverage ratio. .................................................................................................................. 172 Figure 8.6: Effects of flow maldistribution severity on the gas velocity ratio ...... 173 Figure 8.7: Effects of liquid generation rate on the fuel cell output voltage at a current density of 800 mA/cm2. ......................................................................... 174 Figure 8.8: Effects of liquid generation rate on the gas velocity ratio. ............... 175 Figure 8.9: Effects of GDL surface contact angle on the fuel cell output voltage and the GDL water coverage ratio. ................................................................... 176 Figure 8.10: Effects of GDL surface contact angle on the total pressure drop and the gas velocity ratio. ........................................................................................ 177 Figure 8.11: Effects of channel wall contact angle on the fuel cell output voltage and the GDL water coverage ratio. ................................................................... 178 Figure 8.12: Effects of channel wall contact angle on the total pressure drop and the gas velocity ratio. ........................................................................................ 179 xvi  Figure 8.13: Effects of gas stoichiometric flow ratio on the fuel cell output voltage. .......................................................................................................................... 180 Figure 8.14: Effects of gas stoichiometric flow ratio on the GDL water coverage ratio. .................................................................................................................. 181 Figure 8.15: Effects of gas stoichiometric flow ratio on the gas velocity ratio. .. 182 Figure 8.16: Effects of inlet resistance on the fuel cell output voltage............... 184 Figure 8.17: Effects of gas inlet resistance on the total pressure drop.............. 185 Figure 8.18: Effects of gas inlet resistance on the GDL water coverage ratio. .. 186 Figure 8.19: Effects of inlet resistance on the gas velocity ratio. ...................... 187 Figure 9.1: Picture and schematic of the transparent mini-channel. ................. 190 Figure 9.2: Side views of liquid droplet development. ....................................... 191 Figure 9.3: Top views of two-phase flow patterns. ............................................ 193 Figure 9.4: Picture and schematic of two parallel channels. ............................. 194 Figure 9.5: Effects of the degree of flow maldistribution on gas flow rate ratio. 195 Figure 9.6: Effects of the degree of flow maldistribution on normalized pressure drop................................................................................................................... 197 Figure 9.7: Effects of communication channels on gas flow rate ratio between parallel channels. .............................................................................................. 198 Figure 9.8: Effects of communication channels on normalized pressure drops.199  xvii  Lists of symbols and abbreviations Bo  Bond number, dimensionless  C  concentration, mol/m3  Ca  Capillary number, dimensionless  D  diffusion coefficient, m2/s  E  electric potential, V  E0  open circuit potential, V  F  Faraday‟s constant, 96,485.339 C/mol  Fvol  momentum source term, kg/(m2∙s2)  g  gravity, m/s2  G  Gibbs free energy, kJ  H  enthalpy, kJ  H0  Henry‟s constant, Pa∙m3/mol  I  current, A  i  current density, A/cm2  i0  exchange current density, A/cm2  J  molar flux, mol/(m2∙s)  L  characteristic length, m  n  number of electrons  n  unit vector normal, dimensionless  N  mole flow rate. mol/s  P(p)  pressure, Pa  xviii  R  gas constant, 8.314 J/(K ∙ mol)  Re  Reynolds number, dimensionless  S  entropy, kJ/K  SA  area, m2  Sm  mass source term, kg/(m3∙s)  Su  Suratman number, dimensionless  t  unit tangential normal, dimensionless  T  temperature, K  U  local potential, V  Vout  voltage, V  V  velocity, m/s  W  work, W  Y  species mass fraction kg/m3  Greek letters α  transfer coefficient, dimensionless  αhalf  half angle of channel corner  β  resistance coefficient, dimensionless  γ  surface tension, N/m  δ  mass transfer coefficient, m/s  Δ  change in value, dimensionless  ε  volume fraction, dimensionless  η  efficiency, dimensionless  ηc  cathode overpotential, V xix  θ  contact angle, o  κliquid  surface curvature of the liquid droplet, dimensionless  λ  gas stoichiometric ratio, dimensionless  μ  dynamic viscosity, Pa∙s  ρ  density, g/cm3  ζ  conductivity, S/m2  ζsurf  surface tension coefficient, N/m  Acronyms CCD  charge-coupled device  CFD  computational fluid dynamics  FC  fuel cell  FVM  finite volume method  GDL  gas diffusion layer  HOR  hydrogen oxidation reaction  LBM  Lattice Boltzmann method  MEA  membrane electrode assembly  PEM  proton exchange membrane  PISO  pressure-based segregated solver  PTFE  polytetrafluoroethylene  VOF  volume of fluid method  xx  Acknowledgements I would like to express my sincere gratitude to my supervisors, Dr. Xiaotao Bi and Dr. David Wilkinson for their support, professional supervision and great encouragement throughout my PhD program. Their diligence, self-discipline and enthusiasm always inspire me and set me a good example to follow for the rest of my life. I would also like to extend my thanks to my committee members, Dr. James Feng and Dr. Kendal Bushe, as well as Dr. Jürgen Stumper and Dr. Haijiang Wang for their constructive criticism and valuable advice on my research. Thanks to my high quality research team: Dr. Lifeng Zhang, Ryan Anderson, Cyrus Yau. Thank you for all the help you have given to me in the meetings, in my research and in my life. Many thanks go to Gary Wang for assisting me in conducting all the experiments. I‟d also like to appreciate all faculty, staff, and students in our department who all together create a friendly, supportive, and cooperative environment. I really enjoy the time with them. Special thanks go to my close friends, Cai Long, Jun Xu, Wendi Guo, Jianghong Peng, Wei Wang, Xiao He who have made my life in Vancouver exciting, fascinating, and memorable. At last, I would like to thank my mom and my dad for their love, unconditional support, and tireless dedication throughout my life.  xxi  Dedication  Dedicated to my parents for their love, encouragement and endless support  xxii  1. Introduction1 1.1. Background of fuel cell technology Energy supply and use are of fundamental importance to our modern society. However, the three major energy sources nowadays, i.e., petroleum, nature gas and coal are all fossil fuels, which are non-renewable, and cause great impact on environment, such as global warming, air and water pollution. As society further develops, the demand for these energy sources will grow rapidly, which may cause more severe environmental problems, and as the fossil fuels deplete gradually, searching for alternative energy sources becomes more and more urgent. Hydrogen, which is the most abundant chemical element in the universe, is widely regarded as the next generation clean energy source. However, how to convert the hydrogen into power efficiently is still a challenge work. Fuel cells, which directly convert the chemical energy into electricity, show great potential as a high-efficiency energy-conversion device. In recent years, fuel cell research has been receiving much attention, both academic and industries. A number of applications have been developed and demonstrated in the markets, such as automobiles, distributed power generation, back-up power, portable power, space, and airplanes [1]. Several common types of fuel cells, characterized by the different electrolyte used, have been developed, such as polymer-electrolyte membrane (PEM) fuel cell, alkaline fuel cell, phosphoric acid fuel cell, molten carbonate fuel cell and solid oxide fuel cell. Among those types of fuel cells, PEM fuel cells have been  1  Sections of this chapter were published. Anderson R, Zhang L, Ding Y, Blanco M, Bi X, Wilkinson DP. A critical review of two-phase flow in gas flow channels of proton exchange membrane fuel cells. Journal of Power Sources 2010;195:453153.  1  receiving the most attention for their low operating temperatures (60~90℃), simpler design and lower cost. However, the current PEMFCs are still more expensive than conventional fuel combustion engines, thus, in-depth research should be carried out to reduce the cost and increase the power density of PEM fuel cells.  1.2. Proton exchange membrane (PEM) fuel cells 1.2.1. PEM fuel cell components Figure 1.1 shows a typical unit of PEM fuel cell, which consists of three parts the anode side, the cathode side and the electrolyte. Both the anode and cathode sides include gas flow channels, gas diffusion layers and catalyst layers, and the electrolyte is a polymer-electrolyte membrane.  Figure 1.1: Schematic of a PEM fuel cell  Hydrogen flows into the anode gas flow channel, diffuses through the gas diffusion layer, and then reacts on the catalyst layer where hydrogen ions and electrons are generated. The hydrogen ions diffuse through the polymer electrolyte membrane, and the electrons flow through anode gas diffusion layer  2  to the current collector and into the electric load. Electrons enter the cathode side through the cathode current collectors and gas diffusion layer. At the cathode catalyst layer, the electrons, the hydrogen ions and the oxygen react to form water. 1.2.1.1. Gas flow channels In PEM fuel cell applications, bipolar plates are generally used to separate different cell units and the gas flow channels are carved into the bipolar plates. Gas flow channels provide a pathway of reactants. Meanwhile excess liquid water will be expelled into the flow channels. Because gas channel flooding will have a damaging effect on the performance of PEM fuel cell, the design of gas flow channels has been given special considerations in the past. It was also found that the appropriate design of flow channels was the most effective strategy in tackling water flooding issues [2]. Three main types of flow-field developed to date are the straight and parallel flow field, the interdigitated flow field and the serpentine flow field (See Figure 1.2).  Figure 1.2: Flow field designs: (A) straight and parallel; (B) interdigitated; (C) serpentine. (Reprinted from Li et al. [2] with permission from Elsevier)  1.2.1.2. Gas diffusion layer (GDL) The gas diffusion layers, as shown in Figure 1.1, connect catalyst layers and gas channels and provide a pathway for electrons, gases, and liquid water to move to or from the catalyst layers [3]. Two materials are typically employed as gas diffusion layers in PEM fuel cells; carbon cloth and carbon paper (See Figure 1.3).  3  Both materials are porous media and fabricated from carbon fibers. In order to balance membrane hydration and flooding, the gas diffusion layers are typically treated with a PTFE (Teflon) during its manufacturing. PTFE, as a hydrophobic material, can help prevent liquid water flooding but allow water vapour to be transported to the membrane.  Figure 1.3: SEM micrographs of: (a) carbon paper and (b) carbon cloth. (Reprinted from Wang et al. [4] with permission from Elsevier)  1.2.1.3. Catalyst layer Due to the slow electrochemical reactions especially at the cathode side, a catalyst layer is essential to enhance the electrochemical reaction rates. The reactions can only occur at the three-phase boundary of catalyst surface where all the necessary reactants, namely gases, electrons and protons meet. Thus catalyst layers must be installed between the membrane and gas diffusion layers where electrons can travel through the electrically conductive fibers of gas diffusion layers and catalyst itself, reactant gases can travel through the voids of gas diffusion layers and protons can migrate through the membrane. The most commonly used catalyst in PEM fuel cells for both oxygen reduction and hydrogen oxidation reactions is platinum [1]. Typically the platinum is supported on carbon clump forming nano-porous agglomerate structures. The structure of catalyst is probably of a “spaghetti and meatball” structure, where the  4  carbon agglomerates are connected to each other and covered by thin tendrils of membrane [3]. (See Figure 1.4)  Figure 1.4: Schematic of catalyst layer. (Reprinted from Djilali and Sui [5] with permission from Taylor & Francis)  Considering both the transport of protons and the permeation of reactant gases, the catalyst layers should be made reasonably thin to reduce the potential losses. At the same time, Pt particles should be dispersed as finely as possible to maximize the active surface area of catalyst. The catalyst layer in PEM fuel cells is typically 10-30  m thick and can have platinum loading as low as 0.2 mg/cm2, with the diameter of Pt particles being 4 nm or smaller [6]. Besides, the Nafion in catalyst layer is hydrophilic and will absorb and retain liquid water. Thus it is important to remove excess liquid water to prevent its blocking of reaction sites. 1.2.1.4. Membrane The membrane in PEM fuel cells acts as the electrolyte to transport protons and prevents the cross-flow of reactant gases and electrons. The membrane used most often in the past was a melt-extruded membrane manufactured by DuPont and sold under the label Nafion No.117 [7]. This type of membrane consists of a  5  hydrophobic fluoro-carbon backbone attached with hydrophilic sulfonic acid (SO3-) groups. The proton conduction takes place via the sulfonic acid groups with the aid of water, because of their strong attraction to protons. Thus, the protonic conductivity is highly dependent on the fixed-charge concentration (the acidic groups) and the water content in the membrane. It is thus crucial to keep the membrane fully humidified at all times. Another factor that has a great impact on protonic conductivity is the thickness of the membrane. A thinner one reduces the ohmic losses of the membrane, but increases the reactants crossover, which causes the short circuiting of fuel cell. The typical thickness of a membrane is in the range of 5~200  m , and a fully humidified membrane‟s conductivity is usually around 0.1 S/cm [8]. Operating temperature is another important parameter to maintain high performance of membrane. For Nafion, the maximum operating temperature is about 125 oC, above which the chemical and thermal stability – as well as resistance to strong oxidizing or reducing agents – begins to decrease [9]. 1.2.2. PEM fuel cell performance The PEM fuel cell as one kind of energy conversion devices is more efficient than conventional combustors in theory. However, the irreversible process in the PEM fuel cell greatly reduces its efficiency. Thus, it is essential to understand the limitation and influential factors of its performance. The performance of a fuel cell can be captured with the polarization curve, which shows the voltage output for a given current load. A typical polarization curve for a PEM fuel cell is shown in Figure 1.5.  6  Figure 1.5: Typical polarization curve of PEM fuel cells. (Reprint from Larminie et al. [10] with permission from John Wiley and Sons)  An ideal fuel cell would maintain a constant voltage determined by thermodynamics (indicated by the dotted line). Unfortunately, the real voltage of a fuel cell is much lower than thermodynamically predicted one due to its irreversible losses, especially at higher current load. In a fuel cell, these irreversible losses are generally caused by different factors:  Kinetics of the electrochemical reactions  Crossover of reactants and internal currents  Electrical and ionic resistance  Mass transport limitation 1.2.2.1. Thermodynamics The basic electrochemical reactions in a PEM fuel cell are:  7  At the anode:  H2  2H+ +2e-  (1.1)  1 O2 +2H + +2e-  H 2O 2  (1.2)  1 O2 +H 2  H 2O 2  (1.3)  At the cathode:  Overall:  The maximum work obtainable from the reaction corresponds to the Gibbs free energy change of the reaction. According to basic thermodynamics,  ΔG=ΔH-TΔS  (1.4)  It means that not all the energy input can be converted into electricity, and some irreversible losses in energy conversion due to the entropy loss. At an environmental temperature of 25  o  C and ambient pressure, 237.34 KJ/mol  energy can be converted into electrical energy and 48.68 KJ/mol is converted into heat [1]. Meanwhile, the electrical work is equal to:  We =nFE  (1.5)  Where: n = number of electrons per molecule of H2 (= 2 electrons per molecule) F = 96,485 Coulombs/electron-mol (Faraday‟s constant) E = electric potential 8  As mentioned previously,  We,max =-ΔG  (1.6)  Thus, the theoretical potential of a fuel cell is then:  E max =  -ΔG nF  (1.7)  According to equation(1.7), at 25oC and ambient pressure, the theoretical potential of hydrogen/oxygen fuel cell is:  E0,max =  -ΔG 237340 J/mol = =1.229 V nF 2  96485 A/mol  (1.8)  Thus, the theoretical efficiency is (assuming liquid water formation):  ηmax =ΔG/ΔH=237.34/286.02=83%  (1.9)  At a different temperature, the theoretical potential can still be calculated by equation (1.4) and (1.7),  E max =  -ΔG  ΔH TΔS  =-   nF  nF nF   (1.10)  Both ΔH and ΔS are functions of temperature and can be calculated by the integration of heat capacity. However when the temperature is low (below 100 o  C), the enthalpy and entropy can be considered constant [1].  At different pressures or concentrations of gas reactants, the theoretical potential can be calculated by the Nernst equation: 0.5 RT  PH2 PO2  E max =E 0,max + ln   nF  PH2O   (1.11)  where P indicates the partial pressure of the reactants.  9  Although the theoretical efficiency predicted by equation (1.9) is much higher than that in conventional heat engines, it requires “equilibrium” state during the whole process to achieve that efficiency. It also means that infinite time is needed, which is impossible in a real fuel cell. Thus, some irreversible processes are inevitable, which will decrease a fuel cell‟s efficiency. 1.2.2.2. Electrochemical reactions For any electrochemical reaction, the reaction rate is finite due to the activation energy barrier that the charge must overcome. In reality, even the simple hydrogen oxidation reaction (HOR) also consists of a series of basic steps, and the overall reaction rate is limited by the slowest step which is always the electrons transfer between reactants and the electrode surface. On the other hand, electrons may transfer towards both directions, thus, an electrochemical reaction involves both oxidation and reduction of the species:  Ox+ne-  Red  (1.12)  The measured current is the net current, namely, the difference between forward and backward current. Generally, the net current density of any electrochemical reaction can be calculated by the Butler-Volmer equation:     F ( E  Er )   Ox F ( E  Er )   i  i0 exp  Rd  exp    RT RT       (1.13)  where: i0 = exchange current density, equal to the forward (or reverse) current density when the net current density is equal to zero (at equilibrium), α = transfer coefficient, depends on the symmetry of the activation barrier, F = Faraday‟s constant,  10  E = potential, Er = equilibrium potential. The exchange current density depends on the chemical reactions (concentration dependent), electrode catalyst loading, catalyst specific surface area, and also temperature and pressure. It is also a measurement of electrochemical reaction activity. In a PEM fuel cell, the exchange current density at the anode is much higher (usually several orders of magnitudes) than at the cathode. The equilibrium potential at the PEM fuel cell anode is 0 V by definition and at the cathode it is 1.229 V (at 25 0C and atmospheric pressure). 1.2.2.3. Voltage losses As mentioned above, voltage losses are caused by different factors, and some of factors will become predominant when current density is within a certain range. 1.2.2.3.1. Internal currents and crossover losses When the electrical circuit of a fuel cell is not closed, no current will be generated, and the potential of the fuel cell is called the open circuit potential. In principle, this potential should equal the theoretical potential. However, it is always significantly lower. This voltage loss is due to internal currents and reactant gases crossover. Although the membrane of fuel cell is effectively impermeable to reactant gases and electrons, some small amount of hydrogen will diffuse through the membrane and react with oxygen on the cathode side directly, or some electrons pass through the membrane resulting in external circuit (shortcircuiting). In addition, oxygen may also diffuse through the membrane. However, this permeation rate is much lower than that of hydrogen.  11  Reactant gases crossover is a function of membrane permeability, membrane thickness and gases concentration gradient. Thus, when electrical circuit is open or current density is very small, this kind of voltage loss will be predominant. However, as the current density increases, the gas concentration in the catalyst layer decreases. As a result, the gas crossover loss can be neglected. 1.2.2.3.2. Activation losses As shown in equation(1.13), some voltage difference from equilibrium potential, which is called overpotential, is present which is essential as a driving force to keep the electrochemical reaction going. This overpotential is also called activation loss, which indicates the voltage loss due to the slowed electrochemical reactions. Obviously, higher exchange current density results in lower activation loss. Thus, oxygen reduction reaction (ORR) in PEM fuel cell requires much higher overpotential than HOR and has a higher activation loss. According to the Butler-Volmer equation(1.13), at the cathode side, the overpotential is negative, which makes the second term negligible:    F ( Ec  Er ,c )  ic  i0,c exp  c  RT    (1.14)  Thus,  Vact ,c  Er ,c  Ec   RT  i  ln    c F  i0,c   (1.15)  According to equation(1.15), as the current density increases, activation losses first increase significantly at low current densities, and then remain constant. This tendency is the same at the anode side. Thus, this kind of voltage loss is always predominant at low current densities.  12  1.2.2.3.3. Ohmic losses Ohmic losses are due to the resistance when ions pass through the membrane and electrons pass through the electrically conductive parts of PEM fuel cell. It can be expressed by Ohm‟s law:  Vohm  iR  (1.16)  R indicates the total resistance of the fuel cell. Typically, electronic resistance can be neglected because of much higher conductivity of materials. Because ohmic losses are linear to the current density, this kind of voltage loss becomes predominant at high current densities. 1.2.2.3.4. Concentration losses Concentration losses occur when the gas reactants are consumed faster than supplied. Gas consumption depends on current density - higher current density results in lower gases concentration at the catalyst surface. Thus, this kind of voltage losses becomes predominant at very high current densities. On the other hand, gas supply depends on the diffusion rates, which are also functions of gas concentration gradients. Thus, the maximum current density will be reached when the gas concentration at the catalyst surface approaches zero. This current density is also called limiting current density. However, due to the non-uniform structure of porous media, this limiting current density is seldom reached in a real fuel cell. Some areas would reach “equilibrium” while other areas still not.  1.3. Water management in PEM fuel cells In the operation of PEM fuel cells, water management is a key issue to improve the fuel cell performance. Too little water will cause membrane dehydration, which limits proton conductivity, and too much water can flood the fuel cell, 13  causing less reactant to reach active catalyst sites and consequently decreasing the cell performance. Therefore, the objective of water management is to maintain a proper water balance in a fuel cell. Usually the membrane dehydration issue can be addressed by using humidified reactant gas streams. However, on the other side, flooding issue is almost unavoidable, especially when the fuel cell is operating at high current density with humidified reactant gases. To address the flooding issue, several aspects must be fundamentally understood, such as, Liquid water generation at the catalyst layer. Water transport inside the membrane. Water transport inside the GDL. Gas-liquid two-phase flow in the flow field channel. A recent review [2] detailed issues associated with water management, which described the role of each layer of the PEM fuel cell „sandwich‟ and how each area is prone to flooding. Trabold [11] noted the importance of two-phase flow research in PEM fuel cells, and explained that the gas-liquid flow within the flow channels is complex and requires an understanding of electrochemistry, heat and mass transfer, and fluid mechanics. Meanwhile, the first three aspects, i.e., water transport inside the MEA (membrane electrode assembly) which happens at microscale or mesoscale, are beyond the scope of the research. Therefore, in the following literature review, much attention is placed on gas-liquid two-phase flow issues in the flow channel.  14  1.3.1. Experimental visualization techniques 1.3.1.1. Optical visualization Directly observing a transparent fuel cell is probably the most convenient technique to study the two-phase flow in the flow field. Charge-Coupled Device (CCD) camera is often used to capture two-phase flow patterns through a transparent window [12-19]. One concern about this technique is that the materials used for the transparent window may not be able to mimic the commercially used graphite flow field plates. It is found that the hydrophobicity of material has great impact on the two-phase flow. Therefore the flow patterns observed from transparent fuel cells may differ from those in a real fuel cell. Because most of the transparent materials are not electrically conductive, they cannot be used to make contact with GDL as an electron collector, which means only a top view, not a side view can be provided from a transparent fuel cell. Even in ex-situ experiment, while dealing with multiple channels, it is still impossible to observe the two-phase flow in the crossthrough GDL direction. Another important material property is the surface roughness, which influences the total pressure drop [20] or water condensation [21], and in turn may alter the two-phase flow patterns. Another drawback of this technique is that it is not able to provide quantitative information, i.e., volume of liquid water, water coverage ratio on each wall, etc. Some quantities such as droplet velocity can be measured from post-image processing, but the reflective GDL background usually complicates the processing [17]. It is useful to correlate the images with a pressure drop signal or voltage signal, but it still provides only qualitative information about the cell [22]. Despite these drawbacks, optical visualization of fuel cells provides a method to validate existing models and to understand the impacts of key variables, which is particularly important in incorporating two-phase flow into existing models [23].  15  1.3.1.2. Other visualization methods Other visualization methods include neutron radiography [24-31], magnetic resonance imaging (MRI) [32-37] and the X-ray method [38-43]. They are all noninvasive and can quantify liquid water accumulation in operating PEM fuel cells. Neutron radiography uses neutron beams, which can detect high hydrogen content organic materials or water, to obtain 2-D images of liquid water. This technique allows the user to gain greater quantitative information. One drawback of neutron radiography is that it is difficult to distinguish between liquid water and vapor. Also the 2-D images can only provide through-plane or in-plane water distribution [28]. Although reconstructing three-dimensional water distribution by using more sets of neutron beams is possible [30], it has only ever been applied to an inactive PEM fuel cell. Furthermore, neutron imaging typically provides a spatial resolution of only 0.1 mm by 0.1 mm, which is insufficient to analyze the through-plane evolution of water transport [44]. Magnetic resonance imaging (MRI) is widely available, inherently threedimensional, and capable of visualizing water in opaque structures [45]. Unlike optical visualization, this method can detect the liquid water distribution under the gas channel and landing areas. However, limited temporal resolution, limited inplane spatial resolution [33] and limited size of the magnet-core for fuel cell housing restricts its application [44]. The X-ray image technique can give high temporal and spatial resolution. The use of synchrotron radiation makes it capable of reaching higher spatial resolutions coupled with high signal to noise ratios [38]. However, this method is relatively unstable because the X-ray beam is easy to scatter and is absorbed by electrons, a feature which makes X-ray method less popular than the neutron method [46].  16  Recent reviews provides more details and in-depth discussion on water visualization techniques by these methods and others [44-48]. 1.3.2. Gas-liquid two-phase flow models In addition to experimental efforts, many attempts have been made to model and simulate two-phase flow and transport in PEM fuel cells. In particular, computational fluid dynamics (CFD) is considered to be a very powerful tool in fuel cell design and operation optimization. In the literature, the earliest PEM fuel cell models date back to the early 1990s by Springer et al. [49] and Bernardi and Verbrugge [50]. Increasingly complex and detailed models have been developed since then, from one-dimensional, single-phase flow, isothermal, steady state and single layer models to three-dimensional, two-phase flow, non-isothermal, transient and multiple layer models. However, due to the various complicated phenomena in PEM fuel cells, the modeling and simulation of PEM fuel cells still remains a challenge. Complications include two-phase flow, electrochemical reaction, charge transport, diffusion in porous media, and coupling different length scales (such as the nanometer components of catalysts, the micrometer scale heterogeneous pores in the gas diffusion layers (GDLs), and the millimeter scale of the flow field channels). In recent years, several reviews have been published about fuel cell models. Weber and Newman [51] presented various types of transport and corresponding models in each fuel cell layer. However, their review mainly focused on onedimensional models and two-phase flow in gas channels was not taken into account. Wang [52] summarized fundamental models for PEM fuel cell engineering but limited the review to computational fluid dynamic (CFD) methods only. Biyikoglu [53] presented a review about different aspects of modeling and simulation, including CFD modeling and flow field design. Tao et al. [54] presented a comprehensive review of mathematical modeling of PEM fuel cells, which especially focused on model validation and parameter influence. Djilali et al. [5] gave a critical discussion about computational strategies for the polymer 17  electrolyte membrane, porous gas diffusion electrodes, and microchannels. Multiscale strategies were also discussed. Siegel [55] presented a detailed literature overview of PEM fuel cell models with a focus on modeling strategies and commonly used model assumptions. Empirical models, mechanistic models, and computational fluid dynamics (CFD) models have been developed to study the gas-liquid two-phase flow. In PEM fuel cells, gas-liquid two-phase flow occurs simultaneously with mass transfer, heat transfer, and electrochemical reactions, and is affected by the material properties in different components. Therefore, CFD models can be effective tools for the parametric investigation of the effect of two-phase flow on PEM fuel cell performances. Liquid water transport was first incorporated in fuel cell modeling in the early 2000s, treating the liquid water as a solid species that only occupies a certain volume fraction [56] or neglecting the convective transport of liquid water [57]. As computational power increased, more complex two-phase flow models have been applied to the modeling of PEM fuel cells. In this section, several key two-phase flow models applied to PEM fuel cells are reviewed, including the multi-fluid model, mixture model, volume of fluid method (VOF), and Lattice Boltzmann method (LBM). 1.3.2.1. Multi-fluid model The multi-fluid model was first used in PEM fuel cell modeling by Berning and Djilali [58]. In this model, each phase is represented by one complete set of conservation equations (mass, momentum, and energy), and the two phases are coupled by the saturation state. This model has only a few assumptions, but has the highest number of dependent variables to solve, and the coupling of the phases can lead to unstable solutions [6].  18  1.3.2.2. Mixture model The mixture model was first used to model PEM fuel cells by Wang et al. [59] and uses the same equation set as the multi-fluid model. Each phase is modeled using an individual mass conservation equation, but a single momentum equation is solved to obtain the velocity field of the mixture, of which physical properties are the average of the two phases. Each phase velocity can then be extracted from the mixture velocity in post processing. Recently, Gurau et al. [60] commented that the mixture model was limited to flows without phase transitions or phase production because the momentum term due to phase change is neglected. For more complex situations, such as in PEM fuel cells, this model may lead to predictions of unrealistic velocity and scalar fields. Although Wang‟s group [61] responded that the “missing” term was relatively small compared with the Darcy term, Gurau [62] insisted that the missing term possibly had the same order of magnitude as the Darcy term. 1.3.2.3. Volume of fluid method (VOF) The volume of fluid (VOF) method was developed in the 1970s as a flexible way to simulate complicated free boundaries [63] and this method has become popular in simulating gas-liquid flows in fuel cell gas flow channels since it was first incorporated by Quan et al. [64]. The model can simulate immiscible fluids by solving a single set of momentum equations and then tracks the volume fraction of each of the fluids throughout the domain. Due to its capacity to consider surface tension and wall adhesion effects, liquid droplet behaviours can be captured and traced. Thus, this model is especially suited for surface tension dominated flows and flows in channels with different wall materials. However, because the specific structure of the flow domain is required this model has been only applied to gas flow channel simulations, and it is difficult to be coupled to the electrochemical reactions in the fuel cell.  19  There are two similar interface tracking methods dealing with immiscible twophase flow as VOF method, namely, Level Set (LS) Method [65, 66] and Front Tracking (FT) Method [67, 68]. The basic idea of the LS method is to represent the interface by the zero level set of a smooth scalar function Φ(x) [69], and then the position of the interface is known implicitly through the nodal values of Φ. In the FT method, interfacial locations are tracked by a set of Lagrangian marker points [69]. These two methods have been applied to various applications in mini- and microchannels [70-73]. However, they have not been applied to PEM fuel cells because of the large computational demand for simulating large and multiple channels. 1.3.2.4. Lattice Boltzmann method (LBM) Instead of solving the Navier-Stokes equations with traditional CFD methods, the Lattice Boltzmann method models the fluid as fictive particles, which perform consecutive propagation and collision processes over a discrete lattice mesh. In conventional CFD models, it is difficult to implement microscopic interactions, such as interfaces between gas and liquid phases, into the macroscopic Navier– Stokes equation. However, in the LBM, the particulate kinetics provides a relatively easy and consistent way to treat the microscopic interactions by modifying the collision operator [74]. Thus, this method shows great potential to simulate the two-phase flows in PEM fuel cells. Unfortunately, as a mesoscopic model, it is difficult to apply this method to large length scales, and the coupling of this model with heat transfer and reactions is still a challenge. 1.3.2.5. Two-phase flow model coupling with electrochemical reactions The membrane electrode assembly (MEA) also consists of porous media, with a length scale much smaller than those of the GDL. To date, because of the huge diverse length scale ratios between different layers, only macroscopic two-phase models can be employed to integrate the MEA and electrochemical reactions, and the MEA usually was assumed to be a homogenous medium.  20  Le et al. [75] developed a 3D general model of PEM fuel cells, and coupled the VOF method with electrochemical reactions, heat transfer and species transport. Liquid droplets were initially located in the serpentine gas channels of the cathode side. Furthermore, Le et al. [76] applied this model to investigate the flow behaviours of liquid water in serpentine-parallel flow channels and their effects on PEM fuel cell performance. However the VOF method using an interface tracking algorithm was unable to trace the gas-liquid interface without knowing the microstructure of porous media. Therefore, in their model, the GDL, catalyst layer and membrane were all assumed to be homogenous media. This confirms that it is not appropriate to apply the VOF method to simulate this region. 1.3.3. Two-phase flow studies in PEM fuel cells Two-phase flow in fuel cells is different from traditional two-phase flow in other applications [77]. One difference is that the water content changes along the length of the channel as water is introduced into the channels from the GDL after reaction at the catalyst surface. This introduction method means that water droplet generation and removal at the GDL surface into the channel must be considered. This issue is further complicated by the location of the emerging droplet because the water removal process depends on whether the droplet is created on the porous GDL surface towards the center or the wall of the channel [78]. The surfaces of the channel also have dissimilar contact angles, since the flow field plates and the GDL have different contact angles, which may influence droplet behaviour. Once the droplet is removed from the GDL surface, it can coalesce with droplets downstream, forming two-phase flow patterns. It has been found that the fuel cell performance depends on the two-phase flow patterns. Also, the two-phase flow in fuel cell flow channels is characterized by a large gas to liquid flow ratio and a decreasing mole fraction of the reactant gas down the length of the channel due to hydrogen and oxygen consumptions. Therefore, the two-phase flow pattern  21  may vary along the channel, which makes predicting the two-phase flow pattern very challenging. Non-uniform temperature distributions created by local hotspots can also change the amount of water that will remain in the gas phase, which affects the water balance in the flow channel [79]. Non-uniform current distribution, which changes the amount of water being produced, can also lead to non-uniform distribution of the water product in the channels. For multiple parallel channels, non-uniform current distribution also causes non-uniform gas distribution in each channel, i.e., gas flow maldistribution, which aggravates non-uniform water production which should be well avoided. In this section, experimental and numerical studies about those concerns, i.e. droplet behaviour, two-phase flow patterns, and gas flow maldistribution will be reviewed in detail. 1.3.3.1. Droplet behaviour In the fuel cell, water enters the gas flow channels from the GDL and the behaviour of these droplets is important in understanding the development of two-phase flow. Schillberg and Kandlikar [80] provided a detailed review of water droplet detachment mechanisms, summarizing relevant operating variables such as the channel dimensions, droplet sizes, Reynolds number, GDL properties, temperature, water introduction rate, and gas flow conditions. One important characteristic of droplet is the critical detachment diameter, which is the maximum droplet size on a specific surface. The critical detachment diameter on a GDL is a function of the air flow rate, water injection rate, and material contact angles. The GDL material, for instance, a standard Toray carbon fiber paper (without PTFE treatment) is highly hydrophilic, which facilitates water spreading instead of  22  droplet formation. The droplet dynamics under hydrophilic and hydrophobic conditions are schematically illustrated in Figure 1.6.  Figure 1.6: Effect of GDL hydrophobicity on droplet shape.  When the contact angle is less than 90o, the surface is considered as hydrophilic and the droplet wets the surface; when the contact angle is greater than 90o, the surface is hydrophobic and the droplet beads up on the surface. The hydrophobic cases allow a droplet to reach a critical diameter. However, when the surface is hydrophilic, the droplet does not detach and remains on the surface, blocking oxygen diffusion into the GDL and starving the electrochemical reaction. The contact angle of the channel wall is also an important parameter. Theoretically, water film formation along the channel is dictated by the Concus-Finn condition [81]:     half   / 2  (1.17)  where θ is the contact angle of water on the channel and α is the half-angle formed by the channel corner. For a rectangular channel, αhalf is equal to 45o. The wall contact angle has to be lower than 45o in order to achieve film formation along the flow channels. In PEM fuel cells, the channel walls are usually hydrophilic, and more hydrophilic channel surfaces are desired for proper water management since top wall film flow is considered to be a preferable flow pattern for water‟s removal in fuel cells.  23  Ous et al. [78] showed that the air velocity that caused detachment is inversely proportional to the droplet size, i.e., smaller droplets detached at higher velocities. Taller droplets are easier to be removed than flatter droplets due to a greater drag force relative to the surface adhesion force. Polytetrafluoroethylene (PTFE) loading causes the droplet to bead up, which can increase contact angle hysteresis and therefore more deformation occurs. Greater deformation of the droplet decreases the surface tension between the water and carbon fiber paper, leading to detachment. Temperature is also an important variable in the droplet detachment process. As the temperature increases, the surface tension decreases, which allows droplets to be removed from the GDL at lower velocities [82]. The effect of advancing and receding contact angles is also important. The definition of contact angle is in reference to a static droplet. However, when the air flows over the droplet, the droplet deforms and two contact angles are created. These are referred to as advancing and receding contact angles and are also called contact angle hysteresis, which should be considered when modeling droplet detachment from a GDL surface [12]. A schematic of the dynamic behaviour of a droplet with contact angle hysteresis is shown in Figure 1.7.  Figure 1.7: Schematic of contact angle hysteresis (A - advancing contact angle, R - receding contact angle).  24  The capillary number, the ratios the viscous forces exerted on a droplet by the air to the surface tension, has been used to characterize the deformation of a droplet. The capillary number, Ca, is defined as:  Ca   V   (1.18)  where μ is the viscosity and V is the velocity of the continuous phase (in this case air) and γ is the interfacial surface tension. Over the range of Ca from 0.014 to 0.219, droplet deformation was studied numerically on a solid surface and it was found that the deformation was a strong function of Ca when it is large [83]. Droplet detachment can also be characterized by a critical Ca, which corresponds to the point that the advancing and receding contact angles reach an observable limit [82]. These results were studied numerically and experimentally on GDL surfaces relevant to fuel cells. Besides experimental studies, Kumbur et al. [84] developed a macroscopic force balance containing relevant parameters such as the contact angle hysteresis (difference in advancing and receding contact angles), flow velocity, droplet height, chord length, and the channel height. The results of their analysis are in good agreement with experimental results, and one conclusion of this study is that for a constant droplet size and channel width, a lower channel height aids in droplet removal. Chen et al. [85] developed a two-dimensional simplified cylindrical droplet model to predict the instability of a single water droplet based on macroscopic force balances and a droplet-geometry approximation. Their qualitative results indicated that increasing the flow channel length or mean gas flow velocity, decreasing channel height or contact angle hysteresis, or making the GDL surface more hydrophobic would reduce the critical detachment diameter and enhance the removal of droplets. The same model also predicted the critical gas velocity required for a spherical water droplet to detach from the GDL surface [86]. It was found that the critical gas velocity varied inversely with water-droplet size (to the 2/3 power), and decreased with increasing GDL surface  25  hydrophobicity, decreasing contact-angle hysteresis, and lowering the channel height. However, the geometry approximation used in this simplified model would result in an inaccurate drag force on the droplet, especially at high gas flow velocity. Parametric studies with the VOF method [87] also showed that the height of the channel as well as the width of the pore had a significant impact on the detachment of the water droplet. The critical velocity was found to decrease with increasing droplet size and decreasing GDL pore diameter. More hydrophobic GDL surfaces aid in droplet detachment because of lower capillary forces, as shown by He et al. [88] using the multi-fluid model. These results are in agreement with Hao et al. [89], who applied a multiphase LBM approach to show that high gas flow velocities and a more hydrophobic GDL surface were beneficial for the water removal. An analytical model based on a force balance was also developed to predict the droplet detachment size. Zhu et al. [90] found that the critical air velocity for detachment decreased with increasing the hydrophobicity of the surface and increasing the initial size of the droplets. Temperature also has an effect on the droplet‟s detachment [82]. Experimentally measured contact angles and operating conditions were input into a numerical model based on the VOF method where water was injected from a single pipe. The results showed that higher temperatures facilitate the droplet‟s detachment due to lower surface tension and adhesion forces. From simulation, contact angle hysteresis was found to play a major role in droplet detachment dynamics. Fang et al. [91] investigated the effects of contact angle hysteresis on droplet detachment height using the VOF method and showed that contact angle hysteresis impacts slug elongation and detachment. The model agreed well with the deformation of droplets measured in microchannels. Without considering contact angle hysteresis, the droplet‟s detachment height was quite different from what was observed in experiments. These results  26  are shown in Figure 1.8. The results also implied that the contact angle distribution along the droplet can be approximated by piecewise linear functions.  Figure 1.8: Droplet detachment height versus air velocity. Both hysteresis and non-hysteresis results are presented for comparison with experimental results. (Reprinted from Chen et al. [91] with permission from Elsevier)  After droplet detachment, a liquid droplet may have different behaviour while moving along the gas flow channel because of different operating conditions or materials (different wettabilities of channel surfaces). A hydrophobic GDL surface and hydrophilic channel sidewalls, which is a common configuration in PEM fuel cells, turn dispersed droplets into thin water films attached to the channel sidewalls [92]. Shirani et al. [83] used the VOF method to investigate the motion of a liquid droplet. By studying the effects of gas velocity, the density and viscosity of water, and the surface tension on the droplet deformation, it was found that the droplet shape strongly depended on the capillary number when the capillary number was large and poorly correlated with the Reynolds number. Different fabrication techniques can alter the surface roughness, which is often only reported from the manufacturer as an average value. The surface roughness is especially important for channels with small hydraulic diameters (0.62-1.067 mm), as the pressure drop and heat transfer rate can be increased  27  with increased surface roughness [93]. The work of Dunbar et al. [94] suggested that slugs of liquid water in channels often move from surface defect to surface defect. These findings suggest surface roughness plays a role in two-phase flow in PEM fuel cells. The liquid water droplets appear in preferential areas, rather than uniformly along the channel [95]. Once detached, the droplet moves along the flow channel, where it can combine with other droplets to form slugs [96]. The droplet formation from the GDL surface in PEM fuel cells is difficult to study in ex-situ experiments or in simulation. In ex-situ experiments, water is often introduced between the middle and the end of the channels instead of uniformly along the channel. Also, droplets have been identified in two categories: land-touching and non-landtouching, and those that touch the land grow faster and to a larger size [78]. The location where water droplets are first observed also changes over long operating times (3000 hrs). Observed at 161 hrs, 2036 hrs, and 3092 hrs, the emergence of droplets moved toward the exit with time [97]. However, most exsitu two-phase flow studies are carried out over short time periods (typically < 1 hr at each data point), and conclusions drawn from ex-situ studies or computational simulation on droplet dynamics may not accurately capture fuel cells with long expected lifetimes. 1.3.3.2. Two-phase flow patterns In an operating fuel cell, two-phase flow patterns impact the pressure drop and liquid water distribution in the flow channel, which will alter the PEM fuel cell performance. Liquid water holdup is a particular concern for low Bond number (10-4≤Bo≤10-1) and low Suratman number (103≤Su≤105) environments [15]. The Bond number is the ratio of gravitational force (body force) to surface tension for a liquid surface and the Suratman number is the ratio of surface tension to viscous forces:  28   gL2   L Su  2   Bo   (1.19) (1.20)  where ρ is the density difference between phases, g is the acceleration due to gravity, L is a characteristic length such as the drop radius, γ is the surface tension, and μ is the dynamic viscosity. For low Bo and low Su conditions, the noted flow patterns are slug flow, core-annular, and transition flows [15]. Typical flow patterns in operating fuel cells can be seen in Figure 1.9 from the work of Hussaini and Wang [98].  Figure 1.9: Typical flow patterns in PEM fuel cell flow channels. (Reprinted from Hussaini et al. [98] with permission from Elsevier)  Not all two-phase flow studies show the same flow patterns and the lack of consistency highlights the difficulty in understanding and characterizing twophase flow in operational cells. Additionally, the schematic in Figure 1.9 contains stray droplets in the description of single-phase flow, which would be more accurately described as mist flow or as a pseudo-homogenous flow. Mist flow has been identified in an operating fuel cell but at an air stoichiometry of 10, which may be unrealistic for a fuel cell due to high parasitic power losses [14]. Further complicating the identification of flow patterns, fuel cells operate at  29  different relative humidity and temperatures, with different flow channel configurations, different flow rates, and different surfaces (GDL and channel). Individual results are thus noted for specific setups, and little work has been done to determine the effect of such operating conditions on two-phase flow in a general sense. In the ex-situ experiments, flow patterns were investigated under flow conditions related to fuel cell operations. Similar to the flow regimes identified in the in-situ experiments, typical flow regimes of ex-situ two-phase flow relevant to fuel cells are shown schematically in Figure 1.10.  Figure 1.10: Flow patterns in PEM fuel cells.  The gray in the figure represents liquid water while the clear areas represent air. Various researchers have observed the flow regimes shown in Figure 1.10. These include slug flow (Figure 1.10a) where a discrete droplet grows to or close to the size of the channel, blocking gas passage [15], transition flow from slug flow to annular flow (Figure 1.10b) [15], wavy-stratified flow (Figure 1.10c) [94], and stratified (or annular) flow (Figure 1.10d) that occurs at high superficial gas velocities with low pressure drop fluctuations [15]. Lu et al. [99] found that at low gas velocities (typically stoichiometric ratios below 5) slugs or semi-slugs are in dominance. They have also reported a mist flow, which is considered as an effective way to remove water because liquid droplets are dispersed in the gas 30  phase and removed convectively. However, mist flow requires high gas velocities, resulting in high parasitic power losses when applied to an operational fuel cell. Film flow or stratified flow is therefore considered to be the most favorable flow pattern for water removal from the gas flow channels because it requires a minimum gas velocity to achieve the desired flow pattern. Predicting two-phase flow patterns in simulation is challenging. One main problem is the limitation of computational time. As stated above, only VOF method is capable of tracking droplet movement at the gas flow channel scale. In an operating PEM fuel cell, the formation of two-phase flow pattern in gas flow channels usually takes several minutes or even hours due to the large gas to liquid flow ratio. However, the current computational capacity can only simulate up to several seconds. One remedy is to amplify the liquid generation rates [87, 90, 100-105]. Although Quan et al. [103, 104] stated that the unrealities introduced by the amplified water generation rate are not significant, to the author‟s knowledge, there is no proof of that argument in the literature yet. Another remedy is to locate a certain amount of liquid water or droplets in the gas channel prior to the simulation [64, 75, 76, 106]. From prior simulation work, droplet flow, stratified flow, and slug flow were usually identified as in experiment. It was found that the two-phase flow patterns strongly depend on the operating conditions, e.g., gas flow rate, channel configuration, and surface wettability. However, none of those work consider the microstructure of GDL, which determines the location of water emergence sites. Water was usually injected into the channel from one inlet or the whole GDL surface in the first method. And the second method completely neglects the water formation. Therefore, all previous simulation results were not able to capture the two-phase flow pattern accurately. Flow pattern maps are useful because it shows how superficial air and liquid velocity can be exploited to give a particular flow regime. Hussaini and Wang [98] constructed a flow map showing different regions at different superficial gas and  31  liquid velocities. Trabold [11] recommends operating the channels in the annular flow regime, which would require a superficial gas velocity of 5 to 6 m s-1. This regime allows water to be removed on the channel walls while leaving the GDL surface available for gas transport. An ex-situ flow pattern map was given by Lu et al. [99] (Figure 1.11). They found the superficial gas velocity should be more than 3 m s-1 to avoid slug flow. Flow pattern maps have been generated in terms of superficial gas velocities and superficial liquid velocities for other applications as well. The bubbly flow pattern is not observed in PEM fuel cells due to the high ratio of gas flow rates to liquid flow rates required. In micro-channels, the surface tension, inertia, and viscosity are important parameters. Different forces can be combined to form several dimensionless groups as discussed by Akbar et al. [107], which may help create dimensionless flow regime maps with greater relevance to fuel cells.  Figure 1.11: Ex-situ flow patterns in terms of superficial gas and liquid velocities. (Reprinted from Lu et al. [99] with permission from Elsevier)  32  1.3.3.3. Flow maldistribution Flow maldistribution has been indentified from both in-situ and ex-situ PEM fuel cell experiment with multiple channels as shown in Figure 1.12. A uniform distribution of current density is considered to be important in fuel cell operations because it leads to a uniform distribution of temperature and liquid water production, and lower mechanical stresses on the membrane electrode assembly (MEA) [108]. Flow maldistribution arises in PEM fuel cells when common inlet and exit headers are used to supply reactant gases to multiple channels and to collect exhaust gases with the same pressure drop across the flow field. In addition to fuel cell flow fields, flow maldistribution is also commonly encountered in other applications that utilize a large number of parallel channels such as heat exchangers [109, 110] and micro-reactors [111, 112].  Figure 1.12: Flow maldistribution in both ex-situ and in-situ PEM fuel cell experiment. (Reprinted from Kandlikar et al. [113] with permission from Elsevier)  The flow maldistribution occurs in PEM fuel cells due to different flow resistances caused by inherent geometric differences, GDL intrusion during compression [114], and uneven liquid water accumulation in the flow field channels [115]. This maldistribution is problematic because it leaves the flooded channels without sufficient air flow and the unflooded channels with excessive air flow, leading to a non-uniform distribution of current density and membrane hydration. As a result, it lowers fuel cell performance, causes pressure drop and erratic fluctuations in  33  current density, and leads to the loss of cell efficiency and lowered cell durability. It is also considered to be an important factor in reducing the operating lifetime of a fuel cell [116, 117]. Therefore, proper gas reactant distribution is critical to ensure high performance and long lifetime for a PEM fuel cell. The gas flow rate in the entrance region of individual channels can be used as an indicator of flow distribution in the flow fields, since no other measures are available to characterize this effect. Kandlikar et al. [113] developed an entrance region pressure drop measurement technique to determine instantaneous gas flow rates through individual channels. The method was employed in both an exsitu and in-situ experimental set-up, and it was found that a porous GDL backing could lead to severe flow maldistribution compared to impermeable backing for the same channels. At in-situ operating conditions, flow maldistribution was also observed due to water blockage in gas flow channels. One limitation of this method is the difficulty for implementation in operating fuel cells. Nevertheless, it still can provide some valuable information of flow maldistribution in parallel channels related to other fuel cell operating parameters such as current density, gas stoichiometry, and gas humidity in a fuel cell with specially designed gas introduction. Zhang et al. [118] employed a similar method to measure flow distributions in parallel channels and it was found that flow distributions also depend on whether the gas flow rate is changed in an ascending or descending manner. In-situ (electrochemically active) experiments also showed significant maldistribution among parallel channels at 35 oC. Low temperature and low current densities (low gas flow) are particularly relevant conditions for automotive applications [119], and are also the conditions prone to flooding and therefore maldistribution. Realistic air stoichiometry (defined as the ratio of the number of moles of oxygen at the inlet to the theoretically needed number of moles of oxygen), namely 1~5, gives rise to slug/semi-slug flow patterns in the flow channels, leading to severe maldistribution in parallel channels and large fluctuations in the pressure drop [99]. Basu et al. [114] used numerical simulation to investigate maldistribution in 7 parallel channels with GDL intrusion and  34  showed that the edge channels were susceptible to air stoichiometry less than 1, which would be insufficient to sustain the electrochemical reaction. Models have also shown flow maldistribution among multiple channels in interdigitated flow fields, where high GDL permeability tends to allow for greater maldistribution [120]. The maldistribution also varies in parallel channels if the inlet and exit headers are connected in a U-type or Z-type configuration [121]. Basu et al. numerically explored the effect of GDL intrusion and manifold design [122]. It was found that an improved manifold design with a flow splitter could reduce the maldistribution by more than 50%. Multiple bifurcations have been proposed In order to reduce flow maldistribution [123]. It was found that the length and shape of bifurcations has great impact on the flow distribution uniformity. Novel flow field designs such as a bionic flow slab have been proposed numerically to increase flow uniformity while maintaining a reasonable pressure drop [124]. A slotted-interdigitated design has also been optimized to reduce maldistribution but the optimization method did not consider flooding-induced maldistribution [125]. Despite many research efforts on reducing flow maldistribution in either in-situ or ex-situ fuel cell flow fields, little attention has been paid to the effect of communication between parallel channels on flow distributions.  These  communications are important since the requirement of an equal pressure drop and different flow resistance in each parallel channel unavoidably lead to flow communications through the porous gas diffusion layers. In a single-phase flow system, flow distributions only rely on the flow resistance of the gas flow channels and communication in parallel channels enables flow redistribution when there is a pressure difference between adjacent channels. In contrast, in a two-phase flow system (due to the presence of liquid water), both gas and liquid water can be redistributed due to the pressure difference between adjacent channels. The redistribution of two-phase flow can produce two distinct impacts on flow distributions in PEM fuel cells, a beneficial one and an adverse one. The  35  beneficial one is to enable more even distribution of both gas and liquid and the adverse one is to create more severe maldistribution by flooding one channel with a significant presence of liquid water and drying other channels with an excess amount of gas flow. In the second situation, the overall performance of fuel cells can be significantly lowered. Therefore, an understanding of two-phase flow distribution in communicating channels is required to design optimal fuel cell flow fields and to achieve stable and durable operation of PEM fuel cells. The occurrence and the recoverability of channel blockage in parallel channels are directly related to the instabilities of the two-phase flow in parallel flow channels and are affected by many parameters such as the channel flow rate and the wall physical properties such as the contact angle. Theories or models to interpret instability-induced flow maldistribution are still lacking in the open literature. An attempt was made recently to analyze the stability of possible solutions of gas and liquid flow distributions in parallel channels for fuel cells with a one-dimensional momentum balance equation across the channels [126]. All possible combinations of gas and liquid flow distributions must satisfy the equal pressure drop across all channels if they share a common inlet and outlet. Theoretically, even flow distribution is one default solution of the equal pressure drop requirement, but experimental results show that an even distribution is not always observed. Instead, flow maldistribution appears as a stable solution, indicating that flow distributions of two-phase flow in parallel channels depend on not only pressure drop but also flow stability. If the two-phase flow is operated at an unstable condition, a small perturbation will shift the flow to the nearest stable condition. A more rigorous theoretical analysis should be conducted over a wide range of operating conditions in the future.  1.4. Research objectives and thesis layout As stated in the literature review, several issues were identified that need further investigation, i.e., multiple droplets behaviours when GDL microstructure is considered, two-phase flow patterns in gas flow channels, impact of two-phase 36  flow on PEM fuel cell performances, impact of flow mal-distribution on PEM fuel cell performance, and maldistribution mitigation. Therefore, the objectives of this work are to investigate the gas-liquid two-phase flow phenomena in single and multiple gas flow channels and its impacts on the PEM fuel cell performance. To achieve those objectives, following work has been carried out. To study the impact of two-phase flow on PEM fuel cell performances, a socalled 1+3D model is presented in Chapter 2. This model deals with two-phase flow in the cathode side only.  A 1D MEA model is adopted to simplify the  transport phenomena across the GDL, catalyst layer and membrane, and the two-phase flow in the 3D gas flow channel is simulated by the VOF method. Multiple droplets behaviour in a small-scale three-dimensional gas flow channels was firstly modeled by VOF method using a commercial CFD software FLUENT. The non-uniform surface of gas diffusion layers was considered by creating several pores on the GDL surface to represent the real microstructure of GDL surface. Different GDL surface structures were also simulated by changing the pore size and pore number. Parametric studies such as the influences of wettability of channel surfaces, influences of amplified liquid flow rate, were carried out. This part of work will be presented and discussed in Chapter 3. To study the two-phase flow patterns in a large-scale gas flow channel, the work in Chapter 3 is extended by implementing the simplified microstructure to a much larger channel, which is comparable with realistic PEMFC gas flow channels. The effects of liquid injection rates and surface wettability on two-phase flow patterns in the flow channel are investigated and simulation results are presented in Chapter 4. The impact of flow maldistribution in parallel channels on two-phase flow pattern was simulated in Chapter 5. The gas maldistribution was introduced by injecting water into the parallel channels with different flow rates. Two-phase flow patterns  37  under different degrees of gas maldistribution were compared. The impact of GDL surface and channel wall wettability on the gas distribution and GDL water coverage ratio was also simulated. In Chapter 6, in order to develop maldistribution mitigation method, two parallel flow channels connected by several communication channels were designed and tested in simulation. Different number and/or width of the communication channels were simulated under severe maldistribution conditions. Both steadystate and transient flow behaviours were simulated and discussed. All of above work only considered the two-phase flow in gas flow channels without involving electrochemical reactions in the MEA. In Chapter 7, the influence of two-phase flow in a single channel on the fuel cell performance was investigated using the 1+3D model developed in Chapter 2. The PEM fuel cell output voltage and pressure drop were compared under different gas stoichiometric flow ratios, water generation rates and surface contact angles. The effects of flow maldistribution on fuel cell performance was studied in a unit consisting of two parallel channels and presented in Chapter 8. Effects of flow maldistribution and liquid generation rate on PEM fuel cell performance were tested. Several maldistribution mitigation methods were also examined in simulation, including gas flow rates, material surface wettability and adding inlet resistance. To evaluate the simulation results, two experiments were carried out in Chapter 9. In one experiment, the two-phase flow pattern in a transparent gas channel under different GDL surface was recorded. In another experiment, the effects of maldistribution on two-phase flow in two parallel channels, and the impact of communication channels on the gas maldistribution were conducted. The degree of flow maldistribution with or without communication channels was compared.  38  In the last chapter, all the work of this thesis is summarized. Conclusions and recommends for future studies are given based on the findings of this research.  39  2. Development of two-phase flow model of PEM fuel cells The two-phase flow model adapted in this thesis is called 1+3D model, which couples 3D volume of fluid (VOF) model with a 1D MEA model. The VOF part was implemented in commercial software, FLUENT® 6.3.26, and the 1D MEA model was coupled with the FLUENT through its user defined functions written in C programming language. In Chapters 3, 4, 5, 6, only VOF model was used to simulate the droplet behaviour and two-phase flow patterns in a single channel and gas maldistribution in two parallel channels. In chapter 7 and 8, the whole 1+ 3D model was adapted to study the effects of two-phase flow and gas maldistribution on the fuel cell performance.  2.1. 1+3D model The proposed model only deals with two-phase flow in the cathode side. As shown in Figure 2.1, a 1D MEA model is adopted to simplify the transport phenomena across the GDL, the catalyst layer and the membrane, and the twophase flow in the 3D gas flow channel is simulated by the VOF method. The MEA model and the 3D flow channel model are coupled at the GDL surface which lies between the gas flow channel and the MEA. Several pores at the GDL surface are used to represent the large pores inside the GDL where liquid water emerges.  40  Figure 2.1: Illustration of 1+3D cathode side of PEM fuel cell.  Specifically, as shown in Figure 2.2, the temperature, oxygen concentration and gas volume fraction at the bottom wall of the gas channel, i.e. at the GDL surface, are input into the 1D MEA model as its boundary conditions. Then the overall voltage is calculated for the given boundary conditions and other parameters. Finally, the flux of oxygen, liquid water and water vapor are calculated based on the mass balance, and returned to the gas channel as mass sources at the GDL surface.  41  Figure 2.2: 1+3D model coupling method.  2.2. 1D MEA model The 1D MEA model in this work is derived from Berg‟s 1D MEA model [127], except that the membrane is assumed to be fully humidified, which is justified when the purpose of the model is to examine effects outside the membrane (e.g. cathode flooding) [51]. Berg‟s model has been validated against their experimental results by comparing the polarization curves obtained from experiments and simulations. It was proved that the 1D MEA model is capable of accurately predicting the PEM fuel cell performance at given operating conditions. In 1D MEA model, the local potential U is calculated by,  U  E0 c   I  (2.1)  where the anode overpotential is neglected due to its relatively small value compared to the cathode overpotential. From the Butler-Volmer equation, the local current density I is given by,   C cl I  io,c  oref  Co    c F   (1   c ) F   c   exp   c    exp  RT RT        (2.2)  42  At high current densities, the overpotential is an explicit function of the current density and the oxygen concentration in the membrane, which is approximated by its Tafel form.  c   RT  Co ref I ln   c F  io,cCo cl      (2.3)  where Co cl is the oxygen concentration in the catalyst layer, which can be expressed via Henry's law,  Co cl   RT m Co Ho  (2.4)  where Co m is the oxygen concentration at membrane, which can be calculated from the mass balance.  J O2    Do2 eff LG  (Co2  Co2 m )    I 4F  (2.5)  Define the oxygen mass transfer coefficient  as,    LG 4 Do2 eff F  (2.6)  The cathode overpotential then becomes,   H oCo ref I RT  c  ln    c F  RTio,c (Co2   I )   (2.7)  where temperature, T , and oxygen concentration, Co2 , are the boundary conditions from CFD calculations. The overall ohmic loss is calculated by,  43  I  (2.8)  where  is the averaged cell conductivity. Finally, the local potential can be expressed as,   RT  Co ref H o I U  E0  ln    I  c F  io,c RT (Co2   I )   (2.9)  Equation (2.9) gives the relationship between potential and current density at a local position, which varies along the channel due to the variation of oxygen concentration. Therefore the total current density at a given voltage is the integration of local current density weighted by the mesh area. In reality, the total current density is fixed in the experiment instead of the fixed voltage. Thus, a feedback program is added to maintain a specific overall current density. After knowing the local current density, the oxygen and water mole fluxes can be calculated according to the mass balance,  IS A  gas 4F IS  (1   ) A  gas 2F  NO2   N H 2O  (2.10) (2.11)  where S is the area of a local mesh at the GDL surface. Since both liquid water and gas mixture exist in the channel, it is assumed that the liquid covering the GDL surface will block the diffusion path of gases from the channel through to the membrane, by assuming that the gas diffusion through the liquid layers is much slower than through the gas mixtures. Therefore the effective area for the reactants diffusion is reduced, and so is the effective reaction area at the catalyst layer. Another assumption is that all the water generated by the electrochemical reaction and the electro-osmotic flow diffuses back to the cathode side channel only, ignoring the water crossover term, e.g. 44  the diffusion of the liquid water from cathode to anode, since it is relatively small in comparison to the value in equation(2.11).  2.3. 3D VOF method The two-phase flow hydrodynamics in the 3D gas flow channel is simulated by the VOF method. In this method, a single set of momentum equations are solved to obtain the volume fraction of each immiscible phase, with the interface between phases being tracked throughout the domain. Conservation equations include, Mass conservation equation,     (  v)  Sm t  (2.12)  where S m is the mass sources and for gas-liquid flow, the mixture density is defined as,     liquid liquid   gas  gas  (2.13)  where  is the volume fraction of each phase. The volume fraction of each phase is given by the following relationship,   liquid   gas  1  (2.14)  Momentum conservation equation, T  (  v)    (  vv)  p    (  (v  v )   g  Fvol t  (2.15)  where p is the static pressure and μ is the dynamic viscosity given by,     liquid liquid   gas  gas  (2.16)  45  Volume fraction conservation equation is given by, ( liquid liquid ) t     ( liquid liquid v)  0  (2.17)  The surface tension effect is considered by the continuum surface force (CSF) model proposed by Brackbill et al. [128]. In this model, extra surface tension results in a source term in the momentum equation (2.15) given by,  Fvol   surf (   liquid  liquid 1 ( liquid   gas ) 2  )  (2.18)  where  surf is the surface tension coefficient, and  liquid is the surface curvature of the liquid droplet defined in terms of the divergence of the unit normal nliquid , given by,   liquid    nliquid  (2.19)  The unit normal nliquid is computed from local gradients in the surface normal at the interface.  nliquid    liquid  liquid  (2.20)  Wall adhesion effects are considered by adjusting the curvature of the surface near the wall, where the gas-liquid interface meets the solid surface. The local curvature of this interface is determined by the static contact angle,  w , which is the angle between the wall and the tangent of the interface at the wall. The surface normal, n , at the wall is given by,  n  nw cos  w  t w sin  w  (2.21)  46  where n w and t w are the unit vectors normal and tangential to the wall, respectively. Species conservation equations,  ( Yi )    (  vYi )    ( J i )  Si t  (2.22)  where i indicates different species. In the cathode side gas flow channel, it could be water vapor, nitrogen, and oxygen. J i is the diffusion flux of species i , which is given in Maxwell-Stefan equations [129]. Si is the species source term which can be obtained from 1D MEA model as mentioned previously. The temperature in the channel is assumed to be isothermal, therefore the energy conservation equation is not considered in this model. All of the above governing equations are implemented in a commercial software, FLUENT® 6.3.26 [129]. The geometric reconstruction scheme from the work of Youngs [130] is used to represent the interfaces between two fluids, which is one of the accurate methods for interface tracking .The 1D MEA model was coupled in the FLUENT through its user-defined functions written in C programming language.  2.4. Water treatment How liquid water enters the channel is important, since it determines the droplet size and behaviour. According to experimental observations [82], the emergence of droplets from the GDL surface is not homogenous. Large pores in the GDL are preferred sites for droplets. To address this issue, a series of pores are randomly placed at the GDL surface along the channel so as to represents the droplet emerging sites, as shown in Figure 2.1. The liquid water generated by the reaction does not directly enter the channel at its local boundary but from the 47  nearest pore instead. Therefore, the liquid flow rate at each pore is the summation of liquid flux nearby. The water flux calculated in the 1D MEA model contains both liquid and vapor water. Although the heat transfer model is not coupled in the model, the water evaporation rate in the gas channel is usually fast enough that gas is saturated if liquid water is present there [51], and the minimum heat released from electrochemical reaction, as shown in section 1.2.2.1, can provide enough heat required to evaporate all liquid water generated at catalyst layer. Therefore, we assume that the water generated by reaction firstly humidifies the gas to saturation. This amount of water directly enters the channel at its local boundary, and if there is still water left, it will emerge from the nearest pore as liquid. On the other hand, if water vapor somewhere in the channel is oversaturated, the excess water will condense, and be added into the nearest emerging pores as well. Table 2.1 summarizes all key assumptions made in the current 1+3D model.  Table 2.1: Assumptions of 1+3D model. Anode overpotential and hydrogen and water transport are not considered. Membrane is fully humidified. Transport phenomena across the GDL, the catalyst layer and the membrane are simplified into a 1D MEA model. Liquid water emerges from several randomly located pores at the GDL surface. Liquid covering the GDL surface blocks the diffusion path of gases from the channel through to the membrane, giving zero gas diffusion rate through liquid blocked GDL surface. Water crossover is ignored. The temperature in the channel is uniform. Liquid water enters the channel from the nearest pore. Water generated by reaction firstly humidifies the gas to saturation.  Oversaturated water vapor condenses, and is added into the nearest pores.  48  2.5. CFD algorithm and solution procedure In all simulations, the computational grids are generated by using Gambit ® 2.4 mesh generation software. The governing equations are discretized by Finite Volume Method (FVM). The Pressure-Based Segregated Solver and PISO algorithm are adapted for solving the governing equations. Adaptive time step method is utilized in order to keep the global courant number less than 0.5, which ensures that the droplet behaviour in the channel can be well captured. The convergence criterion, mesh dependency and VOF validation were examined in a case from chapter 3. As shown in Figure 2.3, a three dimensional cuboid channel as the base case has a 250μm × 250μm square cross section and a 1250 μm length. Air flows into the channel from one end and liquid water was injected from 16 open pores on the bottom wall. Each pore has a diameter of 50 µm. The gas velocity was set at 10 m∙s-1, and the liquid injection velocity was set at 0.0625 m∙s-1 for all the pores. A total of 82,525 meshes with a uniform size of 0.1 mm are used in the simulation.  49  Figure 2.3: Three dimensional computational domain.  2.5.1. Precision improvement The default residual criterion for each governing equation was set to 10-3, and the program was run in the single-precision mode for all out early cases. To increase the accuracy, a lower residual criterion (10-6) and the double precision solver were tested, respectively. Simulation results (Figure 2.4) showed that using the lower residual criterion or double precision solver does not change the over-all droplet behaviour. However, according to the governing equations and the geometry used, water droplets emerged into a channel is expected to be symmetrical, which is not seen in Figure 2.4.  50  Default Double precision Lower residual criteria Figure 2.4: Effects of solver precision on droplet behaviour.  2.5.2. Mesh dependency Mesh size is a critical parameter in numerical simulation. Small mesh size would increase the accuracy, but also exponentially increase the computational time. Thus, appropriate mesh size should be carefully chosen for both representing the fluid characteristics and saving the computational time. The original mesh size for the base case is 10 µm, and the total number of elements is 82,525. To examine the mesh size and mesh symmetry effects, I first refined the mesh size to obtain a more symmetric profile, and then further refined the mesh size. Figure 2.5 shows different meshes for base case, improved case (X1), and further refined cases (X2 and X4). For improved case (X1), the mesh size is 10 µm, which is the same as base case. For X2, the mesh size is 8 µm and the total meshes are 155,434, almost doubled the number of meshes of the base case. For case X4, the mesh size is 6 µm and the total meshes are 364,812, about four times of the base case.  51  Y  Z  Base case  Y  X  Improved mesh (X1)  Z  X  Refined mesh (X2)  Refined mesh (X4)  Figure 2.5: Top view of meshes for base case, improved case X1, and refined cases X2 and X4  The water coverage ratio on the two side walls for different cases is shown in Figure 2.6. It is clearly seen that with improved meshes (X1), the symmetry is improved, but the difference between the two walls still exists. Meanwhile, with finer meshes (X2, X4), the symmetry is even getting worse. One possible reason is that the water inlets are circular, while the channel wall is rectangle. Thus, it is difficult to make symmetrical meshes for both the rectangle channel and the circular water inlet holes.  52  Base case  Improved mesh (X1)  Refined mesh (X2)  Refined mesh (X4)  Figure 2.6: Time evolution of the water coverage ratio on the two side walls  In order to create completely symmetrical meshes, the shape of water inlet was changed from circle to square (Figure 2.7). For case X1, the mesh size is 10 µm, and the total meshes are 117,600. For case X2, the mesh size is 8 µm and the total meshes are 230,300. For case X4, the mesh size is 5 µm and the total meshes are 940,800.  53  Y  Z  Geometry  X1  X2  X  X4  Figure 2.7: Top view of symmetrical meshes.  The water coverage ratio on the two side walls for different cases is shown in Figure 2.8. It is seen that symmetrical flow distribution along the channel can now be obtained when the meshes are symmetrical for both the channel and the water inlets and also fine enough (case X4). However, it requires extraordinarily large computational time -- four weeks to obtain one result.  54  Case X1  Case X2  Case X4  Figure 2.8: Time evolution of water coverage ratio on the two side walls.  Using finer mesh definitely increases numerical accuracy, but also greatly increases computational time. It was also found that the mesh size at the boundary has more impact on the droplet behaviour than the mesh size in the central. Therefore, in order to reduce the computational time, non-uniform size mesh, which has finer meshes near the wall, but coarser meshes in the center, has been utilized in this study. Figure 2.9 shows the side views of two types of non-uniform meshes, “X1_5” means the total meshes are the same as case X1, and the mesh size adjacent to the wall is now 5 µm. Thus, the computational time is almost the same as case X1.  55  Case X1  Case X1_5  Case X2_2  Figure 2.9: Side views of non-uniform meshes  It was found that for circular water inlet, non-uniform size mesh is able to predict the water distribution and pressure drop as good as the case with much finer mesh, meanwhile, saving the computational time significantly. As discussed above, symmetrical simulation results will be obtained when the meshes are symmetrical and fine enough. However, it takes very long CPU time to run the simulation. It was also found that some quantitative variables, such as total pressure drop, water distribution and water volume fraction do not change too much with mesh refinement. Therefore, to some extent, we can still trust that the simulation results at large mesh sizes are still able to represent the major characteristics of fluid flow. Therefore, non-uniform mesh with moderate size was still used in all simulations in this study, which ensures the computational time is acceptable and the results are creditable. 2.5.3. VOF validation A lot of studies have been done on VOF method validation in various applications. Two very relevant studies have been focused particularly on the water droplet movement and deformation in fuel cell channels [91, 131]. Quantitatively, their results showed good agreement between VOF simulation results and experimental data in both slug curvature and detachment length for different droplet size and different GDL materials, and both of their work considered the contact angle hysteresis.  56  Fang et al. [91] also found that for the non-hysteresis cases, the liquid slug detachment is much earlier than that in the hysteresis simulation as well as in the experiment. The reason is likely that the absence of the hysteresis precludes the asymmetry of the contact angle distribution and the existence of the retentive force as well, leading to the early detachment. As shown in Figure 1.8, the hysteresis model yields much better results (solid line) than the non-hysteresis model (dashed line). Both Theodorakakos et al.[131] and Fang et al. [91] used their own codes to implement VOF method. The commercial FLUENT software VOF code, however, does not include the contact angle hysteresis by default. Figure 2.10 shows that non-hysteresis model results predicted by FLUENT in this study are well in agreement with their results, which suggest that the VOF code in FLUENT performs similarly to the in-house VOF code of Fang, when hysteresis model is not included. However, after searching and trying for several months, unfortunately, it is found impossible to add contact angle hysteresis model in FLUENT unless its source code can be accessed. Therefore, only qualitative results can be obtained from using FLUENT, which is still meaningful when investigating multiple droplet behaviours since the effect of single droplet behaviour to some extent may be canceled out in a large flow channel.  Figure 2.10: Comparison of FLUENT simulation results with literature results. (Modified from Chen et al. [91] with permission from Elsevier)  57  3. Droplet behaviours in gas flow micro-channel1 As discussed in the literature review, most work using VOF method only focused on single droplet behaviour or several droplets initially located in the gas channel; but none has considered the effects of microstructure of the GDL surface. The GDL surface has been treated either as a homogenous surface or a surface with only one pore opening for liquid water injection into the channel. The multiple droplet behaviour in the flow channels of a real operating fuel cell is likely quite different from that determined based on droplets emerging from a uniform GDL surface, or from single droplet behaviour. However, full consideration of detailed microstructure of the GDL, the so-called pore-scale model, needs extremely large computational time, and has been only employed for very small length scale systems [40, 89, 132-136]. The effect of GDL microstructure on the multiple droplet behaviour in the gas flow channel is still poorly understood. In this section, we simplified the microstructure of the GDL by creating a number of representative pores on the 2D GDL surface. Droplet formation and motion were then simulated using the VOF method. Different GDL surface structures were then simulated by changing the pore size and pore number. Parametric studies were also conducted to investigate the effect of GDL and channel walls wettability and liquid velocity on the two-phase flow pattern and pressure drop through the flow channel.  1  A version of this chapter has been published: Y. Ding, H.T. Bi, D.P. Wilkinson. Threedimensional numerical simulation of water droplet emerging from a gas diffusion layer surface in micro-channels. J. Power Sources 2010; 195: 7278-7288. Part of results were presented in the conference: Ding, Y., H.T. Bi, and D.P. Wilkinson. 3-D Numerical Simulation of Water Droplet Emerging from a Non-uniform Gas Diffusion Layer Surface. in 20th International Symposium on Transport Phenomena (ISTP-20). 2009. Victoria, BC.  58  3.1. Computational domain and boundary conditions The three dimensional computational domain is shown in Figure 3.1. This cuboid channel, which was used in all the cases in this chapter, has a 250μm × 250μm square cross section and a 1250μm length with a hydrophobic GDL surface on the bottom. Air flows into the channel from one end and liquid water was injected from several open pores on the GDL surface.  Figure 3.1: Three dimensional computational domain.  For the base case, 16 pores were located within a 250μm×250μm square on the GDL surface, each with a pore diameter of 50µm, which corresponds to typical GDL surfaces with a porosity of 0.5 [6]. At the gas inlet, the velocity of air was set at 10 m∙s-1, which is of the same order of magnitude as flows encountered in automotive fuel cell stacks [90]. The theoretical calculation of liquid generation rate for this channel with 3.125×10-7 m2 of reaction area at a current density of 800 mA∙cm-2 is 0.20 mg∙hr-1. However, in order to shorten the time for water accumulation in the channel, the liquid injection velocity was set at 0.0625 m∙s-1 for all the pores with a corresponding mass flow rate of 7 g∙hr-1. Quan and Lai 59  [103] reported that deviations introduced by amplifying the water generation rate by several orders of magnitude are insignificant due to the large flow rate ratio between gas and liquid water. Laminar flow and no-slip boundary condition were assumed since the Reynolds number of each phase is quite small (Reg=171, Rel=7.8). The static contact angle for the GDL surface and channel wall surfaces were set at 140o and 45o respectively, based on typical values for PTFE treated carbon paper GDL materials [6]. There are two side walls and one top wall in the channel in addition to the bottom GDL surface. The time step for the baseline simulation was set at 10-7s to ensure that the global courant number is less than one. The mesh size was set to 10μm, and there were 82625 meshes in total for the base case. Smaller mesh sizes, e.g. 8μm and 5μm, were also tested, and showed no obvious difference from the base case. Therefore, the mesh size for all the other cases was also set to 10μm to save computational time.  3.2. Results and discussion Due to small length scale of the channel, the surface tension effects were more noticeable. For base case conditions, the capillary number (Ca), the Weber number (We), and the Bond number (Bo) were much less than one, which means the surface tension force was the dominant force in comparison to the viscous force, inertia force, and gravity force in the mini flow channel. 3.2.1. Droplet evolution Figure 3.2 clearly shows three stages of two-phase flow pattern evolution in the channel: emergence and merging on the GDL surface, accumulation on the side wall, and detachment on the top wall. In the emergence and merging stage, small droplets emerging from each pore coalesce on the GDL surface, which generates a liquid film covering the injection area. At the same time, liquid water is continually accumulating in the channel, and due to the wettability difference between the GDL surface and channel walls, liquid water tends to accumulate on the side walls which results in the formation of large droplets. As these droplets  60  “grow” larger and larger, they begin to move slowly along the side walls. Once these droplets reach the top wall, the liquid water rapidly spreads on the top wall, resulting in a thin film on the top wall. As this liquid film moves outward quickly, some of water on the side wall is also dragged away. In the last stage, the top liquid film detaches itself from the droplets on the side walls due to its faster speed, and a new cycle begins. It is worth noting that the emergence and merging stage and the accumulation stage occur continuously in the channel, but the detachment stage only occurs periodically.  (a)  (b)  (c)  Figure 3.2: Three stages of the emerging water droplet into a channel: (a) emergence and merging, (b) accumulation, and (c) detachment.  These three stages of dynamic droplet behaviours can also be identified by the time evolution of the water coverage ratio on the different surfaces of the channel,  61  the water volume fraction and the total pressure drop (Figure 3.3). The water coverage ratio is defined as the ratio of the surface area covered by water to the total surface area. In PEM fuel cells, reactants diffuse through the GDL, hence, the GDL surface water coverage ratio is a key parameter that indicates the negative effects of liquid water on PEM fuel cell performance. The water volume fraction or degree of saturation indicates the degree of channel flooding. The pressure drop, which is another important parameter, indicates the energy loss of fluid flowing through the channel. For operating PEM fuel cells, a lower GDL surface water coverage ratio, a lower water fraction, and a lower pressure drop are preferred.  Figure 3.3: Time evolution of water distribution and pressure drop.  From Figure 3.3, it can be seen that the water coverage ratio on the GDL surface continuously increases at the beginning (up to 4 ms), due to the emerging small droplets. At the same time, liquid water also accumulates on the side wall, and the higher water coverage on the side wall is due to the different wettabilities between the side wall and GDL the surface. The pressure drop also increases,  62  corresponding to the increase in the water volume fraction. After about 4 ms, a sudden drop in the GDL surface water coverage ratio occurs, followed by a steep increase in the top wall water coverage ratio, indicating the droplet's fast spread to the top wall. The pressure drop also decreases rapidly, even though the level of water saturation is still increasing. This is because the buildup of droplets on the hydrophobic GDL surface occupies more cross sectional area. On the other hand, water on the hydrophilic channel walls forms a thin film, which occupies much less cross sectional area, imposing little influence on the pressure drop even though the volume of this liquid film is much higher than the droplets on the GDL surface. The liquid film on the top wall moves faster than the droplets on the side walls, and soon (at 5.5 ms) this film detaches from the droplet, resulting in a maximum top wall water converge ratio and a minimum GDL and side wall water coverage ratio. Once the detached liquid film is flushed out of the channel, a new cycle begins. From the flow pattern shown in Figure 3.3, it can be concluded that the pressure drop does not correspond to the level of water saturation, but to the size of the droplets located on the GDL surface. Thus, these results indicate that most empirical equations [137-141] which correlate pressure drop with water fraction should not be applied to PEM fuel cell channels. 3.2.2. Effects of water inlet structures GDLs, which are typically fabricated from carbon fibers, connect catalyst layers and gas flow field and provide a pathway for electrons, gases, and liquid water to move to and from the catalyst layers. The microstructure of the GDL plays an important role in water management. In this section, several pore structures are simulated to examine the effect of the microstructure of the GDL surface on the flow patterns in order to identify a simple but representative structure for future use in the simulation of large-scale fuel cell stacks. Five structures, (uniform inlet, 1-pore inlet, 4-pore inlet, 16-pore inlet, and 64-pore inlet), with the same open area and liquid flow rate, were selected and are shown in Figure 3.4.  63  Figure 3.4: Different water injection inlet structures representing the GDL: (a) Uniform, (b) 1-pore, (c) 4-pores, (d) 16-pores (base case), and (e) 64-pores.  The uniform inlet case (Figure 3.4a) is a commonly used approximation in many CFD simulations. For this case, the microstructure of the GDL is neglected, and a very distinct two-phase flow pattern is observed in the channel, as shown in the three snapshots of flow pattern evolution in Figure 3.5a. Only one big droplet is formed during the merging stage, and the liquid accumulation occurs simultaneously on the GDL surface and side walls. As a result, more liquid stays on the hydrophobic GDL surface, which slows down the droplet spread on the hydrophilic side walls and makes the liquid droplet detach before it reaches the top wall, leading to the formation of liquid slugs.  64  (a)  (b)  (c)  (d) Figure 3.5: Effects of water inlet configurations on the two-phase flow patterns: (a) Uniform, (b) 1pore, (c) 4-pore, and (d) 64-pore.  The 1-pore case (Figure 3.4b) is another commonly used configuration to investigate the droplet behaviour in PEM fuel cell channels. Similar to the uniform inlet case (Figure 3.5b), only one droplet is formed in the channel. However, this droplet is not big enough to touch either the side walls or the top wall before detachment, leading to droplet flow in the channel. The two-phase flow patterns  65  with 4 pores and 64 pores (Figure 3.5c, d) are very similar to the case with 16 pores, with three stages of water dynamic behaviour present. The effects of pore structures on the flow patterns can also been seen clearly from the time-averaged water distribution shown in Figure 3.6. The uniform inlet and 1-pore cases have a higher water coverage ratio on the GDL surface and a much lower coverage ratio on the channel walls, compare to the 4-pore and 64pore cases which have very similar water distribution with the base case of 16 pores due to the similar two-phase flow pattern in the channel.  Figure 3.6: Effects of water inlet configurations on the time-averaged water distribution.  The effects of pore structures on the time-averaged pressure drop are shown in Figure 3.7. The uniform inlet and 1-pore inlet cases show a much greater pressure drop than the three multi-pore cases which have a lower and similar pressure drop. As discussed previously, the pressure drop is dominated by water residing on the GDL surface, which creates more blockage of the gas channel than water on the other channel walls. Thus, the pressure drop of both the  66  uniform inlet and 1-pore cases is much higher, which results in faster liquid removal, and thus lower average water volume fraction in the flow channel, as shown in Figure 3.6.  Figure 3.7: Effects of water inlet configurations on the time-averaged pressure drop.  In a real fuel cell, droplets emerge from preferential sites on the GDL surface. At the front of the channel, small droplets are formed and the flow pattern is mainly droplet flow as shown in the 1-pore case. As these small droplets move along the channel, droplets coalesce and form larger droplets which are able to touch side walls. As a result, the flow pattern in this region would be similar to that shown by the multiple pore cases. Since the GDL pores are quite small and non-uniform in size and the pores are randomly distributed on the GDL surface, a full consideration of the detailed microstructure of the GDL surface is almost impossible. However, the previous cases demonstrate that when the pore size is small enough (e.g., 4-pore case), the flow pattern changes little with a further increase in the number of pores or a reduction in the pore diameter. Thus, the 4pore case could be considered as a minimum required number of pores in order  67  to capture the two-phase flow pattern in the fuel cell mini channel at a reasonable computational time demand. 3.2.3. Effects of GDL surface contact angles The wettability of the GDL, which is characterized by the surface contact angle, can be altered by varying the PTFE content of the GDL. Although the GDL wettability has been shown to have a significant impact on the liquid water transport inside the GDL [142], its impact on the two-phase flow pattern in the gas channel is still unclear. In this section, the effects of GDL surface wettability on the flow pattern in the mini channel were investigated by varying the water/GDL surface static contact angle from 0o to 180o, i.e., 0o, 30o, 45o, 60o, 90o, 120o, and 180o. Figure 3.8 shows different two-phase flow patterns in channels for different GDL wettabilities. When the GDL surface contact angle is less than that of side walls (Figure 3.8a, b), a thin liquid film forms on the GDL, which covers almost the entire GDL surface. For a hydrophilic GDL surface with a contact angle less than 90o (Figure 3.8c, d, e), a liquid water film still tends to be formed. This liquid film blocks the diffusion pathway of gas reactants to the catalyst layer, leading to decreased fuel cell performance. Since the GDL is hydrophilic, some water may even flow back from the channel towards the catalyst layer, which would further decrease the fuel cell performance. As the contact angle of the GDL surface increases, liquid water begins to accumulate on the side walls and the higher the contact angle, the more water moves from the GDL surface to the side walls. When the GDL surface is hydrophobic with a contact angle greater than 90o (Figure 3.8f, g), liquid water droplets form on the GDL surface, and grow until they are detached. It should be noted that the flow patterns on hydrophobic GDL surfaces are all similar to the base case in section 3.2.1.  68  (a)  (b)  (c)  (d)  (e)  69  (f)  (g) Figure 3.8: Effects of GDL surface wettabilities on the two-phase flow patterns:(a) θ=0o, (b) θ=30o, (c) θ=45o, (d) θ=60o,(e) θ=90o, (f) θ=120o, and (g) θ=180o.  Hydrophobic and hydrophilic GDL surfaces show very distinct liquid distribution profiles as shown in Figure 3.9. For the hydrophilic GDL surface, the water coverage ratio on the side walls changes very little with varying GDL contact angle, and no water is present on the top wall, but the water coverage ratio on the GDL surface decreases significant as the contact angle increases. That is because a higher contact angle lifts the water up, i.e., move from a film to a droplet, which is more conducive to water removal from the GDL surface. For the hydrophobic GDL surface, the water coverage ratio on both the GDL and side/top surfaces changes little with varying the GDL contact angle, and a the stable flow pattern is formed as observed in the base case. It is worth to note that for 90o case, the water coverage ratio on the GDL surface is similar to that of hydrophilic GDL surface but the water coverage ratio on the side wall is close to that of hydrophobic GDL surface, indicating that the flow pattern for this case is in a transition state. The contact angle is big enough to lift up the water, leaving less water on the GDL surface, but the lift is not sufficient to make the liquid to reach the top wall. As a result, the water coverage ratio on the side wall is the highest among all the cases.  70  Figure 3.9: Effects of GDL surface wettabilities on time-averaged water distribution.  The time-averaged water volume fraction is also shown in Figure 3.9. For both hydrophilic and hydrophobic GDL surfaces, the water volume fraction decreases a little with an increase in contact angle. For a hydrophilic GDL surface, the liquid water is more in the form of a film. Increasing the GDL contact angle lifts up the water more from the GDL surface to occupy more cross sectional channel area, hence facilitating liquid water being flushed out of the channel. For a hydrophobic GDL surface, the flow pattern is similar to the base case, where increasing the contact angle helps the developing droplet to reach the top wall, which also facilitates the liquid water being flushed out of the channel. The corresponding time-averaged effects of GDL wettability on pressure drop is shown in Figure 3.10. As discussed previously, the pressure drop is mainly caused by the liquid blockage of the gas channel. Thus, increasing the GDL surface contact angle should always increase the pressure drop. For the hydrophobic GDL surface, however, the pressure drop varies little with increasing the contact angle, since the water coverage ratio on the GDL surface almost remains the same as shown in Figure 3.9.  71  Figure 3.10: Effects of GDL surface wettabilities on time-averaged pressure drop.  It can be concluded from the above observation that increasing the hydrophobicity of the GDL surface is helpful to expel liquid water from the GDL surface and also reduce the water fraction in the channel, although the pressure drop only increases slightly. The wettability of the GDL also affects the water transport inside the GDL. Thus, it is difficult to conclude whether high or low wettability is beneficial for water management in the whole fuel cell system without coupling the mass transfer and electrochemical reactions into the hydrodynamics. 3.2.4. Effects of channel wall surface contact angles Instead of changing the GDL surface wettability, one can also change the channel wall wettability for a given GDL to improve the water management in the gas channel. The effects of channel wall surface wettability were investigated by varying the water/wall surface static contact angle from 0o to 180o. Figure 3.11 shows the effects of different wall surface contact angles, i.e., 0o, 30o, 45o, 60o,  72  90o, 120o, and 180o, on the two-phase flow pattern. For hydrophilic wall surfaces, three-stage flow patterns similar to the base case are observed. More hydrophilic wall surfaces make the droplet on the GDL easier to move to the top wall, and also decrease the thickness of the liquid film on both side and top walls. For hydrophobic wall surfaces, the flow pattern is similar to that of the 1-pore case with the droplet flowing on the GDL surface only and more hydrophobic wall surface prevents the formation of big droplets.  (a)  (b)  (c)  73  (d)  (e)  (f)  (g) Figure 3.11: Effects of channel wall wettabilities on two-phase flow patterns:(a) θ=0o, (b) θ=30o, (c) θ=60o, (d) θ=90o,(e) θ=120o, (f) θ=140o, and (g) θ=180o.  The effects of channel wall surface wettability on water distribution are shown in Figure 3.12.  74  Figure 3.12: Effects of channel wall wettabilities on time-averaged water distribution.  As discussed previously, the more hydrophobic channel walls suppress liquid water removal from the GDL surface. As a result, the water coverage ratio on the top and side walls always decreases with increasing channel wall contact angle. Although the water volume fraction in the channel is higher for the more hydrophilic channel wall, the water coverage ratio on the GDL is smaller, resulting in more area for the gas reactants to diffuse through the GDL, leading to improved fuel cell performance. The pressure drop with different channel wall wettabilities (Figure 3.13) also suggests that a more hydrophilic channel wall should be used in fuel cells to reduce the energy loss.  75  Figure 3.13: Effects of channel wall wettabilities on time-averaged pressure drop.  3.2.5. Effects of liquid flow rates As specified in the boundary conditions section, the liquid flow rate used in the base case simulation, i.e. 7 g∙hr-1, is much higher than the theoretical values according to the reaction rates, e.g. 0.2 mg∙hr-1. The effects of liquid flow rates were thus investigated in order to make sure the previous simulation results reveal the flow patterns in a real PEM fuel cell. In this section, two cases with liquid flow rate reduced to 10% and 1% of base case flow rate are investigated. Simulation results shown in Figure 3.14 indicate the three-stage flow pattern remains the same even when the liquid flow rate is reduced by two orders of magnitude.  76  (a)  (b)  Figure 3.14: Effects of liquid flow rates on two-phase flow patterns: (a) 1/10th, (b) 1/100th.  Water distribution with respect to liquid flow rates in Figure 3.15 also shows that a lower liquid flow rate only results in a lower water coverage ratio on the top and side walls, but has little influence on the water coverage ratio on the GDL surface. This is because at the lower liquid flow rate, the liquid accumulation rate on the side walls is much smaller, and the adhesive force between top wall liquid thin film and side wall droplets becomes weaker, which makes the droplet detachment from the top wall much easier. As a result, the water coverage ratio on the top and side walls decreases rapidly due to the relatively rapid removal of liquid water from the top wall. Correspondingly, the time-averaged water volume fraction and time-averaged pressure drop (Figure 3.16) also become lower.  77  Figure 3.15: Effects of liquid flow rates on time-averaged water distribution.  Figure 3.16: Effects of liquid flow rates on time-averaged pressure drop.  78  3.3. Conclusions The following conclusions can be drawn from this simulation study: 1) For a microstructured GDL surface, three stages of two-phase flow patterns can be identified when water is injected through the GDL into the gas channel, namely emergence and merging of liquid water on the GDL surface, accumulation on the side walls, and detachment from the top wall. The flow patterns can be characterized by the channel crosssectional water distribution and coverage on the different walls. Water on the hydrophobic GDL surface tends to form droplets, while water on any of the hydrophilic surface tends to form a thin liquid film. The total pressure drop in the flow channel is mainly caused by droplet blockage of the channel by the droplet. 2) The microstructure of the GDL surface has a significant impact on the twophase flow patterns in the flow channel. The uniform inlet case and 1-pore case, approximations commonly used in previous studies, showed distinct flow patterns which are quite different from those observed in multiple pores. For the 4-pore, 16-pore and 64-pore cases, the flow patterns are similar, suggesting that the 4-pore case may be a minimum requirement to represent the microstructure of the GDL surface with a reasonable CPU time. 3) The wettability of both the GDL surface and channel walls also has a significant influence on the two-phase flow patterns in the flow channel. It was shown that more hydrophobic GDL surface and/or more hydrophilic channel walls would be helpful to remove the liquid droplets from the GDL surface to the channel walls. Decreasing liquid water coverage of the GDL will facilitate the gas reactants diffusion through the GDL to the catalyst layer and should improve the PEM fuel cell performance.  79  4) A lower liquid water flow rate into the channel can facilitate faster water removal due to the relatively faster detachment of water from the top wall. However, simulation results also showed that the liquid flow rate had little influence on the two-phase flow patterns formed in the channel. Thus, in order to shorten the computational time, the liquid flow rates could be amplified by several orders of magnitude in order to study the flow characteristics qualitatively in operating fuel cells.  80  4. Two-phase flow patterns in gas flow channel1 As reviewed in the chapter 1, the two-phase flow pattern in PEMFC gas flow channels has a great impact on the reactant distribution and the pressure drop. Several flow patterns have been identified in the PEM fuel cell channel both from in-situ [15, 98] and from ex-situ [94, 99] experiments, such as mist flow, droplet flow, film flow, annular flow and slug flow. Hussaini and Wang [98] and Lu et al. [99] have developed two-phase flow regime maps in terms of superficial gas and liquid velocities. These flow regime maps are useful for selecting optimal operating conditions. However, the two-phase flow pattern depends not only on the gas and water flow rate, but also the channel design and operating procedure, material properties, etc. The two-phase flow pattern may also vary along the length of the channel, since water is constantly introduced into the channels from the GDL. Therefore, there is still a lack of consistency on the two-phase flow patterns in PEM fuel cell channels in the literature. The main challenge in simulating flow patterns in fuel cells is how to consider the GDL microstructure. Most of the simulations do not address this issue, and have treated the GDL surface either as a homogenous surface or a surface with only a single open pore for liquid water injection. Some pore-scale models, such as the pore-network model and the Lattice Boltzmann model, are able to consider the detailed microstructure of the GDL, but at the micro-scale, which is thus not suitable for simulating the two-phase flow patterns at the macro-scale that form in the flow channel. In chapter 3, a simplified microstructure for the GDL surface was proposed in order to simulate the multiple droplet behaviour in a microchannel of a PEMFC. It was found that, in the channel width direction, at  1  A version of this chapter has been published. Y. Ding, H.T. Bi, D.P. Wilkinson. Threedimensional numerical simulation of water droplet emerging from a gas diffusion layer surface in micro-channels. J. Power Sources 2010; 195: 7278-7288. Part of results were presented in the conference: 8th International Conference on Nanochannels, Microchannels, and Minichannels, 2010, Montreal, Canada.  81  least 2 pores are required to represent the microstructure of the GDL surface, and the two-phase flow pattern in the channel only changes slightly with further increase in the pore number. In this chapter, the same approach as used in last chapter was extended to a much larger channel by implementing a simplified microstructure, which is comparable with realistic PEMFC gas flow channels. The effects of liquid injection rates and surface wettability on two-phase flow patterns in the flow channel were then investigated.  4.1. Computational domain and boundary conditions The three-dimensional computational domain is shown in Figure 4.1. A cuboid channel, which is used in all the simulations, has a 1.0 mm × 1.0 mm square cross section and is 100 mm in length with a hydrophobic GDL surface on the bottom and three hydrophilic channel walls. There is a 20 mm long entrance region before the water emergence area, which ensures the gas flow being fully developed before contacting the water droplets. In Chapter 3, it was shown that two pores in the width direction are enough to represent the microstructure of the GDL surface in a microchannel. Therefore, in this chapter, the same strategy is employed to simplify the GDL surface structure by opening 320 pores on the GDL surface with the same diameter of 400 µm. The selection of pore diameter ensures the total open area of GDL surface is around 50%, which corresponds to typical GDL surfaces with a porosity of 0.5. Air flows into the channel from one end and liquid water is injected from the multiple pores along the GDL surface.  82  (a)  (b)  (c)  Figure 4.1: Three-dimensional computational domain: (a) overview (b) simplified GDL microstructure (c) Local view of GDL microstructure with meshes.  For the base case, the velocity of air was set at 5 m∙s-1 at the atmospheric pressure, which is of the same order of magnitude as flows encountered in automotive fuel cell stacks [90]. The liquid injection velocity was set at 10-4 m∙s-1 for all the pores, which corresponds to the theoretical liquid generation rate at a current density of 0.5 A∙cm-2. Laminar flow and the non-slip boundary condition are assumed since the Reynolds number of each phase is quite small (Reg=458, Rel=11.9). The static contact angles of GDL surface and channel wall surface were set at 140o and 45o, respectively, based on typical PTFE treated carbon paper GDL materials and carbon plate. Air and water physical properties, i.e., densities, viscosities and surface tension coefficient were all set to a typical PEM fuel cell operating temperature of 70 oC. The model is isothermal, thus, liquid evaporation and condensation were not accounted. The time step for the baseline simulation was set at 10-6 s, which ensures that the global courant 83  number is less than 1. 149,690 meshes in total were used in the simulation with a mesh size of 0.1 mm. Figure 4.1(c) shows the local view of meshes at water inlet and GDL surface. Time step and mesh size independency were examined for the base case, to prove that they were small enough.  4.2. Results and discussion 4.2.1. Two-phase flow patterns in gas flow channels It is worth noting that since the length to height ratio of the channel is 100:1, it is difficult to have an overview of the whole channel. Therefore, for a better view, the length of the channel has been suppressed by a factor of 10 times, making the spherical droplets ellipsoid-like. Meanwhile, the undistorted two-phase flow patterns at certain location are also given by snapshots besides the overall view. As shown in Figure 4.2, for the base case, the three-stage droplet behaviour, which was found in Chapter 3, namely emergence and merging of liquid water on the GDL surface, accumulation on the side walls, and detachment from the top wall, could be identified in the channel. Liquid droplets first emerge from the pores and then coalesce on the GDL surface to form larger droplets. Due to the wettability difference between the GDL surface and channel walls, droplets on the GDL surface tend to attach to side walls (emergence and merging stage). Then, as the liquid water is constantly injected from the bottom, more and more water accumulates on the side wall, resulting in even larger droplets (accumulating stage). Due to the drag force exerted by the gas, large droplets begin to move slowly along the side wall, which forms so-called bottom corner droplet flow on the GDL surface. As these droplets move forward, more and more droplets coalesce and the droplets grow larger. Once a droplet hits the top wall, it rapidly spreads out on the top wall, and due to its much faster speed, it quickly detaches itself from the bottom corner, with some water on the side wall also being dragged away (detachment stage). As a result, corner droplet flow is formed. Since the droplet on the top wall moves much faster, it gives them more  84  opportunities to coalesce with other droplets that are still sitting on the bottom corner. Therefore, droplets flow on the top corner is usually larger than those on the bottom corner.  Figure 4.2: Two-phase flow patterns in the channel.  To quantitatively illustrate the two-phase flow pattern evolution along the channel, the time-averaged water coverage ratio on different surfaces along the channel was calculated as shown in Figure 4.3. The x-axis indicates the distance from the channel inlet, and the water coverage ratio is defined as the ratio of the surface area covered by water to the total surface area. The GDL water coverage ratio is  85  more important since the water on the GDL surface blocks the reactant diffusive pathway. To maintain a good cell performance, GDL water coverage should be kept as low as possible. In this case, the high GDL water coverage ratio means that the flow pattern is mainly bottom corner flow. The fluctuation in water coverage on the GDL and side wall reflects the droplet detachment from the bottom. The top wall water coverage ratio is still increasing, indicating that the flow is still developing along the channel.  Water Coverage Ratio (%)  100 80  Top Wall GDL Side Wall  60 40 20 0 0.00  0.02  0.04  0.06  0.08  0.10  X (m) Figure 4.3: Time-averaged water coverage ratio on different surfaces along the channel.  4.2.2. Effects of liquid flow rates In the base case, the liquid injection rate was set according to the theoretical liquid generation rate by the electrochemical reaction in the cathode catalyst layer using dry air as the reactant gas. Practically, the inlet gas is usually humidified to avoid membrane dehydration. Therefore, the liquid water formation rate in the cathode side channel is usually much higher than that generated by the reaction due to water condensation. Also, in an active PEM fuel cell the gas  86  channel is much longer, and more water tends to accumulate in the channel especially in the downstream section of the channel, resulting in various twophase flow patterns as observed in the literature. However, simulating such a long channel requires an extremely huge computational time. To shorten the water accumulation, we simply increased the liquid injection rates by 10 (x10 case), 100 (x100 case) and 1000 (x1000 case) times, respectively to mimic the flow patterns in the downstream section of a long channel. The two-phase flow pattern for each case is shown in Figure 4.4, and the timeaveraged water coverage ratio on different surfaces along the channel is shown in Figure 4.5. For the x10 case (Figure 4.4.a), droplets emerge into the channel much faster than that in the base case (Figure 4.2), resulting in more and bigger droplets on the GDL surface. Meanwhile, the detachment occurs more frequently, with more water flowing on the top wall. Since the top wall is hydrophilic, the droplets in the top corner tend to form a liquid film. Therefore, the flow pattern for the x10 case is bottom corner droplet flow and top wall film flow. Figure 4.5.a shows that the GDL surface water coverage ratio reaches its highest point at about 0.04 m, indicating that beyond this point, droplets detach faster than they emerge. The top wall water coverage ratio begins to exceed the GDL surface water coverage ratio at about 0.065 m, and it still continues to increase at the end of channel, indicating that the flow is still developing.  87  (a)  (b)  (c)  Figure 4.4: Effects of liquid flow rate on the two-phase flow patterns in the channel: (a) x10 case, (b) x100 case, and (c) x1000 case.  88  Water Coverage Ratio (%)  100 80  Top Wall GDL Side Wall  60 40 20 0 0.00  0.02  0.04  0.06  0.08  0.10  0.08  0.10  X (m) (a)  Water Coverage Ratio (%)  100 80  Top Wall GDL Side Wall  60 40 20 0 0.00  0.02  0.04  0.06  X (m) (b)  Water Coverage Ratio (%)  100 80 60 40  Top Wall GDL Side Wall  20 0 0.00  0.02  0.04  0.06  0.08  0.10  X (m) (c)  Figure 4.5: Time-averaged water coverage ratio on different surfaces along the channel: (a) x10 case, (b) x100 case, and (c) x1000 case.  89  Further increasing the liquid flow, as shown in the x100 case (Figure 4.4.b), results in more water present on all the channel walls. Droplets on the GDL surface merge together and the liquid film on the top wall covers almost all the top wall surface, which makes the flow pattern in the channel similar to annular flow. The water coverage ratio on each wall eventually becomes stable as shown in Figure 4.5.b. Most of the water flows on the top wall, which results in the GDL water coverage ratio being similar to the x10 case. When the liquid flow rate is amplified by 1000 times, quite distinct flow patterns are formed in the channel (Figure 4.4.c). Liquid water covers almost all the channel wall surface, and becomes a continuous phase. An annular flow pattern develops in the entrance section, and then gas bubble and slug flow are identified in the exit section of the channel. Due to the density differences between water and air, discrete gas bubbles can easily escape from the top of the channel. Therefore, the GDL surface water coverage ratio is much higher than that of the top wall (Figure 4.5.c), which indicates a much different water distribution for very high water injection rate. These selected four cases imply that, the flow patterns in a longer PEM fuel cell flow channel may follow a flow pattern evolution along the length of channel, namely, corner droplet flow at the beginning of channel, followed by top wall film flow, annular flow, and finally slug flow in the channel. From a practical point of view, slug flow should be avoided as much as possible, since it causes an extremely high GDL surface water coverage ratio and therefore significant loss in fuel cell performance. The effects of liquid flow rates on the time-averaged water volume fraction (also called water saturation) and total pressure drop were also analyzed as shown in Figure 4.6. The water volume fraction or the water saturation indicates the degree of channel flooding. The pressure drop, which is a key parameter, indicates the energy loss for the fluid flowing through the channel. It is obvious  90  that the pressure drop increases exponentially with the increasing water volume fraction. Therefore, water accumulation in the channel for an operating PEM fuel cell should be avoided as much as possible to reduce parasitic energy loss, reactant flow maldistribution, and poor cell performance.  Figure 4.6: Effects liquid flow rate on the pressure drop and water volume fraction in the channel.  4.2.3. Effects of GDL surface wettability The wettability of the GDL, which is characterized by the surface contact angle, plays a significant role in PEM fuel cell water management. To investigate the impact of GDL surface wettability on the two-phase flow patterns formed in the channel, several GDL surface contact angles, i.e., 45o, 60o, 90o, 120o, and 140o were tested in the simulation. The contact angle of the other walls was kept constant at 45o. The liquid water injection rate was amplified by 10 times as in the x10 case.  91  Figure 4.7 shows quite different two-phase flow patterns in the channel for different GDL contact angles. When the GDL surface contact angle is hydrophilic, i.e., Figure 4.7a, b, c, water on the GDL surface tends to form a liquid film. The lower the GDL contact angle, the higher the GDL area that is covered by the water, which blocks the reactants pathway to the catalyst layer and leads to decreased fuel cell performance. There is no water flowing on the top wall for these hydrophilic cases. Therefore, the flow pattern in the channel is bottom film flow or bottom corner flow. As the contact angle of the GDL surface increases, liquid water begins to move to and accumulate on the side walls and the higher the contact angle, the more water moving from the GDL surface to the side walls. When the GDL surface is hydrophobic, i.e. Figure 4.7 d, e, water begins to form liquid droplets rather than liquid films, and the three-stage droplet behaviour can be identified again. The two-phase flow pattern becomes corner droplet flow as observed in the base case.  (a)  (b)  92  (c)  (d)  (e)  Figure 4.7: Effects of contact angle of GDL surface on the two-phase flow patterns in the channel :(a) θ=45o, (b) θ=60o,(c) θ=90o, (d) θ=120o, (e) θ=140o.  The effects of GDL wettability on the time-averaged water distribution are shown in Figure 4.8. For the hydrophilic GDL surface, the water coverage ratio on the side walls changes very little with varying GDL contact angle due to the balance of liquid droplet moving from the GDL to the side wall and faster expelling out of the channel. No water is present on the top wall, and the water coverage ratio on the GDL surface decreases significantly as the contact angle increases. This is because a higher contact angle lifts the water up, changing the water on the GDL  93  from a liquid film to droplets, thus occupying less GDL surface area. On the other hand, for the hydrophobic GDL surface, further increase in contact angle decreases the water coverage ratio on the GDL surface. Meanwhile, the liquid droplets are able to touch the top wall, resulting in the top corner flow, which decreases the water coverage ratio on the side walls slightly.  Figure 4.8: Effects of GDL surface contact angle on the water distribution in the channel.  The effects of GDL surface wettability on the time-averaged pressure drop and water volume faction are shown in Figure 4.9. Increasing the GDL surface contact angle always increases the total pressure drop, which is consistent with experimental results [143, 144]. This is because higher hydrophobicity of the GDL surface lifts more water up from the GDL surface to occupy more of the cross sectional channel area, which blocks the gas pathway in the channel, resulting in higher pressure drop. In addition, the increased drag forces between the liquid and gas also shorten the water residence time in the channel and thus reduce the water volume fraction.  94  Figure 4.9: Effects of GDL surface contact angle on the pressure drop and water volume fraction in the channel.  It can be concluded from the previous observations that increasing the hydrophobicity of the GDL surface is helpful to expel liquid water from the GDL surface and reduce the water volume fraction in the channel, but the pressure drop also increases slightly. Therefore, in the selection of the GDL contact angles both avoiding flooding and reducing parasitic energy loss should be considered. 4.2.4. Effects of channel wall surface wettability The effects of channel wall surface wettability were also investigated by varying the contact angle of channel walls from 45o to 140o, i.e., 45o, 60o, 90o, 120o, and 140o. The contact angle of the GDL surface was kept at 140o for all these cases. The liquid water injection rate was amplified by 10 times as in the x10 case. Figure 4.10 shows the two-phase flow patterns formed in the channel for different channel wall contact angles. For hydrophilic channel walls (Figure 4.10 a, b, c), the flow pattern is corner droplet flow as observed in the base case. For more  95  hydrophilic wall surfaces, the droplets on the GDL surface move more easily to the top wall, increasing the droplet detachment frequency and forming larger droplets on the top corner. Thus, it facilitates faster water removal from the channel. For hydrophobic channel walls (Figure 4.10 d, e), the flow pattern is droplet flow on the GDL surface, and no water flows on the top wall. All the droplets are expelled from the hydrophobic surfaces, occupying a large cross sectional area of the channel, and thus become easier to be flushed out by the gas. It is also observed that the highly hydrophobic channel walls can further prevent the formation of larger droplets in the channel and reduce the droplet resident time.  (a)  (b)  96  (c)  (d)  (e) Figure 4.10: Effects of contact angle of channel side and top wall surfaces on the two-phase flow patterns in the channel: (a) θ=45o, (b) θ=60o, (c) θ=90o, (d) θ=120o, (e) θ=140o.  The effects of channel wall wettability on the water distribution in the channel are shown in Figure 4.11. It is found that when the channel wall is hydrophilic, increasing the channel wall contact angle significantly decreases the water coverage ratio on the top and side walls, due to the prevention of droplet formation on the side walls. Correspondingly, the GDL surface coverage also decreases slightly, since the larger droplets formed on the GDL surface accelerate their removal from the channel. When the channel wall is hydrophobic,  97  stable droplets are formed, with little water attached onto the top and side walls. Therefore, any further increase in the channel wall contact angle has little impact on the water distribution in the channel.  Figure 4.11: Effects of channel wall contact angle on the water distribution in the channel.  The effects of channel wall wettability on the pressure drop and water volume fraction are shown in Figure 4.12. The pressure drop is seen to increase first with increasing the channel wall contact angle due to the larger water cross sectional area, followed by a slight decrease due to the smaller droplets formed on the GDL surface. Increasing the channel wall contact angle always decreases the water volume fraction in the channel, because of the reduced water resident time. The previous results indicate that using a more hydrophilic channel wall is helpful to reduce the pressure drop but the GDL water coverage ratio also increases slightly. The selection of optimal channel wall wettability thus depends on the specific requirement of fuel cell applications.  98  Figure 4.12: Effects of channel wall contact angle on the pressure drop and water volume fraction in the channel.  4.3. Conclusions The following conclusions can be drawn from this chapter: 1) Three stages of droplet behaviour can be identified, namely, i) emergence and merging of liquid water on the GDL surface, ii) accumulation on the side walls, and iii) detachment from the top wall, For the base case, i.e., at theoretical liquid water production rate, the flow pattern is corner droplet flow, with liquid water mainly presents in the bottom corners. 2) With increasing liquid injection rates, the flow pattern evolves from corner droplet flow to top wall film flow, annular flow, and finally slug flow. 3) To reduce the parasitic energy loss, the water volume faction should be kept as low as possible in the channel, since the pressure drop increases exponentially with it.  99  4) The material wettability has a great impact on the two-phase flow pattern, water distribution and pressure drop. Using a more hydrophobic GDL surface is helpful to expel water from the GDL surface, but also increases the pressure drop. On the other hand, using a more hydrophilic channel wall reduces the pressure drop, but increases the GDL water coverage ratio slightly and the water residence time.  100  5. Impacts of two-phase flow maldistribution on flow behaviours in parallel gas flow channels Flow maldistribution always happens in PEM fuel cells when using common inlet and exit headers to supply reactant gases to multiple channels. This maldistribution is problematic because it leaves the more flooded channels without sufficient air flow and the less flooded channels with excessive air flow. As a result, it lowers fuel cell performance, leads to the loss of cell efficiency, and causes durability problems. As reviewed in Chapter 1, possible reasons for flow maldistribution include different flow resistances caused by inherent geometric differences (e.g., manufacturing tolerances, channel defect or corrosion), GDL intrusion during compression [114], uneven liquid water accumulation in the flow field channels [115], and two-phase flow instabilities [126]. Several analytical models [121, 145-148] and CFD-based simulations [114, 120, 122, 148] have been utilized to study the flow maldistribution in multiple channels. However, most of those works did not consider the existence of liquid water, except that Basu et al. [145] who used their two-phase flow analytical model and found that flow maldistribution is always more severe in a two-phase flow than in a singlephase flow field. They also used the mixture model to study the mitigation of flow maldistribution [114, 122]. However, the mixture model is not capable of capturing the liquid droplet or the slug flow which was found to be mainly responsible for the severe maldistribution in parallel channels from experiments [99]. Thus, the VOF method is probably a more suitable two-phase flow model for investigating the flow maldistribution phenomena in multiple channels. In this chapter, the impact of gas flow maldistribution on two-phase flow pattern, pressure drop and water coverage in two parallel channels is simulated using VOF method. The gas maldistribution is introduced by injecting water into the channel with different flow rates. Effects of liquid injection rates and material wettability on the degree of maldistribution are also discussed. Finally, the  101  potential impact of the GDL microstructure on flow maldistribution in parallel channels is investigated.  5.1. Computational domain and boundary conditions Figure 5.1 shows the three-dimensional computational domain of the two parallel channels, which have a 1.0 mm × 1.0 mm square cross section for each branch and a 100 mm in length. The distance between two parallel channels is 1 mm. Air flows into the channel from the left end and liquid water is injected from two pores at the entrance of each branch on the GDL surface. The setting of twopore-in-the-width-direction ensures that the two-phase flow patterns can be reproduced as discussed in Chapter 3. The distance from the manifold to the water inlets is 5 mm.  Figure 5.1: Three-dimensional computation domain of two-parallel channel.  For all the cases, the velocity of air is set at 4.464 m/s at atmospheric pressure and, for the base case, the water injection velocity is set at 0.1 m/s for all the  102  pores. The operating conditions are set according to a typical operating PEM fuel cell (cathode side) at a current density of 1 A/cm2. Liquid flow rates in each branch are set differently to mimic the severe liquid maldistribution caused by non-uniform water generation. The liquid velocity is amplified by 100 times to shorten the computation time. Same as in Chapter 4, the GDL surface, i.e., the bottom surface is hydrophobic with a contact angle of 140o, and all the other surfaces are hydrophilic with a contact angle of 45o. The total meshes of the geometry are about 220,000 and the time step for all the simulations is less than 10-6 s to ensure the global courant number below 1.  5.2. Results and discussion 5.2.1. Effects of severe flow maldistibution To mimic the severe liquid maldistribution, liquid velocity in the left branch is set at 0.01 m/s, one-tenth of that in the base case. Liquid velocity in right branch is increased to 0.19 m/s to keep the total liquid flow rate the same as in the base case. Figure 5.2 shows the two-phase flow pattern in the two parallel channels with and without maldistribution. The blue color indicates liquid droplets. For the base case (Figure 5.2a) without maldistribution, the flow pattern is the top-wall droplet flow, which is consistent with the case in Chapter 4, and this flow pattern is formed mainly due to the wettability difference between the GDL surface and other channel walls. When more water enters the right branch (Figure 5.2b), larger droplets are formed due to the higher liquid flow rate and the lower gas flow rate. As they flow along the channel, some droplets coalescence and then elongate on the top wall. Especially near the right branch outlet, slug flow sometimes appears. In the left branch, much smaller droplets are formed and quickly removed from the channel due to a high gas flow rate.  103  (a) (b) Figure 5.2: Effects of maldistribution on two-phase flow patterns in parallel channels: (a) w/o maldistribution, (b) with maldistribution  The average gas velocities in each branch are shown in Figure 5.3. The error bar gives the standard deviation of gas velocity after the steady state is reached. Apparently, flow maldistribution not only alters the gas distribution in two branches, but also results in more gas flow fluctuation, which may be detrimental to the fuel cell stable operation. This fluctuation is associated with the formation of slug flow in the right branch. When slug flow forms, it causes a high flow resistance, thus reduces the gas flow in the branch.  104  5  Gas Velocity (m/s)  4  Left Branch Right Branch  3 2 1 0  w/o mal.  with mal.  Cases Figure 5.3: Effects of maldistribution on averaged gas velocity in each branch.  Figure 5.4 shows the average total pressure drop and liquid volume fraction in each branch. As discussed above, the existence of slug flow increases the total pressure drop and its fluctuations. The unstable gas flow in each branch also influences the liquid removal, therefore increases the fluctuation of liquid volume fraction. The total pressure drop with a liquid maldistribution is higher than the base case with uniform distribution, which accelerates the liquid being expelled from the channel. Thus the total volume fraction with maldistribution (9.2%) is smaller than the base case (11.8%).  105  0.30 Pressure Drop Left Branch Right Branch  600  0.25 0.20 0.15  400  0.10 200 0.05 0  w/o mal.  with mal.  Liquid Volume Fraction (-)  Pressure Drop (Pa)  800  0.00  Cases Figure 5.4: Effects of maldistribution on total pressure drop and liquid volume fractions in each branch.  Figure 5.5 shows the effects of maldistribution on wall liquid coverage in each branch. The liquid coverage ratio is defined as the ratio of the surface area covered by liquid to the total surface area. The GDL water coverage ratio is a key parameter to indicate the impact of liquid water on PEM fuel cell performance; because GDL surface coverage by liquid water reduces the reactant gas diffusion into the catalyst layer, where the electrochemical reaction occurs. The total top and side wall liquid coverage ratios with maldistribution are much smaller than the base case. One reason is that in the base case, droplets are smaller and thinner, which are easy to spread on the channel surfaces. Meanwhile, more total liquid volume fraction in the base case further increases the liquid coverage ratio on the top and side walls. Since slug flow formed in the channel with a maldistribution, the total GDL liquid coverage ratio is a little bit higher, which may reduce the fuel cell performance if reactions are included in the simulation.  106  Liquid Coverage Ratio (-)  0.5 Left GDL surface Right GDL surface Left Top Wall Right Top Wall Left Side Walls Right Side Walls  0.4 0.3 0.2 0.1 0.0  w/o mal.  with mal.  Cases Figure 5.5: Effects of liquid maldistribution on wall liquid coverage in each branch.  5.2.2. Effects of the degree of flow maldistribution In operating PEM fuel cells, different degree of flow maldistribution may occur in multiple channels. Some channels may totally be blocked by liquid water, causing more gas flows through other channels, which may lead to membrane dehydration there. To study the effects of the degree of flow maldistribution on fuel cell performance, different liquid injection rate ratios between the two parallel channels were set by changing the water flow rate in one branch while keeping the total water flow rate the same. Figure 5.6 shows the two-phase flow patterns under different liquid injection rate ratios.  107  (a)  (b)  (c)  (d)  (e) (f) Figure 5.6: Effects of the degree of flow maldistribution on two-phase flow patterns.  108  QL and QR indicate the liquid injection rates in the left and right branch, respectively. In general, the two-phase flow pattern mostly remains in the topwall droplet flow under different degrees of flow maldistribution. Fewer droplets present in the left branch and their sizes become smaller when the severity of maldistribution is increased. Meanwhile, in the right branch, droplets first become larger due to their coalescence (Figure 5.6 e,f), then because of less gas flowing through the channel, droplets become longer (Figure 5.6 a,b,c,d) and eventually transform into slugs in the channel. Figure 5.7(a) clearly shows the two tendencies of flow pattern development. When QL/ QR is higher than 0.25, the severity of gas maldistribution changes gradually mainly associated with the changing of liquid volume in the channel. When the liquid flow ratio is lower than 0.25, slug flow begins to appear, and gas maldistibution changes more substantially with the liquid ratio. The gas flow rate ratio‟s variation with liquid flow rate ratio, as shown in Figure 5.7(b), has a hyperbola shape of curve, which is consistent with experimental results [149]. The transformation of two-phase flow pattern well explains the variation of gas flow rate ratio with liquid flow rate ratio. 9 8  Left Branch Right Branch  Gas Flow Rate Ratio (-)  Gas Velocity (m/s)  4  3  2  1  7 6 5 4 3 2 1  0  0 0.0  0.2  0.4  0.6  0.8  Liquid Flow Rate Ratio (-)  1.0  0.0  0.2  0.4  0.6  0.8  1.0  Liquid Flow Rate Ratio (-)  (a) (b) Figure 5.7: Effects of the degree of flow maldistribution on (a) gas velocity (b) gas flow rate ratio.  109  The GDL surface liquid coverage ratios in each branch are shown in Figure 5.8. Although the slug flow does not exist in the left branch, there still is some liquid (about 0.3%) covering the GDL surface due to the emerging droplets from the GDL surface, as shown in Chapter 3, Increasing the gas velocity will elongate the droplets and thus increase the GDL liquid coverage. Therefore, GDL liquid coverage ratio in the left branch firstly increases with increasing the severity of maldistribution (QL/ QR from 1.0 to 0.4), and the opposite trend is followed in the right branch. Increasing gas flow rate also advances the droplet detachment. Therefore, further increasing gas flow rate will decrease the GDL coverage, as observed in the left branch (QL/ QR from 0.4 to 0). The sudden increase (QL/ QR from 0.25 to 0.0) in the right branch is due to the transition to slug flow and the sudden drop (QL/ QR at 0.0) in the left branch is caused by the absence of liquid  GDL Surface Liquid Coverage Ratio (-)  in the channel.  0.008  0.006  Left Branch Righ Branch  0.004  0.002  0.000  0.0  0.2  0.4  0.6  0.8  1.0  Liquid Flow Rate Ratio (-) Figure 5.8: Effects of the degree of flow maldistribution on GDL surface liquid coverage ratio.  The total pressure drop at different degrees of flow maldistribution (Figure 5.9) is complicated. The pressure drop in parallel channels depends on many factors,  110  e.g. two-phase flow patterns, droplet size, liquid distribution and flow instabilities etc. Zhang et al. [118, 126, 150] also recorded such a wavy shape of pressure drop from their experiments, and attributed it to the flow instabilities. But in general the pressure drop seems to be insensitive to the flow maldistribution, which are consistent with the existence of gas and liquid maldistribution in parallel two-phase flow channels.  Pressure Drop (Pa)  450  425  400  375  350  0.0  0.2  0.4  0.6  0.8  1.0  Liquid Flow Rate Ratio (-) Figure 5.9: Effects of the degree of flow maldistribution on pressure drop.  As shown in Figure 5.10, liquid volume fraction in the left branch continuously decreases (QL/ QR from 1.0 to 0) due to the presence of less liquid and more gas. The flow in the right branch is just in contrary, except that the two-phase flow pattern transforms from the droplet flow to slug flow.  111  Liquid Volume Fraction (-)  0.20  0.15  0.10  0.05  0.00  Left Branch Right Branch 0.0  0.2  0.4  0.6  0.8  1.0  Liquid Flow Rate Ratio (-) Figure 5.10: Effects of the degree of flow maldistribution on liquid volume fraction.  5.2.3. Effects of GDL surface hydrophobicity on flow maldistribution In Chapters 3 and 4, the GDL and channel wall hydrophobicity are found to have great impact on droplet behaviours and two-phase flow patterns. However, there is still no open literature on their effects on flow maldistribution. In this section, the effect of GDL surface hydrophobicity on flow maldistribution is investigated based on simulations by varying the GDL surface static contact angle from 45o to 120o, i.e., 45o, 60o, 90o, and 120o. In all those simulated cases, the liquid injection velocities in each branch are set at 0.01 m/s and 0.19 m/s, respectively. The obtained two-phase flow patterns under different GDL surface hydrophobicity are shown in Figure 5.11. As expected, slug flow appears more frequently in the right branch for more hydrophilic GDL surface, e.g. Figure 5.11a, b. Meanwhile, the total liquid volume holdup in the channels also decreases with increasing the hydrophilicity of GDL surfaces.  112  (a)  (b)  (c) (d) Figure 5.11: Effects of GDL surface hydrophobicity on the two-phase flow patterns :(a) θ=45o, (b) θ=60o, (c) θ=90o, (d) θ=120o.  Liquid slugs increase the local flow resistance, meanwhile, accelerate their removal and thus reduce the liquid volume holdup in the channel (Figure 5.12 b). These two mechanisms have opposite effects on the pressure drop, thus resulting in a parabolic shape of the pressure drop curve (Figure 5.12 a). Since there is no slug flow in the left branch, the liquid volume fraction there is dominated by the pressure drop, i.e., higher pressure drop corresponds to less liquid holdup. For hydrophobic GDL surface, e.g., with a contact angle of 120o and 140o, there is almost no slug flow in the right branch. Increasing the hydrophobicity of GDL would decrease the critical droplet detachment diameter,  113  leading to smaller droplets which are easier to be removed from the channel. Therefore, there is a slightly lower liquid content in the right branch when more hydrophobic GDL is used. This conclusion is consistent with previous ones drawn in Chapters 3 and 4. 0.20  Liquid Volume Fraction (-)  Pressure Drop (Pa)  450 430 410 390 370 350 30  60  90  120  150  0.15  0.10  Left Branch Right Branch  0.05  0.00 30  o  GDL Contact Angle( )  60  90  120  150  o  GDL Contact Angle( )  (a) (b) Figure 5.12: Effects of GDL surface hydrophobicity on (a) pressure drop (b) liquid volume fraction.  Using more hydrophilic GDL surface always increases the GDL liquid coverage ratio (Figure 5.13), which is not favored for the fuel cell performance. The flow maldistribution also becomes a little bit more severe when more hydrophilic GDL surface is used (Figure 5.14). Therefore, it can be concluded that more hydrophobic GDL surface should be used in PEM fuel cells to reduce flow maldistribution and increasing fuel cell performance.  114  GDL Liquid Coverage Ratio (-)  0.20  0.15  Left Branch Right Branch  0.10  0.05  0.00 30  60  90  120 o  150  GDL Contact Angle( ) Figure 5.13: Effects of GDL surface hydrophobicity on GDL liquid coverage ratio.  5  Gas Velocity (m/s)  4 3  Left Branch Right Branch  2 1 0 30  60  90  120 o  150  GDL Contact Angle( ) Figure 5.14: Effects of GDL surface hydrophobicity on flow maldistribution.  115  5.2.4. Effects of channel wall hydrophobicity on flow maldistribution In this section, the effects of channel wall hydrophobicity on flow maldistribution are simulated by varying the channel wall contact angle, i.e., 60o, 90o, 120o, and 140o. The settings of injected liquid velocities in each branch remain the same as in the prior section. Figure 5.15 shows that as the channel wall contact angle increases, droplets in the right branch become smaller and fewer. For the hydrophobic channel wall, i.e., Figure 5.15 c, d, slug flow is observed. Changing channel wall wettability also changes the droplet size and shape in the left branch, although its impact is insignificant because of low liquid flow rate.  (a)  (b)  (c) (d) Figure 5.15: Effects of channel wall hydrophobicity on the two-phase flow patterns :(a) θ=60o, (b) θ=90o, (c) θ=120o, (d) θ=140o.  116  Altering channel wall contact angle can exert a great impact on the pressure drop (Figure 5.16a). A decrease of pressure drop as the contact angle varies from 45o to 60o mainly results from the less resident time of liquid in the channel. The sudden pressure drop increase as the contact angle changes from 60o to 120o is related to the more frequent appearance of slugs. Further increasing the channel wall contact angle does not change the slug flow pattern but accelerates the droplet removal, thus decreases the pressure drop. Liquid volume fractions in both branches always decrease with increasing the hydrophobicity of the channel wall (Figure 5.16b), which is consistent with previous simulation results presented in Chapters 3 and 4.  0.20  Liquid Volume Fraction (-)  Pressure Drop (Pa)  600  500  400  300 30  60  90  120  150 o  Channel Wall Contact Angle( )  0.15  Left Branch Right Branch 0.10  0.05  0.00 30  60  90  120  o  150  Channel Wall Contact Angle( )  (a) (b) Figure 5.16: Effects of channel wall hydrophobicity on (a) pressure drop (b) liquid volume fraction.  Using more hydrophilic channel wall is an effective method to reduce the GDL liquid coverage (Figure 5.17), leading to improved fuel cell performance. The high GDL liquid coverage for hydrophobic channel walls indicates the existence of slug flow. In the left branch, however, more hydrophobic channel walls lower the GDL liquid coverage due to the fast droplet detachment.  117  GDL Liquid Coverage Ratio (-)  0.04  0.03  Left Branch Right Branch  0.02  0.01  0.00 30  60  90  120  o  150  Channel Wall Contact Angle( ) Figure 5.17: Effects of channel wall hydrophobicity on GDL liquid coverage ratio.  Channel walls with a contact angle close to 90o give the least flow maldistribution (Figure 5.18). Taking all factors into consideration, it is suggested that a moderate channel wall contact angle, e.g. 80o to 90o, should be used in PEM fuel cells. Within this range of contact angles, pressure drops and the GDL liquid coverage are much lower than in channels with hydrophobic channel walls. The liquid volume faction is also lower and the flow maldistribution is less severe than in channels with more hydrophilic channel walls.  118  5  Gas Velocity (m/s)  4 3 2  Left Branch Right Branch  1 0 30  60  90  120  o  150  Channel Wall Contact Angle( ) Figure 5.18: Effects of channel wall hydrophobicity on flow maldistribution.  5.2.5. Effects of liquid inlet position on flow maldistribution The GDL has a random microstructure as shown in Chapter 1. Therefore, water droplet may emerge from several pores that are randomly located in parallel channels. The random emergence of droplets may cause different flow resistance between channels, leading to flow maldistribution, even though the liquid flow rates in each branch are the same. To examine the effects of liquid inlet position on flow maldistribution, a new parallel channel configuration is considered as shown in Figure 5.19 b where water inlet in the left branch is located to the middle of the channel. The injected liquid velocity is set at 0.1 m/s for both cases.  119  (a) Base case  (b) Mid-inlet Figure 5.19: Comparison of water inlet positions of two cases.  A snapshot of the two-phase flow pattern for the mid-inlet case is shown in Figure 5.20. For both branches, the flow pattern remains in the top wall droplet flow. Since the pressure drop over the two branches must be equal, pressure gradient for each droplet in the left branch must be higher than that in the right branch. Therefore, droplets in the left branch are stretched and thus longer than the droplet in the right branch.  Figure 5.20: Simulated two-phase flow pattern for the mid-inlet case.  120  Figure 5.21 shows that the liquid inlet position has a huge impact on flow maldistribution, and it also causes large gas velocity fluctuations. The gas flow rate ratio (left branch to right branch) is about 1.37, almost the same degree as the ratio at a liquid flow ratio (QL/QR) of 0.67 in the base case. This result implies that the GDL microstructure is probably another factor that may influence the flow maldistribution.  5  Gas Velocity (m/s)  4  Left Branch Right Branch  3 2 1 0  Base Case  Mid-inlet  Cases Figure 5.21: Effects of liquid inlet position on the average gas velocity in each branch.  5.3. Conclusions The following conclusions can be drawn from this chapter‟s simulation study: 1) Severe flow maldistribution caused by uneven liquid flow rates in two parallel channels may transform the two-phase flow patterns from the top wall droplet flow to the slug flow in the higher-liquid-flow-rate-branch. It also leads to severe uneven gas distribution, and increases the pressure drop, gas velocity  121  fluctuations and the GDL surface liquid coverage, which are detrimental to the PEM fuel cell performance. 2) Varying the degree of flow maldistribution does not change the two-phase flow patterns except under severe maldistribution conditions, i.e. at a liquid flow rate ratio less than 0.25, where the slug flow starts to appear, gas maldistibution becomes more severe, and GDL surface liquid coverage soars up. Pressure drops at different degrees of flow maldistribution show a wavy shape that may be associated with the flow instabilities in parallel channels. 3) Using more hydrophobic GDL is helpful to avoid slug flow, reduce the total pressure drop, decrease the GDL surface liquid coverage and improve the gas maldistribution. Meanwhile, a moderate channel wall contact angle, e.g. 80o to 90o, should be used in PEM fuel cells because, within this range, pressure drops and the GDL liquid coverage are lower than in a channel with hydrophobic walls, the liquid volume faction is lower and the flow maldistribution is less severe. 4) Simulation results show that the GDL microstructure is probably another factor that may cause flow maldistribution.  122  6. Communicating  channels  for  flow  mal-  distribution mitigation1 Flow maldistribution in multiple channels of PEM fuel cells is almost unavoidable due to the different flow resistances caused by inherent geometric differences, asymmetrical channel design, two-phase flow instabilities or random droplet emerging. Several researchers have tried different designs to reduce flow maldistribution based on either simulations or experiments. Johnson et al., found that increasing the pressure differences between adjacent gas channels could improve water management or reactant distribution [151]. Based on this principle, several novel channel configurations have been designed, such as channels with varying depth or width [151], wave or sinusoidal shaped channels [152], and crossly connected channels [152]. Basu et al. [122] designed a flow splitter in the manifold, which was found to reduce the flow maldistribution by more than 50% based on simulation. Liu et al. [123] designed multiple bifurcations to distribute gas flow. Wang et al. [124] simulated the gas flow in a novel flow field of a bionic flow slab. Comparing to multiple parallel channels, it increased flow uniformity while maintaining a reasonable pressure drop. Hu et al. [125] designed a slottedinterdigitated channel, which was then optimized to reduce maldistribution. However, none of those works have considered the presence of liquid slugs in the channel or even two-phase flow in their designs. In parallel channels, the different flow resistance also creates a pressure difference between two adjacent channels. Since the GDL is porous media, it allows gas or liquid flows through it, i.e., flow communications. The existence of flow communications may produce two distinct impacts on flow distributions. When liquid flow through the GDL, it reduces the flow maldistribution by evening  1  A version of this chapter will be submitted for publication: Y. Ding, R. Anderson, L.F. Zhang, H.T. Bi, D.P. Wilkinson. Simulations of Two-phase Flow Distribution in Communicating Parallel Channels for a PEM Fuel Cell  123  out the liquid distribution, but if gas flows through the GDL from the high liquid flow channel to the low liquid flow channel, it may make flow maldistribution even worse. However, little attention has been paid to the effect of communication between parallel channels on flow distributions in the open literature. In view of the knowledge gap on the impact of flow communications on flow distributions in interconnected parallel channels for fuel cells, an attempt has been made to simulate two-phase flow in parallel channels with communications in this chapter. The investigated interconnected parallel channels resemble the gas and liquid communication in fuel cell flow fields due to inherent or artificial structures of GDLs. To investigate the effect of communicating channels, two parallel channels are connected at several locations along the length of the channels. The proposed benefits of the across-rib (or landing area) transport have been noted numerically in the study of adverse flow fields [153], but the projected benefits may be unrealistic or lessened since liquid water was not included in the simulation.  6.1. Computational domain and boundary conditions Figure 5.2(a) shows the schematics of simulated parallel communicating channels. Each parallel channel has a 1.0 mm × 1.0 mm square cross section and is 100 mm in length. The distance between the two parallel channels is 1 mm. This geometry has the same dimension as in Chapter 5 except that the two parallel channels are connected by several communication channels along the length of the channels. In order to investigate the effect of communication channels, five different channel designs are tested. There is no communication channel for the base case, and the other four cases have different width of communication channels, i.e., 0.25 mm, 0.5mm, 1 mm and 2 mm, respectively. The depth of the communication channel in those cases is 1 mm. To maintain the overall communication area being the same, the number of communication channels also varies with their width. There are 20 communication channels for the 0.25 mm case, 10 for the 0.5 mm case and 5 for the 1 mm case. For 2 mm 124  case, 5 communication channels are used to test the effects of total communication area. For all cases, communication channels are equally distributed along the two parallel communicating channels.  (a)  (b)  Figure 6.1: Three-dimensional computational domain (a) illustration of parallel communicating channels (b) base case  As shown in Figure 5.2(b), air flows into the channel from one end and liquid water is injected from two pores at the entrance of each branch on the GDL surface. For all the cases, the velocity of air is set at 4.464 m/s at atmospheric pressure, and the water injection velocity is set at 0.01 m/s in the left branch and 0.19 m/s in the right branch. The operating conditions are set according to an operating PEM fuel cell (cathode side) at a current density of 1 A/cm2. Liquid flow  125  rates in each branch are set differently to mimic the severe liquid maldistribution caused by non-uniform water generation. The liquid velocity is amplified by 100 times to shorten the computation time, since the two-phase flow pattern in the channel has been found not to change when the liquid flow rate is amplified several orders of magnitude [154]. The static contact angles of the GDL surface and channel wall surface are set at 140o and 45o, respectively, according to a real PEM fuel cell [155].  6.2. Results and discussion 6.2.1. Effects of communication channels on flow maldistribution Figure 6.2 shows the two-phase flow patterns for five cases simulating the operating conditions for the cathode side of a PEM fuel cell at 1 A/cm2. For the base case, as discussed in Chapter 5, the flow pattern in both branches is mainly top wall droplet flow where the liquid droplets are attached to the top wall. The slug flow also appears in the right branch sometimes. For the 2 mm communication channel case, big liquid slugs are formed in the right branch, which block the whole channel and result in ceased gas flow in the right branch. Some of the water in the right branch is able to flow through the first communication channel (the nearest one to the gas inlet), but it seems that water prefers to stay in the right branch, even flowing back to the entrance of the right branch. For the case with 1 mm size communication channels, the droplets or liquid slugs in the right branch are smaller than those in the 2 mm case, but still bigger than those in the base case. The communication channels allow the gas and liquid to exchange, but the left branch does not appear to gain much liquid via exchange. Liquid droplets in the right parallel channel after the third communication channel form liquid slugs due to the reduced gas flow. For the 0.5 mm case, more liquid begins to move through communication channels from the right to the left branch. All the communication channels are occupied by liquid due to the capillary effect in small channels, which also indicates that gas may not be able to flow through those channels. For the 0.25 mm case, more liquid  126  water is forced to flow from right to the left branch, and no liquid slugs show up in the last two cases.  Figure 6.2: Snapshots two-phase flow patterns in each case.  127  Water saturation is a key parameter, since liquid water formed in the gas flow channel will increase the pressure drop (causing a parasitic power loss) and may reduce PEM fuel cell performance. Figure 6.3 shows the time-averaged water volume fraction, namely water saturation, in both branches for the five cases. It is obvious that using communication channels will increase the water saturation in both branches. Since the liquid injection rates are the same for all those cases, it means that water has a longer residence time in communicating channel branches when there is communication between the two branches. Therefore, it will take longer time to have the liquid removed from the channels. The water volume fraction in the left branch also reflects the amount of liquid transferred through  the  communicating  channels.  Comparing  those  cases  with  communications, the 2 mm case is not recommended since it not only causes longer liquid slugs, but also significantly increases the water holdup in both branches. Meanwhile, smaller size communication channels are able to keep a  0.030 0.025 0.020 0.015 0.010 0.005 0.000  Base Case 2mm  1mm  0.5mm 0.25mm  Right Branch Volume Fraction (-)  Left Branch Volume Fraction (-)  lower liquid water holdup in the right branch and avoid the formation of slug flow. 0.6 0.5 0.4 0.3 0.2 0.1 0.0  Base Case 2mm  1mm  0.5mm 0.25mm  Cases  Cases Figure 6.3: Time-averaged water volume fraction in each branch.  The GDL water coverage ratio, which is defined as the ratio of the GDL surface area covered by water to the total GDL surface area for both parallel channels, is another key parameter to indicate the impact of liquid water on PEM fuel cell performance. In PEM fuel cells, gas reactants diffuse from the flow field channel through the GDLs to the catalyst layer. Once a certain area of GDL surface is  128  covered by liquid water, some of the diffusion pathways will be blocked, resulting in less catalyst sites being used and thus lowering the cell performance. Figure 6.4 shows the time-averaged total GDL water coverage ratio for the five cases. For the base case, the 0.5 mm case and the 0.25 mm case, the flow pattern is in top wall droplet flow with very little GDL surface area being covered by water. Therefore, they have almost the same level of GDL water coverage ratio. For the 1 mm case, due to the larger droplets and longer residence time, the two-phase flow pattern begins to transit from droplet flow to slug flow, which greatly increases the GDL water coverage. For even longer liquid slugs, as formed in the 2 mm case, the GDL water coverage ratio is much higher and keeps growing due to the continuous accumulation of water in the right branch. From these results, it can be concluded that more communication between two branches leads to higher water holdup and higher GDL surface water coverage, with both being harmful to the fuel cell performance. But when the size of communication channels is small enough, it may impose little influence on the fuel cell  Total GDL Water Coverage Ratio (-)  performance due to the unchanged GDL water coverage ratio.  0.30 0.25 0.20 0.15 0.10 0.05 0.00  Base Case 2mm  1mm  0.5mm 0.25mm  Cases Figure 6.4: Time-averaged total GDL water coverage ratio.  129  To quantify the effect of communication on the flow maldistribution in the two parallel channels, the time-averaged superficial gas velocity variation along the left branch is calculated and shown in Figure 6.5. Since the total gas flow rate for those cases are the same, it is not necessary to show the gas velocity in the right branch. For an even distribution, the gas flow velocity should be 2.232 m/s. Therefore, the closer the gas velocity is to 2.232 m/s, the less severe gas flow maldistribution there is in the parallel channel. For the base case, there is no communication between two parallel channels, thus the superficial gas velocity remains the same along the channel. The liquid injection rate difference is the only cause for the flow maldistribution. For other cases, communication channel size has a huge impact on the flow maldistribution. In the 1 mm case, gas flows from the right branch to the left branch at every communication channel, since the superficial gas velocity in the left branch always increases along the channel after each communication channel. Because the parallel channels require the pressure drop across each branch to be identical, the right branch, which has more water than the left branch, must have a lower gas velocity to balance the overall pressure drop across the parallel channel. On the other hand, gas is more mobile than liquid. Therefore, when a pressure difference is formed between the two branches, it is easier for gas to transfer from the higher pressure branch to the lower pressure branch. For the 2 mm case, wider communication channels allow more gas exchange, resulting in no gas passing through the right branch and the formation of long slugs in the right branch. Since liquid slugs cause large pressure drops, the pressure at the exit end of right branch may become less than the left branch, which pushes the gas flow back from the left to the right branch, i.e., the gas flow rate in the right branch increases after the 4th and 5th communication channels. It can be inferred that the high pressure drop caused by the liquid slug also results in less pressure gradient in the rest of the right branch to satisfy the requirement of equal total pressure drop in both branches. Small size communication channels, as in the 0.5 mm and 0.25 mm cases, have completely different hydrodynamics. Due to the capillary force, liquid in those two  130  cases fully fill up the communication channels, which can prevent the gas flowing through them. When there is a pressure difference, only liquid is allowed to exchange between the two branches. With 0.25 mm communication channels, the flow maldistribution is even more reduced comparing to the base case, which proves that the use of small communication channels can be an effective tool in  Superficial Gas Velocity along Left Branch (m/s)  mitigating the flow maldistribution.  5.0 4.5 4.0 3.5 3.0  Base case 0.25mm case 0.5mm case 1mm case 2mm case  2.5 2.0 0.00  0.02  0.04  0.06  0.08  0.10  X (m) Figure 6.5: Time-averaged superficial velocity along the left branch.  The time-averaged pressure drop for the five cases is shown in Figure 6.6. The pressure drops for all those cases almost remain the same and the 0.25 mm case gives the least pressure drop. As discussed in Chapter 5, the pressure drop seems to be insensitive to the severity of flow maldistribution. Therefore, using communication channels will not notably alter the overall pressure drop across the flow field channels.  131  Pressure Drop (Pa)  500 400 300 200 100 0  Base Case 2mm  1mm  0.5mm 0.25mm  Cases Figure 6.6: Time-averaged total pressure drop.  6.2.2. Effects of communication channels on liquid slug removal In an active PEM fuel cell, a liquid droplet or slug usually stays in the channel for a long time before it is removed from the channel by the passing gas, and the droplet or slug removal process is much faster than its accumulation. Therefore, the water removal in channels can be considered to be a transient process in the cell, which means that during the slug removal the water generated from reaction can be neglected. In order to investigate the effects of communicating channels on the water‟s transient behaviour, we started with a 1 cm long slug in the right branch as an initial state. During the whole process of liquid removal, no additional liquid water was injected into the cell. The final states of liquid distribution for those cases are simulated and shown in Figure 6.7. At the final states, liquid distribution will no longer change with time. For the base case, it is seen that all the liquid is flushed out of the channel after a certain time. But in both the 1 mm and 2 mm communication channel cases, some water still remains in the channel, indicating that the gas flow rate is not high enough to  132  expel water out of the channel. Furthermore, the remaining liquid in the 2 mm case forms liquid slugs rather than droplets hanging on the top wall, which has a negative impact on the fuel cell performance and the total pressure drop. Also, stagnant liquid droplets or slugs always hang around the communication channels, indicating that the gas exchange through the communication channels may hinder the liquid movement. In the 0.25 mm and 0.5 mm communication channel cases, almost all the communication channels are filled by liquid. There is still some water left in the parallel channels at the final state, but much less than that in the 1 mm and 2 mm cases, and the water does not cover the GDL surface but hangs on the top wall, which has little impact on the fuel cell performance.  133  Figure 6.7: Snapshots of initial and final water distribution in channels with communications.  Figure 6.8 shows the temporal variation of water volume fraction in both parallel channel branches. Liquid can be completely removed from the right branch after 0.45 s for the base case, but for the 1 mm case after about 0.5 s, there is still 20% of the initial water remaining in the right branch and another 4% in the left branch. Even worse, the 2 mm case appears to have no capability to remove the water, leaving about 90% of the initial water in the channel. Liquid removal in the 0.5 mm case is faster than the 1 mm and 2 mm cases. It takes about 0.45 s for the 0.4 mm case to reach the steady-state, but there is still 15% of the initial water left in the channel (excluding the water hanging in the communication channels). The small communication channel allows more liquid to flow from the right branch to the left branch. But since the liquid in the left branch mainly flows on a hydrophilic top wall, it tends to form a liquid film and takes much longer time to expel out of the channel. As a result, more water stays in the left channel branch at the steady state for the case with small 0.5 mm communication channels.  134  Right Branch Volume Fraction (-)  Left Branch Volume Fraction (-)  0.030 0.025  Base case 0.25mm case 0.5mm case 1mm case 2mm case  0.020 0.015 0.010 0.005 0.000 0.0  0.5  1.0  Time (s)  1.5  2.0  0.10  0.08  Base case 0.25mm case 0.5mm case 1mm case 2mm case  0.06  0.04  0.02  0.00 0.0  0.3  0.6  0.9  1.2  1.5  Time (s)  Figure 6.8: Temporal variation of water volume fraction in each parallel channel branch.  The temporal variation of total GDL water coverage ratio for the three cases is shown in Figure 6.9. It takes about 0.1 s to remove the water from the GDL surface for the base case, about 0.5 s for the 1 mm case and about 0.2 s for the 0.5 mm and 0.25 mm cases. Although there is certain liquid remaining in the channel for the 1 mm, 0.5 mm and 0.25 mm cases, the wettability difference between the GDL surface and the other walls is still large enough to lift the water up from the GDL surface, reducing the GDL water coverage ratio to almost zero. However, for the 2 mm case, too much water remains in the channel. After the liquid passing the second communication channel, most of the gas flows to the left branch, leaving liquid staying in the right channel near the 2nd communication channel to form liquid slugs. Therefore, the GDL water coverage ratio rebounds to about 1.5% after 0.2 s.  135  Total GDL Water Coverage Ratio (-)  0.05  0.04  Base case 0.25mm case 0.5mm case 1mm case 2mm case  0.03  0.02  0.01  0.00 0.0  0.1  0.2  0.3  0.4  0.5  0.6  Time (s) Figure 6.9: Temporal variation of total GDL water coverage ratio.  Figure 6.10 shows the final state of the superficial gas velocity profile along the left channel branch. As discussed previously, gas has a better mobility than liquid and can easily adjust itself to “find” the best way to reduce the total pressure drop. Therefore, when the communication channels are wide (1 mm and 2 mm cases), gas will always avoid passing by the liquid water if an alternative path is available. The velocity profile in Figure 6.10 clearly demonstrates this tendency. Referring to the final state of water distribution in Figure 6.7, more water in one branch always results in a higher gas flow rate in the other branch. When there is less or no water in one channel branch, the gas from the other channel branch will flow through the communication channel to it to balance the pressure drop in the other branch. This result is especially relevant for the 2 mm case, where the whole right branch from the second communication channel to the fifth communication channel is blocked by liquid slugs, causing the most severe gas maldistribution in this section. However, when the communication channels are narrow enough, as in the 0.5 mm and 0.25 mm cases, they are fully filled by  136  liquid which prevents the gas exchange via the strong capillary forces. Therefore,  Superficial Gas Velocity along Left Branch (m/s)  the gas flow rate along the channel remains the same.  5.0 4.5  Base case 0.25mm case 0.5mm case 1mm case 2mm case  4.0 3.5 3.0 2.5 2.0 0.00  0.02  0.04  0.06  0.08  0.10  X Figure 6.10: Final state of superficial gas velocity along the channel length.  From those simulation results, it can be concluded that when a liquid slug passes through parallel channels with communication channels, there will be certain liquid exchange between the two parallel channel branches. The amount of liquid communicated depends on the size of communication channels, i.e., the smaller the channel, the more liquid exchange. Communicating parallel channels always delay the liquid removal from the channel and subsequently liquid expelling from the GDL surface, which would be potentially detrimental to the fuel cell performance.  6.3. Implications of communication channels to PEM fuel cells Communication is always happening inside a PEM fuel cell GDL via in-plane diffusion or convection [156]. While the communication channel geometry may be  137  exaggerated in our simulations, flow fields relying on in-plane GDL diffusion (such as interdigitated flow fields) may need to consider communication-induced maldistribution. This mechanism of maldistribution (in addition to GDL intrusion or water formation in gas flow channels for instance) is new to the literature, and may be a consideration for future experimental and numerical work. High in-plane diffusion may also cause reactant bypass (reducing active area used) in PEM fuel cells with certain flow fields [157]. Serpentine flow field channels experience flow bypass when there is a difference between the local oxygen partial pressure and the total pressure of an adjacent channel. Under these conditions, convective and diffusive transport of gas reactants under the landing area can occur (adjacent channels thus establish communications). While the reactant bypass in the GDL was reported to enhance the electrochemical reaction for the case with a high in-plane permeability [157], the communication-induced maldistribution was ignored in the previous simulation. The negative effects of the maldistribution may offset any electrochemical gains by reducing flooding in the flow field channels. Thus, the results presented here could be incorporated into future simulations for a more accurate understanding of the relationship between water management and PEM fuel cell performance. A slotted-interdigitated flow field plate was recently proposed to eliminate flow maldistribution [158]. These slots connected typical interdigitated channels in a way similar to the communicating channels studied in this paper. An analytical model based on single-phase gas flow was used to determine the optimal slot dimensions, without accounting for the possible presence of liquid water in the flow field channels. Thus, the conclusions obtained from this study may not be applicable to the actual fuel cells when two-phase flow is always present. In the presence of liquid water, we expect that the maldistribution may actually be worse in their proposed geometry since their communication channels are much wider than those of the flow channels.  138  As a potential method to mitigate the flow maldistribution in multiple parallel channels, communicating channels need to be carefully optimized, since they can also cause certain negative effects such as the reduced liquid removal capacity. More experimental or numerical work should be conducted to further determine the optimal size, number and positions of the communication channels in parallel flow channels.  6.4. Conclusions In this study, the VOF method was employed to investigate the performance of communicating  flow  field  channels  on  mitigating  the  two-phase  flow  maldistribution. The following conclusions can be drawn, 1) Communication between parallel channels can have a great impact on the two-phase flow pattern, gas and water distribution and flow maldistribution. When the communication channels are wide, they provide a pathway for gas to shortcircuit the liquid, leading to a worse gas flow distribution. However, when the communication channels are narrow enough, due to the capillary force, liquid will occupy the whole communication channel which prevents the gas exchange between parallel channels. As a result, the small communication channels are helpful for mitigating the flow maldistribution by redistributing the liquid among the parallel flow channels. 2) With communication channels, the liquid water in the flow channel tends to have a longer residence time, which results in a higher water content in both branches, and higher water coverage ratio on the GDL surface, which, however, has little impact on the overall pressure drop. 3) The removal of liquid slugs in the parallel channels is significantly affected by the communication channels as well. With communications, water removal from flow channels or from the GDL surface is delayed, which would decrease the  139  PEM fuel cell performance. Wide communication channels even make water permanently remain in the channel, leading to increased total pressure drop. 4) The current simulation results suggest that wide communication channels between parallel flow channels should be avoided to reduce gas maldistibution and accelerate liquid water removal. Narrow communication channels can be used to mitigate the flow maldistribution but need to be carefully designed. 5) Since the communication between parallel channels naturally occurs in a PEM fuel cell through the porous GDL via in-plane diffusion, further simulations should be carried out to investigate the effects of this inherent communication.  140  7. 3D simulations of the impact of two-phase flow on PEM fuel cell performance1 To better understand the impact of two-phase flow on the PEM fuel cell performance, a two-phase flow model must first be able to reproduce the twophase flow patterns at the channel length scale. As discussed in Chapter 1, the VOF method is perhaps the most potential model because of its capacity to track the gas-liquid interface at large scales. However, coupling the VOF method with electrochemical reaction is still a challenge. Le et al. [75, 76, 106] developed a general model that included electrochemical reactions, heat transfer and species transport in the frame of VOF method. But, they treated the MEA as a homogenous layer without a microstructure, with the liquid droplets being initially located in gas flow channels. Therefore, such an approach can only provide some initial or transient information about the impact of droplets on PEM fuel cell performance. To address this issue, we have developed a 1+3D two-phase flow model as presented in Chapter 2. The VOF method was employed to account for the two-phase flow and species transports in a 3D cathode-side gas flow channel, coupling with a 1D MEA model to account for the electrochemical reactions and species transport inside MEA. In this chapter, the two-phase flow patterns under severe flooding conditions are examined based on simulations. The effects of liquid generation rates, current densities, and gas stoichiometric flow ratios on the PEM fuel cell performance are investigated. In addition, the impact of GDL surface and channel wall wettability on the cell performance is also studied.  1  A version of this chapter has been submitted for publication. Y. Ding, H.T. Bi, D.P. Wilkinson. 3D Simulations of the Impact of Two-phase Flow on PEM Fuel Cell Performance.  141  7.1. Computational domain and parameters A three-dimensional straight single channel was used in all simulations. As shown in Figure 2.1, the channel has a 1.0 mm × 1.0 mm square cross section and 100 mm in length with a hydrophobic GDL surface on the bottom and three hydrophilic channel walls. 20 water emerging pores are located randomly on the GDL surface along the channel with the same diameter of 250 µm, which represents the large pores on the GDL surface. A total of 100,000 meshes with a uniform size of 0.1 mm are used in the simulation, which is fine enough to obtain the stable two-phase flow patterns at the channel scale. Adaptive time step method is utilized in order to keep the global courant number less than 0.5, which ensures that the droplet behaviour in the channel can be well captured.  Figure 7.1: Illustration of 1+3D cathode side of a PEM fuel cell.  Fully humidified air flows into the channel from one end and both liquid and gas exit at the other end. The air flow rate at the inlet is calculated based on the gas stoichiometric ratio  , which is defined as follows:  142    mo2 ,inlet mo2 ,reaction    f o2 ,inlet  mair ,inlet I op  S MEA 4F  (7.23)  where mo2 ,inlet is the number of moles of oxygen at the inlet and mo2 ,reaction is the theoretically needed oxygen calculated by the operating current density and active surface area of MEA. The physical properties of oxygen, nitrogen and water (vapor and liquid), such as density, viscosity and diffusivity, are all calculated based on the kinetic theory at a given operating temperature and pressure [159]. Other operating conditions and model parameters are given in Table 7.1.  Table 7.1: Constants and operating parameters. Constants and operating parameters  Value  Gas constant  R  8.314 J  K 1  mol1  Faraday constant  F  96500 C  mol1  Henry‟s constant for oxygen [127]  H o  20265.0 Pa  m3  mol1  Open circuit potential [127]  E0  0.95 V  Thickness of GDL [127]  LG  2.5 104 m  Cathode transfer coefficient [127]   1 c  1  Cathode oxygen reference concentration [127]  Co ref  40.9 mol  m-3  Oxygen diffusion coefficient [127]  Do2 eff  7.2043 107 m2  s-3  Cathode exchange current density [127]  io,c  426.0 A  m2  Membrane water transfer coefficient [127]  Static contact angle [154]    9.37 106 S  m2  w,side  45o ,w,GDL  140o  Operating pressure  Pop  206703.0 Pa  Operating temperature  T  348.15 K  Operating current density  I op  500 ~ 12000 A  m2  Operating gas stoichiometric ratio    2, 3, 4, 5, 6  Inlet air humidity  100%  Effective conductivity [127]  143  7.2. Results and discussion 7.2.1. Effects of stoichiometric flow ratios for single phase flow A series of base cases were first simulated at gas stoichiometric flow ratios ranging from 2 to 6. Water was assumed to be finely dispersed in the gas phase only, which corresponds to the dispersed two-phase flow, which is analogy to single-phase flow. It is seen in Figure 7.2 that the gas stoichiometric flow ratio or the gas flow rate has great impacts on the performance only in the mass transport limited region at current densities more than 800 mA/cm2 for those simulated cases.  1.0  Voltage (V)  0.8  0.6  stoich.=2 stoich.=3 stoich.=4 stoich.=5 stoich.=6  0.4  0.2  0.0 0  200  400  600  800  1000  1200  2  Current Density (mA/cm ) Figure 7.2: Effects of gas stoichiometric flow ratios for single phase flow.  7.2.2. Effects of two-phase flow on the PEMFC performance Liquid water in the gas flow channel blocks the pathway of reactants diffusing through the GDL, and alters the gas distribution along the channel. The twophase flow patterns are not only determined by how much of water the channel  144  contains, namely the water volume fraction, but also depends on how water distributes in the channel as shown in Chapter 3. Figure 7.3 shows the twophase flow pattern in the cathode channel when the cell is operated at 800 mA/cm2 with a gas stoichiometric flow ratio of 2. Several small droplets can be identified on the GDL surface. However the water volume fraction is only about 0.7% on average, which is so little to give only a decrease of 4.3 mV in the cell voltage compared to the base case with single phase flow in the cathode channel.  Figure 7.3: Two-phase flow pattern in the cathode side channel.  From experiments [15, 94, 98, 99], the water content in the cathode channel is generally much higher than the above value. One possible reason is that the channel length we are simulating is only 10 cm, which is much shorter than a real fuel cell channel which typically has a length scale of 1 m. Therefore, less water is generated and accumulated in the end section of the channel. Another possible reason is that in real fuel cell channels, due to the channel wall roughness, some water may remain in the channel and cannot be readily removed by the passing gas. On contrary, in the current simulation, the channel walls are assumed to be perfectly smooth. Also, it usually takes several minutes  145  or even hours to observe slug flow in the experiment, but it takes only 10 s to reach the steady state in our simulation. It is worthy to note that simulating a longer channel with rough surfaces or simulating over a longer time requires extremely huge computational time, which is unrealistic for the current computational capacity. One alternative way to deal with this issue is to simply increase the water generation rate, which has been used in several simulation studies [100, 154, 160]. It has been proved that at a small scale, the two-phase flow pattern does not change even when the flow rate is amplified by 1000 times [154]. Therefore, we amplified the water generation rated by 10 (X10 case), 100 (X100 case), and 1000 times (X1000 case), respectively, in subsequent simulations presented below, while the oxygen consumption rate is kept the same. The two-phase flow patterns for the three cases are shown in Figure 7.4. For the X10 case, in addition to the appearance of several liquid droplets on the GDL surface, a larger droplet forms and flows on the top wall. This flow pattern, namely top wall droplet flow, is in accordance with the two-phase flow pattern reported in Chapters 3 and 4. When small droplets coalesce to form large droplets, the opportunity for the droplet to touch the hydrophilic side walls is increased. Therefore, once hit the side wall, the droplets will detach them from the hydrophobic GDL surface and then flow on to the top wall. The same mechanism also applies to the X100 case, except that there are more droplets on the top wall, which begin to coalesce and spread on the top wall forming liquid films. When the liquid water generation rate is further increased, as in the X1000 case, the droplets on the top wall become so big that they span to the whole cross sectional area of the channel and finally transform into liquid slugs.  146  (a)  (b)  (c) Figure 7.4: Effects of liquid water generation rates on two-phase flow patterns. (a) 10 times, (b) 100 times, (c) 1000 times.  147  Comparing the fuel cell performance at different liquid generation rate shown in Figure 7.5, it is found that only in the X1000 case, there is a drastic cell voltage drop from the X100 case. The error bar in the figure shows the standard deviation, which indicates the fluctuations over time. The pressure drop also does not change too much in the X1, X10 and X100 cases, since the droplets on the top and side hydrophilic walls tend to form liquid films, which occupy less cross sectional areas than the droplets on the GDL surface. However, the formation of liquid slugs, as in the X1000 case, causes a significant increase of pressure drop. Meanwhile the large fluctuations of pressure drop and cell voltage must have a damaging effect on the cell stability and durability. Therefore, it is suggested that the slug flow should be avoided during PEM fuel cell operation.  400  0.7  300 0.5  Voltage Pressure Drop  200  0.4 0.3  Voltage (V)  Pressure Drop (Pa)  0.6  100 0.2 0  0.1 Single Phase  X1  X10  X100  X1000  Cases Figure 7.5: Effects of liquid water generation rates on the cell voltage and pressure drop.  The liquid water distribution on different walls and liquid volume fraction in the channel are given in Figure 7.6. A water coverage ratio on a surface is defined as the ratio of the surface area covered by liquid water to the total surface area. Clearly, the fuel cell performance is greatly influenced by the water coverage  148  ratio of the GDL surface, because it determines the effective mass transfer area and reaction area in the catalyst layer. For the X1, X10 and X100 cases, the amount of liquid water on the GDL surface changes very little, as shown in Figure 7.6, since most of liquid droplets flow on the top wall. Meanwhile, the water coverage ratio on the top and side walls increases rapidly with increasing the liquid generation rate, indicating that the extra liquid water mainly flows on the top and side walls. Although the water volume fraction increases accordingly, the overall cell voltage and pressure drop do not change much. However, further increase in liquid generation rate makes the flow pattern transition to the slug flow, and the GDL water coverage ratio increases significantly, resulting in a lower cell voltage. All the above results imply that due to the wettability differences between the GDL surface and other walls, the gas flow channel can hold a wide range of water amount without triggering the deterioration of the fuel  Water Coverage Ratio and Volume Fraction (-)  cell performance.  0.4 GDL Water Coverage Ratio Top-wall Water Coverage Ratio Side-wall Water Coverage Ratio Volume Fraction  0.3  0.2  0.1  0.0 X1  X10  X100  X1000  Cases Figure 7.6: Effects of liquid water generation rates on liquid wall coverage and volume fraction.  149  7.2.3. Effects of severe flooding on cell performance Slug flow is inevitable under high loading operating conditions, as observed from the experiment. Thus, it is important to investigate how to mitigate its negative impact on the PEM fuel cell performance. In this section, the fuel cell was operated at different current densities, ranging from 50 mA/cm2 to 1000 mA/cm2, i.e. 50, 100, 200, 400, 600, 800, 1000 mA/cm2. The water generation rate was amplified 1000 times as in the X1000 case in order to study the impact of slug flow, which corresponds to a PEM fuel cell operating under severe flooding which usually occurs at the end section of the flow channel. The fuel cell polarization curve in Figure 7.7 shows that the slug flow only affects the fuel cell performance in the mass transport limited region when the current density is more than 600 mA/cm2. 1.0  Voltage (V)  0.8  0.6  Single Phase Two phase  0.4  0.2  0.0 0  200  400  600  800  1000  2  Current Density (mA/cm ) Figure 7.7: Fuel cell polarization curves under severe flooding conditions.  At 800 mA/cm2, the presence of liquid slugs causes great voltage fluctuations, which are detrimental to the fuel cell power output and its durability. Further increasing the operating current density results in the complete shutdown of the cell, since there is not enough oxygen at the catalyst layer to sustain the reaction.  150  When operating in kinetic and ohmic region, the liquid slug flow has little effects on the fuel cell performance, and the fuel cell output voltage can still remain stable. Figure 7.8 shows the overall pressure drop as a function of current densities. The slug flow causes large pressure drop and fluctuations compared to single-phase flow, even at low current densities. Therefore, to reduce the parasitic energy loss and stabilize fuel cell operation, slug flow should be avoided, especially at high current densities.  Pressure Drop (Pa)  400  Single Phase Two phase  300  200  100  0 0  200  400  600  800  1000  2  Current Density (mA/cm ) Figure 7.8: Effects of operating current densities on overall pressure drop.  Figure 7.9 shows the water distribution on different surfaces and water volume fraction in cathode channels. Since the gas stoichiometric flow ratio is fixed in these cases, namely the ratio between gas flow rate and water generation rate is kept the same at different current densities, changing operating current density does not have great impacts on water wall coverage and distribution among the walls and water volume fraction. At low current densities, increasing the current  151  density results in slightly higher water coverage ratios on top and side walls but lower water coverage ratio on the GDL surface, implying that the increased drag force between two phases, which acts to lift the liquid from the GDL, overweighs the increased water generation rate. One exception is the last point at 1000 mA/cm2, which is due to a reduced water generation rate due to the fuel cell shutdown. Therefore, there should exist an optimal current density, where lowest GDL water coverage ratio can be achieved.  0.5 GDL Water Coverage Ratio Top-wall Water Coverage Ratio Side-wall Water Coverage Ratio Volume Fraction  Ratio (-)  0.4  0.3  0.2  0.1  0.0 0  100  200  300  400  500  600  700  800  900  2  Current Density (mA/cm ) Figure 7.9: Effects of operating current densities on water wall coverage and liquid holdup.  7.2.4. Effects of gas stoichiometric flow ratio on the PEMFC performance To mitigate the negative impact of slug flow, one can increase the gas stoichiometric flow ratio, i.e. the gas flow rate to accelerate the slug flush-out. In this section, several gas stoichiometric flow ratios, i.e. 2, 3, 4, 5 and 6, were used in the simulation at a constant current density of 800 mA/cm2. The two-phase flow patterns for different stoichiometric flow ratios are shown in Figure 7.10. As discussed above, slug flow is identified at a stoichiometric flow ratio of 2.  152  However, when operating at higher stoichiometric flow ratios, the increased drag force between two phases reduces the droplet residence time in the channel, and thus reduces the liquid slug length. Further increasing gas flow rate results in the two-phase flow pattern to shift from slug flow to the top wall droplet flow. The latter flow pattern is helpful to maintain high performance since there is little liquid water left on the GDL surface, and it also helps to reduce the overall pressure drop since the channel is no longer blocked by liquid slugs.  (a)  (b)  (c)  (d)  153  (e) Figure 7.10: Two-phase flow patterns at different gas stoichiometric flow ratios. (a) stoich.=2 (b) stoich.=3 (c) stoich.=4 (d) stoich.=5 (e) stoich.=6.  The fuel cell output voltage and pressure drop over the channel at different gas stoichiometric flow ratios are shown in Figure 7.11. The pressure drop is seen to increase linearly with increasing the gas stoichiometric flow ratio. However, when the gas stoichiometric flow ratio becomes equal to or larger than 3, the fuel cell voltage does not change much, and the channel remains in slug flow. The reason is that more gas flow into the channel not only affects the two-phase flow pattern, but also provides more oxygen so that the fuel cell is not operated in the mass transport controlled region any more. Therefore, in this specific case, to maximize the fuel cell efficiency, the optimal gas stoichiometric flow ratio is around 3, since using larger stoichiometric ratios only causes more parasitic energy loss with little gain in voltage output.  154  800  0.8  600  0.6  400  0.4  200  0.2  0  Voltage (V)  Pressure Drop (Pa)  Voltage Pressure Drop  0.0 1  2  3  4  5  6  7  Stoich. (-) Figure 7.11: Effects of gas stoichiometric flow ratio on the fuel cell voltage output and pressure drop.  The water coverage and distribution on channel walls and water volume fraction at different gas stoichiometric flow ratios are given in Figure 7.12. As the stoichiometric ratio increases, the increase in gas flow rate reduces the droplet residence time, leading to lower water coverage ratios on the walls and the water volume fraction. One exception is at the stoichiometric ratio of 6 where a higher top wall water coverage ratio indicates that the flow pattern begins to transform from slug flow to the top wall droplet flow.  155  0.5 GDL Water Coverage Ratio Top-wall Water Coverage Ratio Side-wall Water Coverage Ratio Volume Fraction  Ratio or Fraction (-)  0.4  0.3  0.2  0.1  0.0 1  2  3  4  5  6  7  Stoich. (-) Figure 7.12: Effects of gas stoichiometric flow ratio on water wall coverage ratio and volume fraction.  According to the previous single-phase flow results, the voltage differences in both kinetic and ohmic regions at different gas flow rate should be very small. Fuel cell operation at different gas stoichiometric flow ratios and different current densities were simulated only in the mass transport limited region, ranging from 800 mA/cm2 to 1200 mA/cm2. As shown in Figure 7.13, increasing the gas flow rate significantly extends the ohmic region, as in the single-phase flow cases, and the effect becomes less significant when the gas flow rate is further increased at high stoichiometric ratios.  156  0.8  Voltage (V)  0.6  0.4  0.2  0.0 400  stoich.=2 stoich.=3 stoich.=4 stoich.=5 stoich.=6  600  800  1000  1200  2  Current Density (mA/cm ) Figure 7.13: Effects of gas stoichiometric flow ratio on polarization curves.  7.2.5. Effects of GDL surface wettability on fuel cell performance To investigate the impact of GDL surface wettability on the two-phase flow patterns and fuel cell performance under severe flooding, several GDL surface contact angles, i.e., 45o, 60o, 120o, and 140o were tested in the simulation. The contact angle of the other walls was kept constant at 45o. While the cell is operated at 800 mA/cm2, it is found that the two-phase flow pattern remains in slug flow with different GDL surface contact angles as shown in Figure 7.14. However, with a hydrophobic GDL surface, (Figure 7.14 c, d) the liquid slug formation is later than with a hydrophilic GDL surface (Figure 7.14 a, b). There is less reaction occurring at the end of the channel due to the reactants depletion, indicating that increasing GDL surface contact angle is an effective method to yield better cell performance.  157  (a)  (b)  (c) (d) Figure 7.14: Effects of GDL surface contact angle on the two-phase flow patterns in the channel at 800 mA/cm2: (a) θ=45o, (b) θ=60o, (c) θ=120o, (d) θ=140o.  Changing the GDL surface wettability can effectively alter the water distribution in the fuel cell as shown in Figure 7.15. The GDL water coverage ratio significantly decreases with increasing GDL hydrophobicity. More water is lifted up due to the wettability difference between the GDL and other walls and flows on the top wall, leading to a higher top wall water coverage ratio. For more hydrophilic GDL surface, there are more liquid slugs formed, which accelerates the liquid removal from the cell and leads to the slightly lower water volume fraction.  158  0.5 GDL Water Coverage Ratio Top-wall Water Coverage Ratio Side-wall Water Coverage Ratio Volume Fraction  Ratio or Fraction (-)  0.4  0.3  0.2  0.1  0.0 40  60  80  100  120  140 o  GDL Surface Contact Angle ( ) Figure 7.15: Effects of GDL surface contact angle on water wall coverage and holdup at 800 mA/cm2.  As discussed in previous chapters, the pressure drop depends on not only the liquid holdup, but also the flow pattern and water distribution. When the GDL is hydrophilic, increasing its contact angle accelerates the liquid slug removal, resulting in lower pressure drops. But when the GDL is hydrophobic, increasing its contact angle pushes droplet to the hydrophilic channel wall causing longer water residence time, and thus increases the pressure drop.  159  Pressure Drop (Pa)  400  300  200  100  0 40  60  80  100  120  140 o  GDL Surface Contact Angle( ) Figure 7.16: Effects of GDL surface contact angle on the fuel cell pressure drop at 800 mA/cm2.  The fuel cell output voltage under three current densities, i.e., 600 mA/cm2, 800 mA/cm2 and 1000 mA/cm2 with different GDL surface contact angles is shown in Figure 7.17. Those three current densities represent the cell operations at ohmic and mass transfer limited region. Clearly, using a more hydrophobic GDL surface is highly recommended because it can significantly extend the ohmic region and increase the fuel cell performance. As we have found in previous chapters, hydrophobic GDL is helpful to lift the droplet up, leaving more area for reactant diffusing through it, therefore is capable to improve the fuel cell performance.  160  0.8 0.7  Voltage (V)  0.6 0.5  @600 mA/cm2 @800 mA/cm2 @1000 mA/cm2  0.4 0.3 0.2 0.1 0.0 40  60  80  100  120  140 o  GDL Surface Contact Angle ( ) Figure 7.17: Effects of GDL surface contact angle on the fuel cell voltage output at three different current densities.  7.2.6. Effects of channel wall wettability on fuel cell performance The impact of channel wall wettability on the two-phase flow patterns and fuel cell performance under severe flooding conditions were simulated by changing the channel wall contact angles from 45o to 140o, i.e., 45o, 60o, 120o, and 140o. The contact angle of GDL surface was kept constant at 140o. It is found that using more hydrophobic channel walls can reduce the liquid slug size but increase the number of slugs (Figure 7.18 c, d). Top wall droplet flow almost disappears in the channel with hydrophobic channel walls.  161  (a)  (b)  (c) (d) Figure 7.18: Effects of channel wall contact angle on the two-phase flow patterns in the channel at 800 mA/cm2: (a) θ=45o, (b) θ=60o, (c) θ=120o, (d) θ=140o.  Using more hydrophobic channel walls also reduces liquid holdup and the top and side wall water coverage ratios (Figure 7.19), indicating that it is easier for the liquid to be removed from the channel. For the hydrophobic channel wall, the GDL water coverage ratio is higher since there is not enough wettability difference to life the liquid up. But for the hydrophilic channel wall, increasing GDL contact angle can slightly reduce the GDL water coverage ratio, which may be attributed to the faster liquid removal from the channel.  162  0.5 GDL Water Coverage Ratio Top-wall Water Coverage Ratio Side-wall Water Coverage Ratio Volume Fraction  Ratio or Fraction (-)  0.4  0.3  0.2  0.1  0.0 40  60  80  100  120  140 o  Channel Wall Surface Contact Angle ( ) Figure 7.19: Effects of channel wall contact angle on water wall coverage and holdup at 800 mA/cm2.  The pressure drop over the channel with hydrophilic channel walls is higher than that with hydrophobic channel walls due to the longer liquid resident time, but is insensitive to the channel wall contact angle as long as the channel walls remain hydrophilic (Figure 7.20). Since there is a flow pattern transition from the top wall droplet flow to the slug flow when a channel wall changes from hydrophobic to hydrophilic, the pressure drop over the channel with hydrophilic walls has higher fluctuations than that with hydrophobic walls.  163  Pressure Drop (Pa)  400  300  200  100  0 40  60  80  100  120  140 o  GDL Surface Contact Angle( )  Figure 7.20: Effects of channel wall contact angle on the fuel cell pressure drop at 800 mA/cm2.  As shown in Figure 7.21, the impact of channel wall wettability on the fuel cell voltage output is not as remarkable as the GDL surface wettability. In the ohmic region, i.e., 600 mA/cm2, the effect is insignificant. In the mass transport limited region, i.e., 800 mA/cm2, a hydrophobic channel wall increases the GDL surface water coverage thus reduces the cell output voltage, and is thus not recommended.  164  0.8 0.7  Voltage (V)  0.6 0.5 0.4 0.3  @600 mA/cm2 @800 mA/cm2 @1000 mA/cm2  0.2 0.1 0.0 40  60  80  100  120  140 o  Channel Wall Surface Contact Angle ( ) Figure 7.21: Effects of channel wall contact angle on the fuel cell voltage output at three different current densities.  7.3. Conclusions The following conclusions can be drawn from this simulation work: 1) Three flow patterns are observed in an operating fuel cell, namely, droplet flow, top wall droplet/film flow and slug flow. To obtain slug flow in the simulation, liquid generation rate is amplified by 1000 times, which represents the severe flooding observed at the end section of a long flow channel or after an extended operating period of time. 2) Under severe flooding, compared to single-phase flow, the slug flow pattern only affects fuel cell voltage output in the mass transport limited region, i.e. at high current densities, and causes large cell voltage decrease and fluctuations, which are detrimental to the fuel cell stability and durability. However, the slug flow always creates high pressure drops and large pressure fluctuations, causing  165  extra energy loss. Therefore, the slug flow should be always avoided in operating PEM fuel cells. 3) Increasing gas stoichiometric flow ratios can effectively extends the stable fuel cell operating range, and may alter the two-phase flow pattern in the channel. However, when operating at a certain current density, too large a stoichiometric ratio is not recommended, because it may cause more parasitic energy loss with little gain in voltage output. 4) Altering GDL surface wettability can significantly change the water coverage ratio and distribution on channel walls and the fuel cell performance. Using more hydrophobic GDL surfaces is recommended to extend the ohmic region and increase the fuel cell output voltage. Although the impact of channel wall wettability on the cell performance is not as remarkable as the GDL surface wettability, increasing the contact angle of a hydrophilic channel wall can slightly increase the fuel cell output voltage because of the reduced GDL water coverage ratio.  166  8. Effects of two-phase flow maldistribution on PEM fuel cell performance In this chapter, the 1+3D two-phase flow model is applied to a fuel cell consisting of two parallel flow channels to investigate the effects of flow maldistribution on the fuel cell performance. The total water flow rate is determined by the current density, and the flow maldistribution is introduced by varying the water distribution in the two branches. The effects of flow maldistribution severity and liquid magnification are discussed in detail. Furthermore, several maldistribution mitigation methods are tested, such as changing material wettability, increasing gas stoichiometric flow ratio and adding inlet resistances.  8.1. Computational domain and parameters The three-dimensional computational domain of the two parallel channels is shown in Figure 2.1. Each channel branch has a 1.0 mm × 1.0 mm square cross section and is 50 mm in length. The distance between two branches is 1 mm. Similarly to the base case in last chapter, there are a hydrophobic GDL surface with a contact angle of 140o on the bottom and three hydrophilic channel walls with a contact angle of 45o. 10 water emerging pores in each branch are located randomly on the GDL surface along the channel with the same diameter of 250 µm, which represents the large pores on the GDL surface. The position of pores in each branch is symmetrical in order to avoid geometry-introduced flow maldistribution as discussed in Chapter 5. Fully humidified air flows into the channel from one end and both liquid and gas exit at the other end. The gas stoichiometric ratio is set at 2 for the base case, and all other parameters, such as reactant physical properties, operating conditions and model parameters are the same as in the last chapter. The total water flow rate is calculated by the operating current density, but for the base case, it is amplified by 1000 times in order to represent the severe flooding at the end section of a long channel or after an extended operating period of time. To introduce flow maldistribution,  167  different amount of water is injected into branches A and B. For example, the „25/75 case‟ indicates 25 percent of total water generated by reaction being injected into branch A and the other 75% into branch B. The „0/100 case‟ gives the extreme flow maldistribution with no water entering branch A, and the „50/50 case‟ represents the even distribution of liquid water. A total of 116,000 meshes with a size of 0.1 mm are used in the simulation, which is fine enough to obtain the stable two-phase flow patterns at the channel scale. Adaptive time step method is utilized in order to keep the global courant number less than 0.5, which ensures that the droplet behavior in the channel can be well captured.  Figure 8.1: Three-dimensional computation domain of two parallel channels.  8.2. Results and discussion 8.2.1. Effects of flow maldistribution severity on fuel cell performance Figure 8.2 shows the two-phase flow patterns under different flow maldistribution severity and at different current densities. As discussed in previous chapters, the  168  two-phase flow pattern in typical PEM fuel cell channels evolves from droplet flow to corner flow, top wall film flow, and finally slug flow. Since the liquid generated rate is amplified by 1000 times, the main two-phase flow pattern in the channel is the slug flow, but other two-phase flow patterns still can be identified in Figure 8.2.  (a)  (b)  (c)  Figure 8.2: Effects of flow maldistribution severity on the two-phase flow patterns in the channel under three different current densities: (a) 100 mA/cm2, (b) 400 mA/cm2, (c) 800 mA/cm2.  Obviously, more severe flow maldistribution leads to longer liquid slug, and in some cases, i.e., the „0/100 case‟ at 400 and 800 mA/cm2, and the „25/75 case‟ at 800 mA/cm2, the slug even occupies the entire branch, which will completely  169  shut down the reaction in it. For low current density, since the gas flow rate is low, severe flow maldistribution results in back flow in the branch B as in the „0/100 case‟ at 100 mA/cm2. The fuel cell polarization curves under different degree of flow maldistribution are shown in Figure 8.3. For moderate flow maldistribution, namely the 25/75 case, only mass transport limited region is affected by the flow maldistribution. Severe flow maldistribution, namely the 0/100 case, results in a huge voltage decrease in both ohmic and mass transport limited regions. Therefore, from the fuel cell performance point of view, flow maldistribution should be well avoided.  0.9 0.8 0.7  Voltage (V)  0.6 0.5 0.4  50/50 25/75 0/100  0.3 0.2 0.1 0.0 0  200  400  600  800  1000  2  Current Density(mA/cm ) Figure 8.3: Effects of flow maldistribution severity on the fuel cell polarization curve.  The pressure drop differences caused by flow maldistribution at the same current density are not obvious at low to moderate maldistribution especially at low current densities, as shown in Figure 8.4. At severe maldistribution with the cell operated at high current densities, the pressure drop becomes higher than the base case with uniform flow maldistribution.  170  Pressure Drop (Pa)  180  135  90  50/50 25/75 0/100  45  0 0  200  400  600  800  1000  2  Current Density(mA/cm ) Figure 8.4: Effects of flow maldistribution severity on the total pressure drop.  Figure 8.5 shows that the flow maldistribution also increases the total GDL water coverage ratio in the whole range of current densities, which is the main cause for poor fuel cell performances. In 4 simulated cases almost 50% GDL surface is covered by the water, indicating the likelihood of branch B being fully flooded.  171  GDL Water Coverage Ratio (-)  0.6 0.5 0.4  50/50 25/75 0/100  0.3 0.2 0.1 0.0 0  200  400  600  800  1000  2  Current Density(mA/cm ) Figure 8.5: Effects of flow maldistribution severity on the GDL water coverage ratio.  The gas velocity ratio in Figure 8.6 is defined by the gas velocity in branch B over that in branch A. A value of 1.0 represents the even gas flow distribution, and there is no gas flow through branch B at a gas velocity ratio of 0.0. Apparently, the 0/100 case leads to the most severe gas flow maldistribution. For the 50/50 case, there is a small deviation from 1.0 when the operating current density is higher than 100 mA/cm2. However, at 50 mA/cm2, the deviation is substantial, which may be due to the flow instabilities when the pressure drop is sensitive to the two-phase flow patterns at low gas velocity. For example, once a liquid slug is firstly formed in one branch, it increases the flow resistance there, and the gas flow rate in the other branch will increase, which thus reduces the chance of liquid slug formation in it. Especially at low current densities, there is not enough liquid in the channel to form slug flow in both branches. Therefore, even though the liquid flow rates are same in each branch, flow maldistribution still occurs. The flow maldistribution at low current densities for the „0/100 case‟ and „25/75 case‟ is less severe than at high current densities owing to the back flow shown in Figure 8.2.  172  Gas Velocity Ratio (B to A) (-)  1.6 1.4 1.2 1.0 0.8  50/50 25/75 0/100  0.6 0.4 0.2 0.0 0  200  400  600  800  1000  2  Current Density(mA/cm ) Figure 8.6: Effects of flow maldistribution severity on the gas velocity ratio  8.2.2. Effects of liquid generation rate on flow maldistribution In the previous section, the liquid generation rate was amplified by 1000 times, which only represents the severe flooding at the end section of a long channel or after an extended operating period of time. To investigate the effects of flow maldistribution on fuel cell performance under low liquid generation rate, we simulated the cell performance with the liquid generation rate amplified by 10 and 100 times, namely the X10 case and the X100 case. The output voltage at a current density of 800 mA/cm2 for the three cases with different liquid generation rates and different degrees of flow maldistribution is shown in Figure 8.7. Without flow maldistribution, namely the 50/50 case, there is a small voltage drop from the X100 case to the X1000 case, indicating that the fuel cell enters the mass transport limited region at high liquid generation rate. Meanwhile, with the flow maldistribution, namely both the 25/75 and the 0/100 cases, the fuel cell almost shuts down, even at low liquid generation rates, which  173  proves that the flow maldistribution in parallel channels has a detrimental impact on the fuel cell performance and should be well avoided.  0.6 0.5  Voltage (V)  0.4 0.3  50/50 25/75 0/100  0.2 0.1 0.0 X10  X100  X1000  Cases Figure 8.7: Effects of liquid generation rate on the fuel cell output voltage at a current density of 800 mA/cm2.  Gas velocity ratios at different liquid generation rates are shown in Figure 8.8. It is found that the gas flow maldistribution is slightly less severe at lower liquid generation rates, indicating that at higher current densities, the gas distribution in parallel channels is more sensitive to the liquid distribution. Even slight liquid uneven distribution may cause great gas flow maldistribution.  174  Gas Velocity Ratio (B to A) (-)  1.2 1.0 0.8 0.6  50/50 25/75 0/100  0.4 0.2 0.0 X10  X100  X1000  Cases Figure 8.8: Effects of liquid generation rate on the gas velocity ratio.  8.2.3. Effects of GDL surface wettability on flow maldistribution In previous chapters we have found that the GDL and channel wall wettability have great impacts on two-phase flow patterns, flow maldistribution and fuel cell performances. To further study their impact on flow maldistribution with reaction, the effects of GDL surface wettability was firstly simulated in this section by varying the GDL surface static contact angle from 45o to 120o, i.e., 45o, 60o, 90o, and 120o. The case with an operating current density of 600 mA/cm2 and 25/75 case was selected for simulation since it represents a critical operating point according to the polarization curve in Figure 8.3. As shown in Figure 8.9, using more hydrophobic GDL can effectively expel liquid water from the GDL surface, and thus reduces the GDL water coverage ratio, which is of beneficial to the fuel cell performance because of the better liquid water removal from the GDL surface as observed in other chapters.  175  0.75  0.60  Voltage (V)  0.6  0.45 0.4  GDL water coverage ratio Voltage  0.30  0.2 0.15  0.0  GDL Water Coverage Ratio (-)  0.8  0.00 45  60  90  120  140 o  GDL Contact Angles( ) Figure 8.9: Effects of GDL surface contact angle on the fuel cell output voltage and the GDL water coverage ratio.  As shown in Figure 8.10, using either extremely hydrophilic or hydrophobic GDL surface can slightly reduce the pressure drop, which is consistent with the observation reported in Chapter 5. However, the total pressure drop over the channel with different GDL surface contact angles does not change much since the two-phase flow pattern in the channel remains in the slug flow. Using more hydrophobic GDL also mitigates the gas flow maldistribution in parallel channels because of a less flooded channel B. It can thus be concluded that a more hydrophobic GDL is recommended to be employed in the fuel cell in order to improve the fuel cell performance, reduce total pressure drop and mitigate the gas flow maldistribution.  176  0.4  Gas velocity ratio Pressure drop  150  0.3  100  0.2  50  0.1  0  0.0 45  60  90  120  Gas Velocity Ratio (B to A)  Pressure Drop (Pa)  200  140 o  GDL Contact Angles( ) Figure 8.10: Effects of GDL surface contact angle on the total pressure drop and the gas velocity ratio.  8.2.4. Effects of channel wall wettability on flow maldistribution Furthermore, the effects of channel wall wettability on flow maldistribution with reaction were simulated by varying the channel wall contact angle from 60o to 140o, i.e., 60o, 90o, 120o, and 140o, meanwhile keeping the GDL surface contact angle at 140o.Similarly to the last section, the operating current density was set at 600 mA/cm2, and only the 25/75 case was simulated. Changing channel wall wettability has less significant impact on the fuel cell performance than changing the GDL surface wettability. As shown in Figure 8.11, increasing the channel wall contact angle suppresses the expelling of liquid water from GDL surface to side walls, thus slightly increases the GDL water coverage ratio. However, on the other hand, the hydrophobic channel walls also reduce the residence time of liquid slugs in the channel. As a result, further increasing the channel wall contact angle leads to a small decrease in the GDL water coverage ratio. But in general, the fuel cell performance does not change too much while varying the channel wall wettability.  177  0.30 0.25  Voltage (V)  0.6 0.20 0.4  0.15  GDL water coverage ratio Voltage  0.10  0.2 0.05 0.0  GDL Water Coverage Ratio (-)  0.8  0.00 45  60  90  120  140 o  Channel Wall Contact Angles( ) Figure 8.11: Effects of channel wall contact angle on the fuel cell output voltage and the GDL water coverage ratio.  Since the slug flow is the dominating two-phase flow pattern in the channel, accelerating the liquid slug‟s removal will effectively reduce the total pressure drop. Therefore, as shown in Figure 8.12, increasing channel wall contact angle decreases the pressure drop. The use of more hydrophobic channel walls is also found to create less severe gas flow maldistribution in the cell. This is because that hydrophobic channel walls increase the number of liquid slugs in branch A, thus balance the flow resistance in the two branches. In conclusion, when the slug flow pattern is dominant in the parallel channels, hydrophobic channel walls give lower pressure drops and less severe gas flow maldistribution than hydrophilic walls, while the fuel cell performance remains unchanged.  178  0.8  Gas velocity ratio Pressure drop  150  0.6  100  0.4  50  0.2  0  Gas Velocity Ratio (B to A) (-)  Pressure Drop (Pa)  200  0.0 45  60  90  120  140 o  Channel Wall Contact Angles( ) Figure 8.12: Effects of channel wall contact angle on the total pressure drop and the gas velocity ratio.  8.2.5. Effects of cathode side gas stoichiometric flow ratio on flow maldistribution In Chapter 7, it is found that increasing the gas stoichiometric flow ratio can accelerate the removal of liquid slugs from the channel and thus increase the fuel cell performance. In this section, several gas stoichiometric flow ratios, i.e. 2, 3, 4, and 5, were used in the simulation at a constant current density of 600 mA/cm2 and under different degrees of flow maldistribution to investigate its impact on the fuel cell performance and gas flow distribution. Figure 8.13 shows the output voltage of the fuel cell when operated at different gas stoichiometric flow ratios. It is seen that under the severe liquid flow maldistribution condition, namely the 0/100 case, increasing the cathode gas flow rate does not improve the fuel cell performance. Because once branch B is completely flooded, the additional gas input will only flow through the other branch. As a result, the total effective mass transfer and reaction area will remain the same, which is supported by the GDL water coverage ratio results shown in Figure 8.14.  179  0.8  Voltage (V)  0.6  50/50 25/75 0/100  0.4  0.2  0.0 2  3  4  5  Stoich. Number Figure 8.13: Effects of gas stoichiometric flow ratio on the fuel cell output voltage.  At severe liquid flow maldistribution, the GDL water coverage ratio remains at about 0.5, indicating a completely flooded branch B. However, when there is no flow maldistribution or only a moderate flow maldistribution, increasing the gas flow rate can effectively lower the GDL water coverage ratio and improve the fuel cell performance. It is worth to be noted that although there is little noticeable voltage increase with increasing the gas flow rate for the 50/50 and 25/75 cases in Figure 8.13, the large drop in GDL water coverage ratio can greatly extend the fuel cell ohmic region, so that the fuel cell can be operated at a higher current density.  180  GDL Water Coverage Ratio(-)  0.6 0.5 0.4  50/50 25/75 0/100  0.3 0.2 0.1 0.0 2  3  4  5  Stoich. Number Figure 8.14: Effects of gas stoichiometric flow ratio on the GDL water coverage ratio.  Since increasing the gas flow rate cannot flush out the liquid slug in branch B when it is already fully flooded, the gas flow maldistribution in the 0/100 case does not get improved with increasing the gas flow rate (Figure 8.15). For the 25/75 case, the gas flow maldistribution becomes slightly less severe when the cell is operated at higher gas stoichiometric flow ratio due to the shift of the twophase flow pattern from slug flow to top wall droplet flow. However, increasing the gas stoichiometric flow ratio seems not to be an effective method for mitigating the flow maldistribution.  181  Gas Velocity Ratio (B to A) (-)  1.2 1.0 0.8  50/50 25/75 0/100  0.6 0.4 0.2 0.0 2  3  4  5  Stoich. Number Figure 8.15: Effects of gas stoichiometric flow ratio on the gas velocity ratio.  8.2.6. Effects of inlet resistances on flow maldistribution Gas flow restrictor/distributor is commonly used in multiphase flow systems to improve the gas and liquid flow distribution and to prevent the gas maldistribution caused by the difference in flow resistances in multiple channels. Bi et al. [149] theoretically analyzed the impact of gas distributor on the gas-liquid two-phase flow maldistibution in two parallel straight channels. They found that the creation of a pressure drop over a gas inlet distributor at an order of 30–50% of the pressure drop over the flow channel might be beneficial for suppressing the gas flow maldistribution. In this section, the same idea was examined by adding an orifice at the gas inlet of each branch. The pressure drop over the orifice is estimated by  1 U p   g  g 2  Cd    1  4   2         2  (8.24)  182  where  is the resistance coefficient defined as the ratio of the orifice diameter to the channel diameter, U g is the superficial gas velocity passing through the orifice,  g is the gas density, and Cd is the frictional coefficient with a value of 0.68 [149]. To investigate the impact of inlet resistance, simulations with several values of resistance coefficient, i.e., 1.0, 0.4, 0.15 and 0.1, were conducted with results being compared. Here a higher resistance coefficient represents a lower inlet resistance, and there is essentially no inlet resistance when its value is 1. The operating current density was set at 600 mA/cm2, and only the 0/100 and 25/75 cases were simulated since it is not necessary to add an orifice when the flow is uniformly distributed among the channels. As shown by the fuel cell voltage output in Figure 8.16, even under severe liquid flow maldistribution, namely the 0/100 case, the fuel cell performance can be recovered significantly when a high flow resistance orifice is used (i.e. at the resistance coefficient of 0.1 and 0.15). The small voltage difference between the 25/75 case and the 0/100 case with high inlet resistances suggests that the cell is operated in the ohmic region.  183  0.7 0.6  Voltage (V)  0.5 0.4  25/75 0/100  0.3 0.2 0.1 0.0 0.1  0.15  0.4  1  Resistance Coefficient (-) Figure 8.16: Effects of inlet resistance on the fuel cell output voltage.  Apparently, adding an inlet flow resistance to the flow channel will add an extra pressure loss to the overall pressure drop, as observed in Figure 8.17. However, under the severe flow maldistribution (the 0/100 case), the total pressure drop remains almost the same when the resistance coefficient is equal to or greater than 0.15. This is because that without the inlet resistance, branch B is fully flooded, and all the gas flows into branch A. The introduction of an inlet resistance decreases the gas flow into branch A, thus decreases the pressure drop over the branch A. If the reduced pressure drop over the branch straight section is balanced by the increased pressure loss over the inlet orifice, the total pressure drop can remain the same.  184  Pressure Drop (Pa)  400  25/75 0/100  300  200  100  0 0.1  0.15  0.4  1  Resistance coefficient (-) Figure 8.17: Effects of gas inlet resistance on the total pressure drop.  Adding gas inlet resistance can also effectively decrease the GDL water coverage ratio even under severe flow maldistribution conditions as shown in Figure 8.18, and the sharp decrease of the GDL water coverage ratio, as discussed in the last section, also extends the ohmic region. However, when the gas flow resistance is high enough, further increase in the resistance can no longer lower the GDL water coverage ratio, Instead, it will cause huge increases in the total pressure drop, as shown in Figure 8.17.  185  GDL Water Coverage Ratio (-)  0.5  25/75 0/100  0.4  0.3  0.2  0.1  0.0 0.1  0.15  0.4  1  Resistance Coefficient (-) Figure 8.18: Effects of gas inlet resistance on the GDL water coverage ratio.  As anticipated, the existence of an inlet resistance significantly reduces the gas flow maldistribution as shown in Figure 8.19. For a resistance coefficient of 0.1, nearly uniform gas distribution is achieved, which proves that adding a gas inlet resistance is a more effective maldistribution mitigation method than changing the wall and GDL material wettability or increasing gas flow rate.  186  Gas Velocity Ratio (B to A) (-)  1.0  0.8  0.6  25/75 0/100  0.4  0.2  0.0 0.1  0.15  0.4  1  Resistance Coefficient (-) Figure 8.19: Effects of inlet resistance on the gas velocity ratio.  8.3. Conclusions The following conclusions can be drawn from simulations presented in this chapter: 1) Slug flow patterns are mainly observed in the flow channel under different flow maldistribution conditions and at different current densities. More severe liquid flow maldistribution results in longer liquid slugs and the liquid slug can eventually occupy the entire branch, leading to significant voltage decrease in both the ohmic and the mass transport limited regions. There are no clear pressure drop differences caused by liquid flow maldistribution at the same current density. However, the flow maldistribution significantly increases the GDL water coverage ratio over the whole range of simulated current densities, which directly leads to poor fuel cell performance.  187  2) Even at low liquid generation rates, the flow maldistribution still can cause poor fuel cell performance and severe non-uniform gas distribution. As a result, liquid and gas flow maldistribution in parallel channels should be completely avoided if possible over the whole range of operations. 3) Using more hydrophobic GDL surfaces is recommended in order to improve the fuel cell performance, reduce total pressure drop and mitigate the gas flow maldistribution.  Meanwhile,  using  more  hydrophobic  channel  wall  is  recommended since it can lower the pressure drop and improve the gas flow maldistribution, while keeping the fuel cell performance unchanged. 4) Increasing the gas stoichiometric flow ratio is not an effective method to mitigate the gas flow maldistribution, although it can extend the fuel cell ohmic region, which makes the fuel cell capable of being operated at a higher current density. 5) Adding a gas inlet resistance is found to be the most effective maldistribution mitigation method among those evaluated in this study including changing the wall and GDL material wettability and increasing gas flow rate. With a carefully selected value of flow resistance coefficient, for example 0.15 for the simulated fuel cell in this chapter, both the fuel cell performance and the gas flow distribution can be significantly improved. Meanwhile, the total pressure drop can still remain the same, since the extra pressure loss caused by the inlet resistance can be cancelled out by the pressure drop decrease due to a more uniform gas distribution among channels.  188  9. Experimental observation and evaluation All the conclusions presented in previous chapters are based on the simulation results from a gas-liquid two-phase flow and electrochemical reaction model using the VOF method. To evaluate some of those simulation results, we have conducted three experiments. In the first experiment, the droplet behaviour and two-phase flow patterns were observed in a transparent mini-channel. Both the top and side views of the two-phase flow inside the channels were captured and used to qualitatively evaluate the three stages of droplet development as observed in the simulation as presented in Chapter 3 and the two-phase flow pattern evolvement as discussed in Chapter 4. The second experiment examined the gas flow maldistribution induced by injecting water into two parallel channels with different flow rates. The total pressure drop and the gas flow ratio were compared with the simulation results. Furthermore, two parallel flow channels with different sizes of communication channels were tested in the third experiment. Since it is difficult to track the gas or liquid communication in those channels, we only recorded the total pressure drop and the gas flow ratio at the entrance section of the channel before the first communication channel.  9.1. Droplet development in a transparent mini-channel In Chapter 3, three stages of droplet development were identified when water was injected through a GDL into the gas channel, namely emergence of liquid droplets on the GDL surface, droplet accumulation on the side walls, and droplet detachment from the top wall. In Chapter 4, we found that as the liquid injection rate increased, the two-phase flow pattern evolved from the top wall corner flow to top wall film flow, annular flow, and finally slug flow. To confirm these simulation-based findings, droplet behaviour and two-phase flow patterns were observed experimentally by a CCD camera in a transparent mini-channel.  189  Figure 2.1 shows the picture and schematic of the Plexiglas mini-channel. The length from gas inlet to outlet is 100 mm with a 1.0 mm × 1.0 mm square cross section. Toray® TGPH-030 GDL with 60% PTFE content was placed between the water reservoir and flow channel.  Figure 9.1: Picture and schematic of the transparent mini-channel.  Compressed air flowed into the channel from one end and liquid water was injected into the reservoir which then flowed through the GDL into the gas channel. Gas flow rate was set at 0.2 SLPM and liquid injecting rate was set at 200 mL/h, which correspond to 100 times of theoretical liquid generation rate at a  190  current density of 0.5 A∙cm-2. The contact angles of Plexiglas and the GDL are about 70o and 110o [161], respectively. Two CCD cameras were located above and by the side of the channel to take pictures of both top and side views of the channel. The dimension of the channel and operating conditions are the same as those used in the simulations in Chapter 4, which enables us to evaluate the simulation results both qualitatively and quantitatively. Water emerging method also mimics that in practical PEM fuel cells, and the transparent channel design can provide different views of two-phase flow in the channel, which helps to better understand droplet behaviors and flow regime transition in a real fuel cell flow channel. Figure 9.2 shows a typical droplet development on the GDL surface from side views of the channel.  Figure 9.2: Side views of liquid droplet development.  191  Similar to the simulation results shown in Chapter 3, a liquid droplet first emerges from the GDL surface and grows in size with time. Once it hits the top wall, due to the wettability difference between the GDL surface and channel walls, the droplet rapidly spreads out on the top wall, and detaches itself from the GDL surface. During this process, the droplet does not move along the channel, which implies that all liquid water flows downstream along the channel walls rather than along the GDL surface. The accumulation of droplets on the side walls as identified in the simulation was not clearly observed from the experiment, because once the droplet hit the side wall, it has already been big enough to reach the top wall immediately. Liquid droplets usually emerge from large pores of a real GDL. Unlike in the simulation presented in Chapter 3, the distance between adjacent droplets is too far for droplets merging on the GDL surface to collide, and the rough GDL surface also prevents the liquid droplet from moving on the GDL surface. Therefore, the droplets only merge together on the side or top walls. Multiple droplet behaviors on the GDL surface simulated in Chapter 3 may not accurately represent what happens with a real GDL, although the simulation results capture the single droplet behaviors in the channel. Figure 9.3 shows the snapshots of two-phase flow patterns as observed from the top of the channel. Several droplets can be identified on the GDL surface, but they do not move along the channel until they touch the side or top channel wall. At the beginning, the two-phase flow pattern is in the top wall corner flow (the first three pictures). As more droplets detach themselves from the GDL surface to the channel wall, the top wall film is eventually formed in the channel as shown in the last picture. Slug flow was not observed in the experiment at the given operating conditions, but it is most likely to appear in a longer channel.  192  Figure 9.3: Top views of two-phase flow patterns.  9.2. Two-phase flow maldistribution in parallel channels In Chapter 5, the impact of gas flow maldistribution on two-phase flow pattern, pressure drop and water distribution in two parallel channels was investigated. To evaluate the simulation results, similar two parallel channels were designed and constructed as shown in Figure 9.4.  193  Figure 9.4: Picture and schematic of two parallel channels.  The total channel length in Figure 9.4 is 30 cm, and the section with two parallel channels is 10 cm long with a 1.0 mm × 1.0 mm square cross section for each branch. The distance between two parallel channels is 1.0 mm. The channel plate was made of Plexiglas, and a hydrophobic Teflon tape was placed on the bottom to prevent the inherent communication through the GDL. Air flows into the channel from the left end at a velocity of 4.464 m/s. Liquid water was injected into each branch separately with a total water flow rate of 129.6 mL/h. Different liquid injection rate ratios between the two parallel channels were tested by changing the water flow rate in one branch while the total water flow rate was kept at a constant. All the operating conditions are the same and the channel dimension is similar as the cases simulated in Chapter 5, which makes sure the simulation and experimental results are comparable. The transparent top wall provides a top view of two-phase flow patterns in the parallel channels, which can be further analyzed to estimate the water volume fraction in each branch. Three pressure transducers were used to record the total pressure drop (from the air inlet to the outlet) and the pressure drops PL and PR over the entrance sections of the  194  bifurcating channels as shown in Figure 9.4. Since only air flows through the entrance section, the pressure drop, PL and PR, is proportional to the air velocity only. Because the two branches has exactly the same dimensions and are composed of the same materials, the pressure drop ratio is expected to be equal to the gas velocity ratio for laminar flow in the mini gas flow channels, which enables us to obtain the gas flow rate ratio and then compare it with simulation results. The total pressure drop P in the experiment is much higher than the one obtained from simulation due to the longer channel used in the experiment. To properly compare the two parameters, the pressure drop was normalized by the channel length. Figure 9.5 compares the gas flow rate ratios between experimental and simulation results.  Gas Velocity Ratio (-)  10 8  Sim. Exp.  6 4 2 0  0.0  0.2  0.4  0.6  0.8  1.0  Liquid Flow Rate Ratio (-) Figure 9.5: Effects of the degree of flow maldistribution on gas flow rate ratio.  Both results show a similar trend, i.e., when liquid flow rate ratio is higher than 0.25, the severity of gas maldistribution gradually increases with decreasing the  195  liquid flow ratio which can be mainly attributed to the uneven liquid volume fractions in the two channels. When the liquid flow ratio is lower than 0.25 with severe non-uniform distribution, gas flow maldistibution increases more substantially with decreasing the liquid flow ratio due to the appearance of the slug flow in the right branch as explained in Chapter 5. There are several factors that may cause the difference between experimental and simulated values. One possible factor is the wall roughness which was not considered in the simulation. The rough wall surface increases the liquid residence time, which may alter the two-phase flow patterns in both branches. Another possible cause is the difference in the material contact angle. In the experiment, the contact angles of Plexiglas and the Teflon tape are about 70o and 110o [161], respectively, while in the simulation the wall contact angle and the GDL contact angle were set at 45o and 140o, respectively. As we learnt in Chapter 5, the material wettability may also slightly influence the flow distribution. Finally, the inherent channel geometry variation due to imperfect manufacturing may also induce flow maldistribution even at a uniform water injection rate into the two parallel channels. The measured and simulated normalized total pressure drops are shown in Figure 9.6. Obviously, a general trend is followed by both measured and simulated pressure drops as the liquid flow ratio varies. From those results, we can conclude that the general characteristics of two-phase flow distribution in two parallel channels can be well captured from the simulation using the VOF method.  196  5000  P/L (Pa/m)  4500  4000  Sim. Exp.  3500  3000  0.0  0.2  0.4  0.6  0.8  1.0  Liquid Flow Rate Ratio (-) Figure 9.6: Effects of the degree of flow maldistribution on normalized pressure drop.  9.3. Effects of communication channels In Chapter 6, we found that communication between parallel channels can have a great impact on the two-phase flow pattern, gas and water distribution along the  channels  and  flow  maldistribution  between  the  channels.  Wide  communication channels provide a pathway for gas to avoid passing and contacting the liquid, which thus worsens the flow maldistribution between the parallel channels. However, narrow communication channels can prevent the gas exchange between parallel channels, which are helpful for mitigating the flow maldistribution. To evaluate above conclusions, two designs of parallel channels with different size of communication channels were tested in the experiment. One design has 5 communication channels, which are 1 mm in width and 1 mm in depth, evenly distributed along the parallel channels. Another design has 17 evenly distributed communication channels along the parallel channels, which are 0.25 mm in depth and 1 mm in width. The liquid flow ratio was set at 0.053 for all the tests, which is the same as simulated in Chapter 6.  197  Unlike the simulations, it is difficult to track the gas or liquid communication through each communication channel in the experiment. Comparing the PL and PR only yields the gas flow rate ratios before the first communication channel which are shown in Figure 9.7. From the simulation, it was found that without communication channels, the gas flow rate ratio was the highest since there was no gas or liquid exchange throughout the parallel channels. The configuration with 1 mm communication channels showed the least gas flow ratio due to the gas exchange at each communication channel. The same tendency is also identified from the experimental results, which implies that similar phenomena as predicted by the simulation may have occurred in the experiment.  Gas Flow Rate Ratio (-)  8  Sim. Exp.  6  4  2  0  No Comm.  0.25mm  1mm  Cases Figure 9.7: Effects of communication channels on gas flow rate ratio between parallel channels.  The normalized pressure drop is another parameter which can be used to compare the simulated and measured two-phase flow patterns inside the parallel channels. From the simulation, it was found the 1 mm communication channels case had the highest total pressure drop since liquid slugs were easily trapped in one branch after the first communication channel. Meanwhile, the 0.25 mm case  198  shows the lowest pressure drop due to the more uniform gas and liquid distribution. The experimental results, as shown in Figure 9.8, also reveal the same trend as predicted by the simulation, which to some extent evaluate s the simulation results.  P/L (Pa/m)  6000  5000  4000  Sim. Exp.  3000  2000  No Comm.  0.25mm  1mm  Cases Figure 9.8: Effects of communication channels on normalized pressure drops.  9.4. Conclusions The following conclusions can be drawn from the experimental work presented in this chapter. 1) The simulated droplet development and two-phase flow patterns were qualitatively evaluated in a transparent mini-channel. Two of the three stages of droplet development, namely droplet emergence from the GDL surface and detachment from the top wall were clearly observed in the experiment and the two main flow patterns, namely top wall corner flow and top wall film flow were identified in the experiment.  199  2) The effects of non-uniform liquid distribution on the gas flow maldistribution in two parallel channels were quantitatively evaluated by comparing the gas flow ratios and the total pressure drop under different liquid injection rate ratios. The experimental results confirm that the characteristics of two-phase flow maldistribution in two parallel channels can be well captured from the simulation. 3) Two designs of communicating channels were examined in experiments to evaluate the effectiveness of communication channels on improving the flow distribution. By comparing the total pressure drop and the gas flow ratios at the entrance section of the channel, it was found that the effectiveness of the communication channels strongly depended on the size of the communication channels, which is similar to what was predicted from simulations and the results to some extent validate the simulation results.  200  10. Conclusions and recommendations 10.1. Conclusions This study focused on understanding the water management issues in gas flow channels of PEM fuel cells. Such issues include droplet behaviour, two-phase flow pattern development, impact of two-phase flow on PEM fuel cell performance, flow maldistribution in multiple parallel channels and its impact on the cell performance, and mitigation of flow maldistribution. The following conclusions were drawn from each chapter of this thesis. 10.1.1. Droplet behaviour Three stages of droplet behaviour in a typical PEM fuel cell flow channel can be identified when water is injected through the GDL into the gas channel, namely emergence of water droplet on the GDL surface, accumulation on the side walls, and detachment from the top wall. These three stages were observed in almost all the cases both numerically and experimentally. Changing gas or liquid flow rate does not alter the three stages of droplet behaviour. However, the microstructure of the GDL surface was found to have a significant impact on the multiple droplet behaviour in a mini-channel. The uniform liquid inlet case and 1-pore case showed distinctly different flow patterns. However, for the 4-pore, 16-pore and 64-pore cases, similar flow patterns are observed, suggesting that a minimum of 4 pores may be required to represent the microstructure of the GDL surface in a fuel cell mini-channel. However, from the experiment observation, it was found that the distance between each emerging pore may be too far for droplets to coalesce on the GDL surface. Instead, they only merge together on the side or top walls. Multiple droplet behaviors on the GDL surface in simulation may not accurately represent the case with a real GDL, but the simulation properly captures single droplet behaviours in a gas flow channel.  201  10.1.2. Two-phase flow pattern development In a large-scale gas flow channel, at theoretical liquid water production rate, the two-phase flow shows the corner droplet flow pattern, with liquid water mainly present in the bottom corners. With increasing liquid injection rates or over a long gas flow channel, the flow pattern may evolve from the corner droplet flow to corner flow, top wall film flow, annular flow, and eventually slug flow. The two main flow patterns, namely the corner flow and top wall film flow, were also observed in the experiment, which qualitatively evaluated the capability of the VOF method in capturing the gas-liquid flow in fuel cell flow channels. 10.1.3. Impact of two-phase flow on PEM fuel cell performance Under severe flooding, e.g., at the end section of a long flow channel or after an extended period of operation, the slug flow will affect fuel cell voltage output in the mass transport limited region at high current densities, causing a large cell voltage decrease and fluctuations, which are detrimental to the fuel cell stability and durability. The slug flow also creates high pressure drops and large pressure fluctuations, causing extra energy loss. Therefore, it should be always avoided in operating PEM fuel cells. Increasing gas stoichiometric flow ratios can effectively extend the stable fuel cell operation range by altering the two-phase flow pattern in the channel. However, when operating at a certain current density, too large a stoichiometric ratio is not recommended, because it may cause more parasitic power loss with little gain in voltage output. 10.1.4. Flow maldistribution and its impact on the cell performance Severe liquid flow maldistribution caused by uneven liquid flow rates in two parallel channels may transform the two-phase flow patterns from the top wall droplet flow to the slug flow in the branch that has a higher liquid flow rate. Severe uneven gas distribution also increases the pressure drop, gas velocity fluctuations and the GDL surface liquid coverage.  202  Varying the degree of liquid flow maldistribution does not change the two-phase flow patterns except under severe maldistribution conditions, i.e. at a liquid flow rate ratio less than 0.25, where slug flow starts to appear, gas maldistibution becomes more severe, and GDL surface liquid coverage soars up. Pressure drops at different degrees of liquid flow maldistribution show a wavy shape that may be associated with the flow instabilities in parallel channels. The effects of non-uniform liquid distribution on the gas flow maldistribution in two parallel channels were quantitatively evaluated by comparing the measured and predicted gas flow ratios and the total pressure drop under different liquid injection rate ratios, and it is proven that the characteristics of two-phase flow maldistribution in two parallel channels can be well captured by the current simulation. Different flow resistances are caused by inherent geometric differences, GDL intrusion during compression [114], uneven liquid water accumulation in the flow field channels [115], or two-phase flow instabilities [126]. Our simulation results also showed that the GDL microstructure was probably also responsible for the flow maldistribution due to the different flow resistance it causes. In two parallel channels with electrochemical reactions, slug flow is the dominant two-phase flow pattern observed. More severe flow maldistribution results in longer liquid slugs and the liquid slug can eventually occupy the entire branch, leading to a substantial voltage decrease in both the ohmic and the mass transport limited region. There is little change in the pressure drop when there is a flow maldistribution for the cell operated at the same current density. However, the flow maldistribution significantly increases the GDL water coverage ratio over the whole range of current densities, which may lead to poor fuel cell performance. Even at low liquid generation rates, the flow maldistribution can still cause poor fuel cell performance and severe non-uniform gas distribution.  203  10.1.5. Impact of material wettability Liquid on a hydrophobic GDL surface tends to form droplets, occupying a large cross sectional area of the channel and causing high flow resistance. But those droplets are easy to detach from the surface. On the other hand, liquid on a hydrophilic GDL surface tends to attach on it and form a thin liquid film, which has a longer residence time compared to the liquid on a hydrophobic surface. For those reasons, hydrophobic GDL and hydrophilic channel wall are commonly used in typical PEM fuel cells, aiming to reduce the parasitic energy loss and improve the fuel cell performance. However, our simulation results showed that these two objectives were always in conflict with each other at certain operating conditions or for a given channel configuration. Using a more hydrophobic GDL always reduces the GDL water coverage and thus improves the fuel cell performance. Usually the pressure drop also increases because more channel cross sectional area is occupied by the droplet on the GDL surface. One exception is that when slug flow appears in the channel, increasing the GDL hydrophobicity will accelerate the water removal which may reduce the pressure drop as shown in Section 7.2.5. In parallel channels, increasing the GDL hydrophobicity leads to the transformation of two-phase flow patterns from the GDL droplet flow or slug flow to the top wall droplet flow, thus reducing the total pressure drop. Depending on the two-phase flow pattern in the channel, more hydrophilic channel walls may be favourable for reducing the GDL water coverage and increasing the fuel cell output voltage in cases when there are no liquid slugs in the channel. Too hydrophilic a channel wall increases the liquid resident time and liquid on a hydrophobic channel wall occupy more cross sectional area, therefore liquid slugs are easily formed under those two circumstances, thus increasing the pressure drop. However, further increase in the channel wall hydrophobicity facilitates the slug‟s removal. As a result, the pressure drop showed a wave shape with varying the channel wall wettability, as shown in Chapters 3 and 5. 204  10.1.6. Maldistribution mitigation methods Increasing the gas stoichiometric flow ratio is not an effective method to mitigate the gas flow maldistribution, since the additional gas input will only flow through the branch with less water. However, it can extend the fuel cell ohmic region, which makes the fuel cell capable of being operated at a higher current density. Using a more hydrophobic GDL surface or more hydrophobic channel walls can slightly improve the gas flow maldistribution, but the effect is quite limited. Using narrow communication channels can effectively mitigate gas flow maldistribution by redistributing the liquid among the parallel flow channels. However, when the communication channels are not narrow enough, they will provide a pathway for the gas to short-circuit the liquid, leading to a poorer gas flow distribution. One side-effect of using communication channels is the delayed water removal from flow channels, which would decrease the PEM fuel cell performance. Wide communication channels even make water permanently remain in the channel, leading to higher total pressure drops. Communication channels should not be applied to serpentine channels, since adding communication channels will lead to gas short circuit and thus lower the reactant conversion. Adding gas inlet resistances is found to be the most effective method in mitigation maldistribution. With a carefully selected value of resistance coefficient both the fuel cell performance and gas flow maldistribution can be significantly improved. Meanwhile, the total pressure drop still remains more or less the same, since the extra pressure loss caused by the gas inlet resistance can be cancelled out by the pressure drop decrease due to a more uniform gas distribution.  205  10.2. Suggestions for ideal gas flow channel design Based on simulation results in this work, several suggestions are given for better flow channel designs,   More hydrophobic GDL surface should be used to expel liquid water from GDL surface. Depending on the amount of liquid water in the flow channel, more hydrophilic channel wall is preferred in order to form liquid film flow on the top wall, but when slug flow is inevitable, more hydrophobic channel wall is suggested to accelerate the slug removal.    Droplets in a channel with lower height to width ratio or a trapezoidal channel may easily detach on the top wall rather than on the GDL surface, which helps to form top wall film flow.    When using communicating channels in multiple parallel channels, it is recommended that the width of communicating channel should be less than 0.5 mm.    Adding gas distributor at the inlet of parallel channels is highly recommended to mitigate the flow maldistribution, and the value of resistance coefficient should be selected based on the overall system pressure drop and the degree of the flow maldistribution.  10.3. Recommendations for future work This thesis investigated several water management issues in PEM fuel cells based on CFD simulations. It is worth noting that all the conclusions obtained from simulation are only validated in limited operating range. The following recommendations are suggested for modeling development, experimental validation, and further simulation/experiment work.  206  10.3.1. Modelling development   Couple contact angle hysteresis model with current VOF method. With contact angle hysteresis, droplets tend to stay longer on a surface, which may change the two-phase flow pattern in the channel. For example, slug flow is more easily formed, which may affect the pressure drop, flow maldistribution and fuel cell performances.    Add wall roughness model into the VOF model, since rough wall delays droplet removal, the slug flow is expected to be more easily formed.    Add the gas diffusion layer in the simulation and consider its microstructure, so that the gas and liquid transport inside GDL can be well captured. The calculation of fuel cell output voltage will become more accurate since the gas reactant concentration at the catalyst layer can be explicitly accessed. Flow maldistribution caused by non-uniform water transport inside GDL and gas/liquid communication between parallel channels though GDL can be simulated.    Add membrane layer and anode side channel into the current 1+3D model, so that the water crossover can be calculated.    Add heat transport into the 1+3D model to investigate the heat management coupling with water management in PEM fuel cells.    Mesh refinement for more accurate results.  10.3.2. Experimental validation   Two-phase flow in a channel using materials with different wettability should be performed. In addition to flow pattern visualization, it would be  207  useful if quantitative variables such as water holdup, water coverage ratio on different surfaces can be obtained.   In-situ experiments should be performed to validate the 1+3D model developed in this thesis.  10.3.3. Further simulation/experiment work   In order to validate the 1+3D model, reaction parameters should be first obtained from experiment with the cell operated in single-phase flow. Then, the polarization curve and total pressure drop from experiments under two-phase flow conditions will be compared with simulation results.    Two-phase flow patterns and performance of PEM fuel cells with different channel designs, such as serpentine or interdigitated channels, can be simulated using the 1+3D model in the future study.    More simulations should be performed in large-scale fuel cells or fuel cell stacks with more channels, where liquid flooding and flow maldistribution are expected to be more severe issues.    More parametric studies should be conducted in the simulation. Such variables include gas humidity, operating pressure and temperature.    It would be interesting to investigate other factors that may lead to flow maldistribution using VOF method, e.g., GDL intrusion and channel manifold design. Also other maldistribution mitigating methods, e.g., new channel designs, can be examined in simulation.    It is necessary to test the maldistribution mitigating methods proposed in this thesis in a real PEM fuel cell. Special attention should be given to the  208  two methods, namely using communication channels and adding inlet resistance, which could greatly reduce the flow maldistribution.   It is better to avoid liquid water‟s formation in the first place by having a more complicated humidity and heat management system, more simulation and experimental work should be carried out to design such a system.  209  References 1.  Barbir, F., PEM Fuel Cells: Theory and Practice2005: Academic Press.  2.  Li, H., et al., A review of water flooding issues in the proton exchange membrane fuel cell. Journal of Power Sources, 2008. 178(1): p. 103-117.  3.  Weber, A.Z. and J. Newman, Modeling transport in polymer-electrolyte fuel cells. Chemical Reviews, 2004. 104(10): p. 4679-4726.  4.  Wang, Y., C.-Y. Wang, and K.S. Chen, Elucidating differences between carbon paper and carbon cloth in polymer electrolyte fuel cells. Electrochimica Acta, 2007. 52(12): p. 3965-3975.  5.  Djilali, N. and P.C. Sui, Transport phenomena in fuel cells: from microscale to macroscale. International Journal of Computational Fluid Dynamics, 2008. 22(1-2): p. 115-133.  6.  Litster, S. and N. Djilali, Two-phase transport in porous gas diffusion electrodes, in Transport phenomena in fuel cells, B. Sunden and M. Faghri, Editors. 2005, WIT Press. p. 175-213.  7.  EG&G Technical Services, I., Fuel Cell Handbook. 7 ed2004: CRC Press.  8.  Berning, T., Three-dimensional computational analysis of transport phenomena in a PEM fuel cell, 2002, University of Victoria (Canada).  9.  Harvey, D.B.P., Three-dimensional CFD model for PEMFC cathodes: Application to serpentine flow fields, 2006, Queen's University at Kingston (Canada).  10.  Larminie, J. and A. Dicks, Fuel Cell Systems Explained, Second Edition. 2003: John Wiley and Sons.  11.  Trabold, T.A., Minichannels in Polymer Electrolyte Membrane Fuel Cells. Heat Transfer Engineering, 2005. 26(3): p. 3-12.  12.  Tüber, K., D. Pócza, and C. Hebling, Visualization of water buildup in the cathode of a transparent PEM fuel cell. Journal of Power Sources, 2003. 124(2): p. 403-414.  13.  Weng, F.-B., et al., Study of water-flooding behaviour in cathode channel of a transparent proton-exchange membrane fuel cell. Journal of Power Sources, 2006. 157(2): p. 674-680.  14.  Zhang, F.Y., X.G. Yang, and C.Y. Wang, Liquid Water Removal from a Polymer Electrolyte Fuel Cell. Journal of the Electrochemical Society, 2006. 153(2): p. A225-A232-A225-A232.  210  15.  Liu, X., et al., Flow dynamic characteristics in flow field of proton exchange membrane fuel cells. International Journal of Hydrogen Energy, 2008. 33(3): p. 1040-1051.  16.  Masuda, H., et al., Comparison between numerical simulation and visualization experiment on water behavior in single straight flow channel polymer electrolyte fuel cells. Journal of Power Sources, 2008. 177(2): p. 303-313.  17.  Spernjak, D., 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, 2007. 170(2): p. 334-344.  18.  Ge, S. and C.-Y. Wang, Liquid Water Formation and Transport in the PEFC Anode. Journal of the Electrochemical Society, 2007. 154(10): p. B998-B1005-B998-B1005.  19.  Ma, H.P., et al., Diagnostic tool to detect liquid water removal in the cathode channels of proton exchange membrane fuel cells. Journal of Power Sources, 2006. 162(1): p. 469-473.  20.  He, G., Y. Yamazaki, and A. Abudula, The effect of wall roughness on the liquid removal in micro-channels related to a proton exchange membrane fuel cell(PEMFC). Journal of Power Sources, 2010. 195(6): p. 1561-1568.  21.  Chen, Y.P., et al., Condensation in microchannels. Nanoscale and Microscale Thermophysical Engineering, 2008. 12(2): p. 117-143.  22.  Cho, J., H.-S. Kim, and K. Min, Transient response of a unit protonexchange membrane fuel cell under various operating conditions. Journal of Power Sources, 2008. 185(1): p. 118-128.  23.  Chang, P.A.C., et al., Flow distribution in proton exchange membrane fuel cell stacks. Journal of Power Sources, 2006. 162(1): p. 340-355.  24.  Geiger, A.B., et al., In Situ Investigation of Two-Phase Flow Patterns in Flow Fields of PEFCrsquos Using Neutron Radiography. Fuel Cells, 2002. 2(2): p. 92-98.  25.  Pekula, N., et al., Study of water distribution and 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, 2005. 542(1-3): p. 134-141.  26.  Hickner, M.A., et al., Real-Time Imaging of Liquid Water in an Operating Proton Exchange Membrane Fuel Cell. Journal of the Electrochemical Society, 2006. 153(5): p. A902-A908-A902-A908.  27.  Owejan, J.P., et al., In situ investigation of water transport in an operating PEM fuel cell using neutron radiography: Part 2 - Transient water 211  accumulation in an interdigitated cathode flow field. International Journal of Heat and Mass Transfer, 2006. 49(25-26): p. 4721-4731. 28.  Park, J., et al., Neutron imaging investigation of liquid water distribution in and the performance of a PEM fuel cell. International Journal of Hydrogen Energy, 2008. 33(13): p. 3373-3384.  29.  Trabold, T.A., et al., 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, 2006. 49(25-26): p. 4712-4720.  30.  Satija, R., et al., In situ neutron imaging technique for evaluation of water management systems in operating PEM fuel cells. Journal of Power Sources, 2004. 129(2): p. 238-245.  31.  Wang, Y. and K.S. Chen, Through-Plane Water Distribution in a Polymer Electrolyte Fuel Cell: Comparison of Numerical Prediction with Neutron Radiography Data. Journal of the Electrochemical Society, 2010. 157(12): p. B1878-B1886-B1878-B1886.  32.  Tsushima, S., 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, 2004. 7(9): p. A269-A272A269-A272.  33.  Tsushima, S., et al., Water content distribution in a polymer electrolyte membrane for advanced fuel cell system with liquid water supply. Magnetic Resonance Imaging, 2005. 23(2): p. 255-258.  34.  Minard, K.R., et al., Magnetic resonance imaging (MRI) of PEM dehydration and gas manifold flooding during continuous fuel cell operation. Journal of Power Sources, 2006. 161(2): p. 856-863.  35.  Dunbar, Z. and R.I. Masel, Quantitative MRI study of water distribution during operation of a PEM fuel cell using Teflon® flow fields. Journal of Power Sources, 2007. 171(2): p. 678-687.  36.  Bedet, J., et al., Magnetic resonance imaging of water distribution and production in a 6 cm2 PEMFC under operation. International Journal of Hydrogen Energy, 2008. 33(12): p. 3146-3149.  37.  Tsushima, S., et al., Investigation of Water Distribution in a Membrane in an Operating PEMFC by Environmental MRI. Journal of the Electrochemical Society, 2010. 157(12): p. B1814-B1818-B1814-B1818.  38.  Albertini, V.R., et al., In Situ XRD Studies of the Hydration Degree of the Polymeric Membrane in a Fuel Cell. Electrochemical and Solid-State Letters, 2004. 7(12): p. A519-A521-A519-A521.  212  39.  Manke, I., et al., Investigation of water evolution and transport in fuel cells with high resolution synchrotron x-ray radiography. Applied Physics Letters, 2007. 90(17): p. 174105-174105.  40.  Sinha, P.K., P.P. Mukherjee, and C.-Y. Wang, Impact of GDL structure and wettability on water management in polymer electrolyte fuel cells. Journal of Materials Chemistry, 2007. 17(30): p. 3089-3103.  41.  Lee, S.J., et al., X-ray imaging of water distribution in a polymer electrolyte fuel cell. Journal of Power Sources, 2008. 185(2): p. 867-870.  42.  Hartnig, C., et al., High-resolution in-plane investigation of the water evolution and transport in PEM fuel cells. Journal of Power Sources, 2009. 188(2): p. 468-474.  43.  Krüger, P., et al., Synchrotron X-ray tomography for investigations of water distribution in polymer electrolyte membrane fuel cells. Journal of Power Sources, 2011. 196(12): p. 5250-5255.  44.  Bazylak, A., Liquid water visualization in PEM fuel cells: A review. International Journal of Hydrogen Energy, 2009. 34(9): p. 3845-3857.  45.  Ji, M. and Z. Wei, A Review of Water Management in Polymer Electrolyte Membrane Fuel Cells. Energies, 2009. 2(4): p. 1057-1106.  46.  Jiao, K. and X. Li, Water transport in polymer electrolyte membrane fuel cells. Progress in Energy and Combustion Science, 2011. 37(3): p. 221291.  47.  St-Pierre, J., PEMFC In Situ Liquid-Water-Content Monitoring Status. Journal of the Electrochemical Society, 2007. 154(7): p. B724-B731-B724B731.  48.  Tsushima, S. and S. Hirai, In situ diagnostics for water transport in proton exchange membrane fuel cells. Progress in Energy and Combustion Science, 2011. 37(2): p. 204-220.  49.  Springer, T.E., Polymer Electrolyte Fuel Cell Model. Journal of the Electrochemical Society, 1991. 138(8): p. 2334-2334.  50.  Bernardi, D.M. and M.W. Verbrugge, Mathematical model of a gas diffusion electrode bonded to a polymer electrolyte. AIChE Journal, 1991. 37(8): p. 1151-1163.  51.  Weber, A.Z. and J. Newman, Modeling transport in polymer-electrolyte fuel cells. Chemical Reviews (Washington, DC, United States), 2004. 104(10): p. 4679-4726.  52.  Wang, C.-Y., Fundamental Models for Fuel Cell Engineering. Chemical Reviews (Washington, DC, United States), 2004. 104(10): p. 4727-4766.  213  53.  BiyIkoglu, A., Review of proton exchange membrane fuel cell models. International Journal of Hydrogen Energy, 2005. 30(11): p. 1181-1212.  54.  Tao, W.Q., et al., Parameter sensitivity examination and discussion of PEM fuel cell simulation model validation: Part I. Current status of modeling research and model development. Journal of Power Sources, 2006. 160(1): p. 359-373.  55.  Siegel, C., Review of computational heat and mass transfer modeling in polymer-electrolyte-membrane (PEM) fuel cells. Energy, 2008. 33(9): p. 1331-1352.  56.  Baschuk, J.J. and X. Li, Modelling of polymer electrolyte membrane fuel cells with variable degrees of water flooding. Journal of Power Sources, 2000. 86(1-2): p. 181-196.  57.  Natarajan, D. and T. Van Nguyen, A two-dimensional, two-phase, multicomponent, transient model for the cathode of a proton exchange membrane fuel cell using conventional gas distributors. Journal of the Electrochemical Society, 2001. 148(12): p. A1324-A1335-A1324-A1335.  58.  Berning, T. and N. Djilali, A 3D, Multiphase, Multicomponent Model of the Cathode and Anode of a PEM Fuel Cell. Journal of the Electrochemical Society, 2003. 150(12): p. A1589-A1598-A1589-A1598.  59.  Wang, Z.H., C.Y. Wang, and K.S. Chen, Two-phase flow and transport in the air cathode of proton exchange membrane fuel cells. Journal of Power Sources, 2001. 94(1): p. 40-50.  60.  Gurau, V., et al., A look at the multiphase mixture model for PEM fuel cell simulations. Electrochemical and Solid State Letters, 2008. 11(8): p. B132-B135-B132-B135.  61.  Wang, C.-Y., Comment on ``A Look at the Multiphase Mixture Model for PEM Fuel Cell Simulations'' [Electrochem. Solid-State Lett., [bold 11], B132 (2008)]. Electrochemical and Solid-State Letters, 2009. 12(2): p. S2S3-S2-S3.  62.  Gurau, V., Response to ``Comment on `A Look at the Multiphase Mixture Model for PEM Fuel Cell Simulations' '' [Electrochem. Solid-State Lett., [bold 11], B132 (2008)]. Electrochemical and Solid-State Letters, 2009. 12(2): p. S4-S6-S4-S6.  63.  Hirt, C.W. and B.D. Nichols, Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries. Journal of Computational Physics, 1981. 39: p. 201-225.  214  64.  Quan, P., et al., Water behavior in serpentine micro-channel for proton exchange membrane fuel cell cathode. Journal of Power Sources, 2005. 152: p. 131-145.  65.  Osher, S. and J.A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 1988. 79(1): p. 12-49.  66.  Osher, S. and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces2003.  67.  Unverdi, S.O. and G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows. Journal of Computational Physics, 1992. 100(1): p. 25-37.  68.  Tryggvason, G., et al., A Front-Tracking Method for the Computations of Multiphase Flow. Journal of Computational Physics, 2001. 169(2): p. 708759.  69.  Wörner, M., Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications. Microfluidics and Nanofluidics, 2012. 12(6): p. 841-886.  70.  Xia, Z., et al., Flow Patterns in the Sedimentation of an Elliptical Particle. Journal of Fluid Mechanics, 2009. 625: p. 249-272.  71.  Mehrabian, H., P. Gao, and J.J. Feng, Wicking flow through microchannels. Physics of Fluids, 2011. 23(12): p. 122108-122108-14122108-122108-14.  72.  Vanderlei, B., J.J. Feng, and L. Edelstein-Keshet, A Computational Model of Cell Polarization and Motility Coupling Mechanics and Biochemistry. Multiscale Modeling & Simulation, 2011. 9(4): p. 1420-1443.  73.  Vélez-Cordero, J.R., et al., Hydrodynamic interaction between a pair of bubbles ascending in shear-thinning inelastic fluids. Journal of NonNewtonian Fluid Mechanics, 2011. 166(1–2): p. 118-132.  74.  Chen, H., S. Chen, and W.H. Matthaeus, Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method. Physical Review A, 1992. 45(8): p. R5339-R5342-R5339-R5342.  75.  Le, A.D. and B. Zhou, A general model of proton exchange membrane fuel cell. Journal of Power Sources, 2008. 182(1): p. 197-222.  76.  Le, A.D. and B. Zhou, Fundamental understanding of liquid water effects on the performance of a PEMFC with serpentine-parallel channels. Electrochimica Acta, 2009. 54(8): p. 2137-2154.  215  77.  Kandlikar, S.G., Microscale and Macroscale Aspects of Water Management Challenges in PEM Fuel Cells. Heat Transfer Engineering, 2008. 29(7): p. 575-587.  78.  Ous, T. and C. Arcoumanis, Visualisation of water droplets during the operation of PEM fuel cells. Journal of Power Sources, 2007. 173(1): p. 137-148.  79.  Wang, M., H. Guo, and C. Ma, Temperature distribution on the MEA surface of a PEMFC with serpentine channel flow bed. Journal of Power Sources, 2006. 157(1): p. 181-187.  80.  Schillberg, C. and S. Kandlikar. in Proceedings of the Fifth International Conference on Nanochannels, Microchannels and Minichannels. June, 2007. Puebla, Mexico.  81.  Concus, P. and R. Finn, On the behavior of a capillary surface in a wedge. Appl. Math Sci. , 1969. 63(2): p. 292-299.  82.  Theodorakakos, A., et al., Dynamics of water droplets detached from porous surfaces of relevance to PEM fuel cells. Journal of Colloid and Interface Science, 2006. 300(2): p. 673-687.  83.  Shirani, E. and S. Masoomi, Deformation of a Droplet in a Channel Flow. Journal of Fuel Cell Science and Technology, 2008. 5(4): p. 041008-8.  84.  Kumbur, E.C., K.V. Sharp, and M.M. Mench, Liquid droplet behavior and instability in a polymer electrolyte fuel cell flow channel. Journal of Power Sources, 2006. 161(1): p. 333-345.  85.  Chen, K.S., M.A. Hickner, and D.R. Noble, Simplified models for predicting the onset of liquid water droplet instability at the gas diffusion layer/gas flow channel interface. International Journal of Energy Research, 2005. 29(12): p. 1113-1132.  86.  Chen, K.S., Predicting Water-Droplet Detachment from GDL/Channel Interfaces in PEM Fuel Cells. ECS Transactions, 2007. 11(1): p. 715-724.  87.  Zhu, X., P.C. Sui, and N. Djilali, Dynamic behaviour of liquid water emerging from a GDL pore into a PEMFC gas flow channel. Journal of Power Sources, 2007. 172(1): p. 287-295.  88.  He, G., et al., A two-fluid model for two-phase flow in PEMFCs. Journal of Power Sources, 2007. 163(2): p. 864-873.  89.  Hao, L. and P. Cheng, Lattice Boltzmann simulations of anisotropic permeabilities in carbon paper gas diffusion layers. Journal of Power Sources, 2009. 186(1): p. 104-114.  216  90.  Zhu, X., P.C. Sui, and N. Djilali, Numerical simulation of emergence of a water droplet from a pore into a microchannel gas stream. Microfluidics and Nanofluidics, 2007. 4(6): p. 543-555.  91.  Fang, C., et al., 3-D numerical simulation of contact angle hysteresis for microscale two phase flow. International Journal of Multiphase Flow, 2008. 34(7): p. 690-705.  92.  Cai, Y.H., et al., Effects of hydrophilic/hydrophobic properties on the water behavior in the micro-channels of a proton exchange membrane fuel cell. Journal of Power Sources, 2006. 161(2): p. 843-848.  93.  Kandlikar, S.G. and W.J. Grande, Evolution of Microchannel Flow Passages--Thermohydraulic Performance and Fabrication Technology. Heat Transfer Engineering, 2003. 24(1): p. 3-17.  94.  Dunbar, Z.W. and R.I. Masel, Magnetic resonance imaging investigation of water accumulation and transport in graphite flow fields in a polymer electrolyte membrane fuel cell: Do defects control transport? Journal of Power Sources, 2008. 182(1): p. 76-82.  95.  Taitel, Y., et al., Flow distribution of gas and liquid in parallel pipes. International Journal of Multiphase Flow, 2003. 29(7): p. 1193-1202.  96.  Kimball, E., et al., Drops, slugs, and flooding in polymer electrolyte membrane fuel cells. AIChE Journal, 2008. 54(5): p. 1313-1332.  97.  Murahashi, T., H. Kobayashi, and E. Nishiyama, Combined measurement of PEMFC performance decay and water droplet distribution under low humidity and high CO. Journal of Power Sources, 2008. 175(1): p. 98-105.  98.  Hussaini, I.S. and C.-Y. Wang, Visualization and quantification of cathode channel flooding in PEM fuel cells. Journal of Power Sources, 2009. 187(2): p. 444-451.  99.  Lu, Z., et al., 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, 2009. 34(8): p. 3445-3456.  100.  Zhu, X., P.C. Sui, and N. Djilali, Three-dimensional numerical simulations of water droplet dynamics in a PEMFC gas channel. Journal of Power Sources, 2008. 181(1): p. 101-115.  101.  Zhu, X., et al., Numerical investigation of water droplet dynamics in a lowtemperature fuel cell microchannel: Effect of channel geometry. Journal of Power Sources, 2010. 195(3): p. 801-812.  217  102.  Zhu, X., et al., Dynamics of Emerging Water Droplet Subjected to Sidewall with Different Wettabilities in a Fuel Cell Cathode Channel. Fuel Cells, 2011. 11(3): p. 404-412.  103.  Quan, P. and M.-C. Lai, Numerical study of water management in the air flow channel of a PEM fuel cell cathode. Journal of Power Sources, 2007. 164(1): p. 222-237.  104.  Quan, P. and M.-C. Lai, Numerical Simulation of Two-Phase Water Behavior in the Cathode of an Interdigitated Proton Exchange Membrane Fuel Cell. Journal of Fuel Cell Science and Technology, 2010. 7(1): p. 011017-14.  105.  Wang, X. and B. Zhou, Liquid water flooding process in proton exchange membrane fuel cell cathode with straight parallel channels and porous layer. Journal of Power Sources, 2011. 196(4): p. 1776-1794.  106.  Le, A.D. and B. Zhou, A generalized numerical model for liquid water in a proton exchange membrane fuel cell with interdigitated design. Journal of Power Sources, 2009. 193(2): p. 665-683.  107.  Akbar, M.K., D.A. Plummer, and S.M. Ghiaasiaan, On gas–liquid twophase flow regimes in microchannels. International Journal of Multiphase Flow, 2003. 29(5): p. 855-865.  108.  Al-Baghdadi, M.A.R.S. and H.A.K.S. Al-Janabi, Effect of operating parameters on the hygro–thermal stresses in proton exchange membranes of fuel cells. International Journal of Hydrogen Energy, 2007. 32(17): p. 4510-4522.  109.  Kim, N.-H. and T.-R. Sin, Two-phase flow distribution of air-water annular flow in a parallel flow heat exchanger. International Journal of Multiphase Flow, 2006. 32(12): p. 1340-1353.  110.  Marchitto, A., et al., Experiments on two-phase flow distribution inside parallel channels of compact heat exchangers. International Journal of Multiphase Flow, 2008. 34(2): p. 128-144.  111.  Delsman, E.R., et al., Microchannel Plate Geometry Optimization for Even Flow Distribution at High Flow Rates. Chemical Engineering Research and Design, 2004. 82(2): p. 267-273.  112.  Commenge, J.M., et al., Optimal design for flow uniformity in microchannel reactors. AIChE Journal, 2002. 48(2): p. 345-358.  113.  Kandlikar, S.G., et al., Measurement of flow maldistribution in parallel channels and its application to ex-situ and in-situ experiments in PEMFC water management studies. International Journal of Heat and Mass Transfer, 2009. 52(7-8): p. 1741-1752.  218  114.  Basu, S., J. Li, and C.-Y. Wang, Two-phase flow and maldistribution in gas channels of a polymer electrolyte fuel cell. Journal of Power Sources, 2009. 187(2): p. 431-443.  115.  Hamilton, P.J. and B.G. Pollet, Polymer Electrolyte Membrane Fuel Cell (PEMFC) Flow Field Plate: Design, Materials and Characterisation. Fuel Cells, 2010. 10(4): p. 489-509.  116.  Knights, S.D., et al., Aging mechanisms and lifetime of PEFC and DMFC. Journal of Power Sources, 2004. 127(1-2): p. 127-134.  117.  Meyers, J.P. and R.M. Darling, Model of Carbon Corrosion in PEM Fuel Cells. Journal of the Electrochemical Society, 2006. 153(8): p. A1432A1442-A1432-A1442.  118.  Zhang, L., et al., Gas flow rate distributions in parallel minichannels for polymer electrolyte membrane fuel cells: Experiments and theoretical analysis. Journal of Power Sources, 2010. 195(10): p. 3231-3239.  119.  Owejan, J.P., et al., Water management studies in PEM fuel cells, Part I: Fuel cell design and in situ water distributions. International Journal of Hydrogen Energy, 2009. 34(8): p. 3436-3444.  120.  Shyam Prasad, K.B., S. Maharudrayya, and S. Jayanti, Flow maldistribution in interdigitated channels used in PEM fuel cells. Journal of Power Sources, 2006. 159(1): p. 595-604.  121.  Zhang, W., et al., Analysis and optimization of flow distribution in parallelchannel configurations for proton exchange membrane fuel cells. Journal of Power Sources, 2009. 194(2): p. 931-940.  122.  Basu, S., C.-Y. Wang, and K.S. Chen, Two-Phase Flow Maldistribution and Mitigation in Polymer Electrolyte Fuel Cells. Journal of Fuel Cell Science and Technology, 2009. 6(3): p. 031007-11.  123.  Liu, H., P. Li, and J.V. Lew, CFD study on flow distribution uniformity in fuel distributors having multiple structural bifurcations of flow channels. International Journal of Hydrogen Energy, 2010. 35(17): p. 9186-9198.  124.  Wang, C.T., et al., Fuel Cell Bionic Flow Slab Design. Journal of Fuel Cell Science and Technology, 2010. 7(1): p. 011009-5.  125.  Hu, P., et al., Optimization design of slotted-interdigitated channel for stamped thin metal bipolar plate in proton exchange membrane fuel cell. Journal of Power Sources, 2009. 187(2): p. 407-414.  126.  Zhang, L., et al., Gas-liquid two-phase flow distributions in parallel channels for fuel cells. Journal of Power Sources, 2009. 189(2): p. 10231031.  219  127.  Berg, P., et al., Water Management in PEM Fuel Cells. Journal of the Electrochemical Society, 2004. 151(3): p. A341-A353-A341-A353.  128.  Brackbill, J.U., D.B. Kothe, and C. Zemach, A continuum method for modeling surface tension. Journal of Computational Physics, 1992. 100(2): p. 335-354.  129.  FLUENT 6.3 User's Guide, Fluent Inc., 2006.  130.  K.W. Morton, Numerical Methods for Fluid Dynamics. Academic Press; 1983.  131.  Theodorakakos, A., et al., Dynamics of water droplets detached from porous surfaces of relevance to PEM fuel cells. Journal of Colloid and Interface Science, 2006. 300(2): p. 673-687.  132.  Markicevic, B., 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, 2007. 171(2): p. 706717.  133.  Gostick, J.T., et al., Pore network modeling of fibrous gas diffusion layers for polymer electrolyte membrane fuel cells. Journal of Power Sources, 2007. 173(1): p. 277-290.  134.  Sinha, P.K. and C.-Y. Wang, Pore-network modeling of liquid water transport in gas diffusion layer of a polymer electrolyte fuel cell. Electrochimica Acta, 2007. 52(28): p. 7936-7945.  135.  Lee, K.J., J.H. Nam, and C.J. Kim, Pore-network analysis of two-phase water transport in gas diffusion layers of polymer electrolyte membrane fuel cells. Electrochimica Acta, 2009. 54(4): p. 1166-1176.  136.  Schladitz, K., et al., Design of acoustic trim based on geometric modeling and flow simulation for non-woven. Computational Materials Science, 2006. 38(1): p. 56-66.  137.  Müller-Steinhagen, H. and K. Heck, A simple friction pressure drop correlation for two-phase flow in pipes. Chemical Engineering and Processing: Process Intensification, 1986. 20(6): p. 297-308.  138.  Mishima, K. and T. Hibiki, Some characteristics of air-water two-phase flow in small diameter vertical tubes. International Journal of Multiphase Flow, 1996. 22(4): p. 703-712.  139.  Tran, T.N., et al., Two-phase pressure drop of refrigerants during flow boiling in small channels: an experimental investigation and correlation development. International Journal of Multiphase Flow, 2000. 26(11): p. 1739-1754.  220  140.  Garimella, S., A. Agarwal, and J.D. Killion, Condensation Pressure Drop in Circular Microchannels. Heat Transfer Engineering, 2005. 26(3): p. 28-35.  141.  Sun, L. and K. Mishima, Evaluation analysis of prediction methods for twophase flow pressure drop in mini-channels. International Journal of Multiphase Flow, 2009. 35(1): p. 47-54.  142.  Cindrella, L., et al., Gas diffusion layer for proton exchange membrane fuel cells--A review. Journal of Power Sources, 2009. 194(1): p. 146-160.  143.  Bevers, D., R. Rogers, and M. von Bradke, Examination of the influence of PTFE coating on the properties of carbon paper in polymer electrolyte fuel cells. Journal of Power Sources, 1996. 63(2): p. 193-201.  144.  Lim, C. and C.Y. Wang, Effects of hydrophobic polymer content in GDL on power performance of a PEM fuel cell. Electrochimica Acta, 2004. 49(24): p. 4149-4156.  145.  Basu, S., C.-Y. Wang, and K.S. Chen, Analytical model of flow maldistribution in polymer electrolyte fuel cell channels. Chemical Engineering Science, 2010. 65(23): p. 6145-6154.  146.  Maharudrayya, S., S. Jayanti, and A.P. Deshpande, Flow distribution and pressure drop in parallel-channel configurations of planar fuel cells. Journal of Power Sources, 2005. 144(1): p. 94-106.  147.  Wang, J., Pressure drop and flow distribution in parallel-channel configurations of fuel cells: Z-type arrangement. International Journal of Hydrogen Energy, 2010. 35(11): p. 5498-5509.  148.  Mohan, G., et al., Analysis of Flow Maldistribution of Fuel and Oxidant in a PEMFC. Journal of Energy Resources Technology, 2004. 126(4): p. 262270.  149.  Bi, H.T., P. Sauriol, and J. Stumper, Two-phase flow distributors for fuel cell flow channels. Particuology, 2010. 8(6): p. 582-587.  150.  Zhang, L., et al., Gas-liquid two-phase flow patterns in parallel channels for fuel cells. Journal of Power Sources, 2008. 183(2): p. 643-650.  151.  Johnson, M.C., et al., United States Patent: 6586128 - Differential pressure fluid flow fields for fuel cells, I. Ballard Power Systems, Editor 2003.  152.  Wilkinson, D.P., O.R. Vanderleeden, and J. Zimmerman, United States Patent: 6541145 - Flow fields for supporting fluid diffusion layers in fuel cells, I. Ballard Power Systems, Editor 2003.  221  153.  Zhang, H., et al., The conception of in-plate adverse-flow flow field for a proton exchange membrane fuel cell. International Journal of Hydrogen Energy, 2010. 35(17): p. 9124-9133.  154.  Ding, Y., H.T. Bi, and D.P. Wilkinson, Three-dimensional numerical simulation of water droplet emerging from a gas diffusion layer surface in micro-channels. Journal of Power Sources, 2010. 195(21): p. 7278-7288.  155.  Sund n, B., Transport phenomena in fuel cells2005, Southampton,Boston: WIT Press.  156.  Feser, J.P., A.K. Prasad, and S.G. Advani, Experimental characterization of in-plane permeability of gas diffusion layers. Journal of Power Sources, 2006. 162(2): p. 1226-1231.  157.  Sun, W., B.A. Peppley, and K. Karan, Modeling the Influence of GDL and flow-field plate parameters on the reaction distribution in the PEMFC cathode catalyst layer. Journal of Power Sources, 2005. 144(1): p. 42-53.  158.  Peng, L., et al., Design, Optimization, and Fabrication of SlottedInterdigitated Thin Metallic Bipolar Plates for PEM Fuel Cells. Journal of Fuel Cell Science and Technology, 2011. 8(1): p. 011002-8.  159.  FLUENT 6.3 User's Guide,Fluent Inc.2006.  160.  Quan, P. and M.C. Lai, Numerical study of water management in the air flow channel of a PEM fuel cell cathode. Journal of Power Sources, 2007. 164(1): p. 222-237.  161.  http://www.accudynetest.com/polytable_03.html?sortby=contact_angle [Accessed: May, 2012].  222  

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
China 58 10
Romania 12 0
United States 11 1
India 5 0
Canada 4 3
Algeria 4 5
Russia 4 0
Japan 3 1
Iran 2 0
Ukraine 2 0
United Arab Emirates 1 0
Portugal 1 0
Germany 1 5
City Views Downloads
Beijing 47 2
Unknown 26 6
Shanghai 5 0
Annaba 4 5
Shenzhen 4 2
Ashburn 3 0
Vancouver 3 2
Odesa 2 0
Washington 2 0
Wuhan 1 0
Moscow 1 0
Shahrak-e Pars 1 0
Severn 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

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-0059060/manifest

Comment

Related Items