Two-Phase Flow in Microchannels with Application to PEM Fuel Cells

by

Te-Chun Wu
B.Sc., Chung Yuan Christian University, 1993
M.Sc., National Cheng Kung University, 1995
Ph.D., National Cheng Kung University, 2000

A Dissertation Submitted in Partial Fulfillment of the

Requirements for the Degree of

DOCTOR OF PHILOSOPHY

in the Department of Mechanical Engineering

 Te-Chun Wu, 2015

University of Victoria

All rights reserved. This thesis may not be reproduced in whole or in part, by photocopy
or other means, without the permission of the author.
 ii
Two-Phase Flow in Microchannels with Application to PEM Fuel Cells

by

Te-Chun Wu
B.Sc., Chung Yuan Christian University, 1993
M.Sc., National Cheng Kung University, 1995
Ph.D., National Cheng Kung University, 2000

Supervisory Committee

Dr. Ned Djilali, Supervisor

(Department of Mechanical Engineering, University of Victoria)

Dr. Rustom Bhiladvala, Departmental Member

(Department of Mechanical Engineering, University of Victoria)

Dr. Alexandra Branzan Albu, Outside Member

(Department of Electrical and Computer Engineering, University of Victoria)
 iii

Supervisory Committee

Dr. Ned Djilali, Supervisor

(Department of Mechanical Engineering, University of Victoria)

Dr. Rustom Bhiladvala, Departmental Member

(Department of Mechanical Engineering, University of Victoria)

Dr. Alexandra Branzan Albu, Outside Member

(Department of Electrical and Computer Engineering, University of Victoria)

ABSTRACT
The performance of PEM fuel cells (PEMFC) relies on the proper control and
management of the liquid water that forms as a result of the electrochemical process,
especially at high current densities. The liquid water transport and removal process in the
gas flow channel is highly dynamic and many of its fundamental features are not well
understood. This thesis presents an experimental and theoretical investigation of the
emergence of water droplets from a single pore into a microchannel. The experiments are
performed in a 250 µm × 250 µm air channel geometry with a single 50 µm pore that
replicates a PEMFC cathode gas channel. A droplet manipulation platform is constructed
using a microfluidics soft lithographic process to allow observation of the dynamic nature
of the water droplets. Flow conditions that correspond to typical operating conditions in a
PEMFC are selected. A test matrix of experiments comprised of different water injection
velocities and air velocities in the gas microchannel is studied. Emergence, detachment
and subsequent dynamic evolution of water droplets are analyzed, both qualitatively and
quantitatively. Quantitative image analysis tools are implemented and applied to the time-
resolved images to document the time evolution of the shape and location of the droplets,
 iv
characteristic frequencies, dynamic contact angles, flow regime and stability maps.
Three different flow regimes are identified, slug, droplet, and film flow. The effects of
the air flow rate and droplet size on the critical detachment conditions are also
investigated.
Numerical simulations using Volume-of-Fluid method are presented to investigate the
water dynamics in the droplet flow. The focus of the modeling is on methods that account
for the dynamic nature of the contact line evolution. Results of different approaches of
dynamic contact angle formulations derived empirically and by using the theoretically
based Hoffmann function are compared with the static contact angle models used to date.
The importance of the dynamic formulation as well as the necessity for high numerical
resolution is highlighted. The Hoffmann function implementation is found to better
capture the salient droplet motion dynamics in terms of advancing and receding contact
angle and periodicity of the emergence process.
To explore the possibility of using the pressure drop signal as a diagnostic tool in
operational fuel cells that are not optically accessible, a flow diagnostic tool was
developed based on pressure drop measurements in a custom designed two-phase flow
fixture with commercial flow channel designs. Water accumulation at the channel outlet
was found to be the primary cause of a low-frequency periodic oscillation of pressure
drop signal. It is shown that the flow regimes can be characterized using the power
spectrum density of the normalized pressure drop signal. This is used to construct a flow
map correlating pressure drop signals to the flow regimes, and opens the possibility for
practical flow diagnostics in operating fuel cells.
 v

Table of Contents

Supervisory Committee ...................................................................................................... ii

ABSTRACT ....................................................................................................................... iii
Table of Contents ................................................................................................................ v
List of Tables ................................................................................................................... viii
List of Figures .................................................................................................................... ix
Acknowledgments............................................................................................................ xiii
Dedication ........................................................................................................................ xiv
1 Introduction ................................................................................................................. 1
1.1. Background and Motivation ............................................................................... 1
1.2. Scope and Organization of the Thesis ................................................................ 3
2 Literature Review........................................................................................................ 5
2.1. Dynamic Behavior of Liquid Water ................................................................... 5
2.2. Flow Regime ....................................................................................................... 7
2.3. Pressure Drop ...................................................................................................... 9
2.4. Effect of Channel Geometry and Surface Properties ........................................ 12
2.5. Status of Prediction Method.............................................................................. 15
2.5.1. Force balance method on droplet detachment .......................................... 15
2.5.2. Volume-of-Fluid (VOF) method............................................................... 18
2.5.3. Level set method (LSM) ........................................................................... 20
2.5.4. Lattice Boltzmann method (LBM)............................................................ 21
2.5.5. Modeling of the dynamic contact angle .................................................... 24
3 Experimental Investigation of Water Droplet Emergence ........................................ 29
3.1. Method and Apparatus ...................................................................................... 29
3.1.1. Microchannel design ................................................................................. 29
3.1.2. Measurement apparatus ............................................................................ 31
3.1.3. Flow conditions ......................................................................................... 33
 vi
3.2. Results and Discussions .................................................................................... 34
3.2.1. Flow regimes ............................................................................................. 35
3.2.2. Droplet emergence frequency ................................................................... 38
3.2.3. Further image analysis of droplet emergence ........................................... 41
3.2.4. Dynamic contact angle .............................................................................. 42
3.2.5. Droplet dynamic at the onset of detachment............................................. 43
3.3. Summary ........................................................................................................... 47
4 Numerical Simulation Using VOF Method .............................................................. 49
4.1. Volume of Fluid Method .................................................................................. 50
4.2. Implementation of Dynamic Contact Angle ..................................................... 54
4.2.1. CFD-ACE+ ............................................................................................... 55
4.2.2. ANSYS FLUENT ..................................................................................... 60
4.3. Simulation Domain and Mesh, Boundary and Initial Conditions ..................... 61
4.3.1. CFD-ACE+ ............................................................................................... 61
4.3.2. ANSYS FLUENT ..................................................................................... 62
4.4. Results and Discussions .................................................................................... 67
4.4.1. CFD-ACE+ ............................................................................................... 67
4.4.2. ANSYS FLUENT ..................................................................................... 73
4.5. Summary ........................................................................................................... 83
5 Pressure Signature and Diagnostic Tool ................................................................... 85
5.1. Introduction ....................................................................................................... 85
5.2. Method and Apparatus ...................................................................................... 86
5.3. Results and Discussions .................................................................................... 91
5.3.1. Characterization of mean pressure drop ................................................... 91
5.3.2. Dynamic characteristics ............................................................................ 96
5.4. Summary ......................................................................................................... 100
6 Conclusion and future works .................................................................................. 101
6.1. Conclusion ...................................................................................................... 101
6.2. Future works ................................................................................................... 102
7 Bibliography ........................................................................................................... 104
8 Appendix A – Microfluidic chip fabrication .......................................................... 112
 vii
9 Appendix B – MATLAB source code .................................................................... 113
10 Appendix C – FLUENT UDF source code for DCA .............................................. 116
11 Appendix D – FLUENT UDF source code for droplet impact ............................... 120
12 Appendix E – Pressure drop measurement data ..................................................... 121
 viii

List of Tables
Table 2.1. Slip models examined by Shikhmurzaev [76] ................................................. 25
Table 3.1. Flow inlet conditions. ...................................................................................... 34
Table 4.1. Properties of liquid and flow conditions. ......................................................... 64
Table 5.1. Test conditions matrix ..................................................................................... 90
Table 5.2. Pressure drop measurement of Plastic-0A using single water injection. ......... 93
Table 12.1. Pressure drop measurement of Plastic-0A using dual water injection......... 121
Table 12.2. Pressure drop measurement of Plastic-02 using single water injection. ...... 122
Table 12.3 Pressure drop measurement of CFP-010A using single water injection. ..... 122
Table 12.4 Pressure drop measurement of CFP-0D using single water injection. ......... 122
 ix

List of Figures
Figure 1.1. An emerged water droplet from the GDL entering the cathode gas flow
channel (reproduced from [89] with permission of Journal of Power Sources). ................ 2
Figure 2.1. Typical flow patterns in gas channel of PEMFC (reproduced from [15] with
permission of Journal of Power Sources). .......................................................................... 8
Figure 2.2. Correlation between fluctuations in cathode P signal and cell voltage
(reproduced from [3] with permission of Heat Transfer Engineering). ............................ 10
Figure 2.3. Two-phase friction multiplier versus the superficial air velocity under
different superficial water velocities (reproduced from [24] with permission of Chemical
Engineering Progress). ...................................................................................................... 11
Figure 2.4. Cross sectional view of typical flow channel. r, rib width; c, channel width; d,
channel depth; α, wall angle (reproduced from [26] with permission of International
Journal of Hydrogen Energy)............................................................................................ 12
Figure 2.5. Schematic of channel design (left) and integration into PEMFC (right)
(reproduced from [30] with permission of Sensors and Actuators A). ............................. 14
Figure 2.6. Schematic of (a) droplet height and contact angle, (b) droplet subjected to a
shear flow with resulting deformation and dynamic contact angles (c) spherical droplet
geometry (d) control volume (reproduced from [45] with permission of ASME
Conference Proceedings). ................................................................................................. 16
Figure 3.1. (a) PDMS chip for droplet manipulation. (b) Cross sectional view of chip. (c)
Field of view in microscope. ............................................................................................. 31
Figure 3.2. Schematic diagram of experimental apparatus. .............................................. 32
Figure 3.3. Typical flow regime in microchannel............................................................. 35
Figure 3.4. Flow map of water emergence phenomena in a model PEMFC cathode gas
microchannel. .................................................................................................................... 37
Figure 3.5. (a) Time domain signal and (b) frequency distribution of droplet emergence
process............................................................................................................................... 39
 x
Figure 3.6. Emergence frequency in droplet flow regime under different flow
conditions. ......................................................................................................................... 40
Figure 3.7. Time resolved images of water emerging from a 50 µm square pore in a 250
µm square gas microchannel with the flow condition of Case 2. a) t = 1 ms, b) t = 3 ms,
c) t = 5 ms, d) t = 10 ms, e) t = 15 ms, f) t = 20 ms, g) t = 25 ms, h) t = 45 ms, i) t = 65
ms, j) t = 75 ms ................................................................................................................. 41
Figure 3.8. Dynamic contact angle evolution through an emergence cycle (13.2 Hz under
Case 2 flow conditions). Period I: surface tension force dominant. Period II: transition.
Period III: drag force dominant......................................................................................... 43
Figure 3.9. Effect of air flow velocity on characteristic droplet size (chord C and height
H) at detachment. .............................................................................................................. 44
Figure 3.10. Effect of air flow on contact angle hysteresis. ............................................. 45
Figure 3.11. Contact angle interpretations and effect of airflow on droplet aspect ratio at
the onset of detachment. ................................................................................................... 46
Figure 4.1. Image of water droplet subjected to the air flow stream. Points A and B are
the receding and advancing points in a 2-D plane of view, whereas r and a designate the
receding and advancing contact angle, respectively. ........................................................ 56
Figure 4.2. Position of advancing point (XB), velocity function and DCA distribution for
droplet emergence cycle Case 2. ....................................................................................... 57
Figure 4.3. Velocity dependent contact angle function from droplet emergence
experiment......................................................................................................................... 58
Figure 4.4. Schematic diagram of capillary rise experiment. ........................................... 59
Figure 4.5. Velocity dependent contact angle function. ................................................... 59
Figure 4.6. Three-dimensional domain and mesh for the numerical simulations of droplet
emergence using CFD-ACE+. .......................................................................................... 61
Figure 4.7. Illustration of numerical grids used for the droplet impact computations. The
region of an adaptive refinement is presented. ................................................................. 63
Figure 4.8. Computation domain and mesh illustration for the droplet emergence study,
presenting the region of an adaptive refinement. .............................................................. 65
Figure 4.9. (a) Top view. (b) Side view of the instant at droplet detachment using SCA. 67
 xi
Figure 4.10. Comparison of contact angle evolution of VOF simulation (SCA, mesh
12.5um) and experiments. ................................................................................................. 68
Figure 4.11. (a) Side view of the instant at droplet detachment using DCA Eq. (4.17). (b)
Comparison of contact angle evolution of VOF simulation (DCA, method 1, mesh 12.5
m) and experiment. ......................................................................................................... 69
Figure 4.12. (a) Side view of the instant at droplet detachment using DCA Eq. (4.19) and
[3]. (b) Comparison of contact angle evolution of VOF simulation (DCA, method 2, mesh
12.5 m) and experiment. ................................................................................................. 71
Figure 4.13. (a) Side view of the instant at droplet detachment using DCA Eq. (4.19) and
[3]. (b) Comparison of contact angle evolution of VOF simulation (DCA, method 2, mesh
6.25 m) and experiment. ................................................................................................. 72
Figure 4.14. Comparison of time sequence of water droplet impact onto wax surface (We
= 90), experiment (left) and numerical (right) (reproduced from [96] with permission of
Experimental Thermal and Fluid Science). ...................................................................... 73
Figure 4.15. Time series images during the spreading phase for SCA modeling of droplet
impact on wax surface. ..................................................................................................... 74
Figure 4.16. Time series images during the recoiling phase for SCA modeling of droplet
impact on wax surface. ..................................................................................................... 75
Figure 4.17. Time series images during the spreading phase for DCA modeling of droplet
impact on wax surface. ..................................................................................................... 76
Figure 4.18. Time series images during the recoiling phase for DCA modeling of droplet
impact on wax surface. ..................................................................................................... 77
Figure 4.19. Schematics of spreading diameter and apex height of drop impacts............ 78
Figure 4.20. Numerical simulation of the temporal evolution of the spread diameter in
comparison with the results of Sikalo et. al [81]. ............................................................. 78
Figure 4.21. Numerical simulation of the temporal evolution of the apex height in
comparison with the results of Šikalo et. al [99]. ............................................................. 79
Figure 4.22. One cycle of droplet emergence time resolved images of numerical results
using static contact angle (SCA, s = 110) approach in FLUENT. Air inlet velocity, Va =
10m/s; water inlet velocity, Vw = 0.04 m/s. ...................................................................... 80
 xii
Figure 4.23. One cycle of droplet emergence time resolved images of numerical results
using dynamic contact angle (DCA) approach in FLUENT. Air inlet velocity, Va = 10
m/s; water inlet velocity, Vw = 0.04 m/s. .......................................................................... 81
Figure 4.24. Comparison of experimental result and the evolution of dynamic contact
angle using FLUENT. ....................................................................................................... 82
Figure 5.1. Experimental setup in pressure drop measurement. ....................................... 87
Figure 5.2. Illustration of the pressure measurement and water injection ports. .............. 87
Figure 5.3. Flow channel plates (Ballard Power System). Unit: mm ............................... 88
Figure 5.4. Schematics and dimensions of the flow channels. ......................................... 88
Figure 5.5. Real time signal of flow rate versus pressure drop in flow plate Plastic-0A. 89
Figure 5.6. Pressure drop calibration curve for various flow plates. ................................ 90
Figure 5.7. Pressure drop signal of Plastic-0A using single injection under the Case 11
test condition. .................................................................................................................... 92
Figure 5.8 The effects of air flow (ReA) and flow channel geometry on the difference
between wet and dry pressure drop under single water injection. (a) Plastic-0A and
Plastic-02, (b) CFP-010A and CFP-0D. ........................................................................... 94
Figure 5.9. The effects of air flow (ReA) and flow channel geometry on the friction
multiplier under single water injection. (a) Plastic-0A and Plastic-02, (b) CFP-010A and
CFP-0D. ............................................................................................................................ 95
Figure 5.10. Real-time pressure drop signal of Plastic-0A channel using single water
injection under Case 34 flow conditions. .......................................................................... 96
Figure 5.11. Snap shot images of water accumulation and drainage process. .................. 97
Figure 5.12. Typical real time pressure drop signal of CFP-010A channel using single
water injection under Case 14 flow condition. ................................................................. 97
Figure 5.13. Dominant frequency in two-phase flow regime of CFP-010A. ................... 98
Figure 5.14. Power spectrum of Case 11, 14, 41 and 44 of CFP 010A. ........................... 99
Figure 5.15. Flow regime identification using DFT (Discrete Fourier Transform) power
of the pressure drop signal. ............................................................................................... 99
 xiii

Acknowledgments

I would express my thankfulness to my supervisor, Dr. Ned Djilali, who supported me

through these years and offered invaluable guidance in many aspects not only as an

advisor but also as a friend.

My special thanks go to Dr. Jay Sui who spent countless time with me in the Energy

Systems and Transport Phenomena (ESTP) lab and his office for the supporting

discussions and friendly encouragement during my research work.

I would like to thank the group members in ESTP and Institute for Integrated Energy

Systems (IESVic). I am fortunate to work with such many talented students, researchers

and staff.

I would also like to thank my family in Taiwan, including my parents and parents-in-

law, for their everlasting encouragement and inspiration through my adventures.

Last but not least, this milestone wouldn't complete without the deepest love from my

wife - Anita Fang, my children - Vicky and Brian. Your overflowing support and

understanding builds up every tiny step in our life to fulfill a huge dream for the future.
 xiv

Dedication

To Anita, Vicky and Brian

It’s simply not possible to continue without your inspiration.
Chapter 1

1 Introduction

1.1. Background and Motivation

PEM fuel cells (PEMFCs) have been the focus of intense research in the last decade
because of their potential to produce clean electricity efficiently for vehicles and
stationary applications as well as their central role in the hydrogen economy. Successful
deployment of PEMFCs into automobile and residential applications will depend on
successfully addressing a number of issues, including water management, more
specifically, the liquid water transport and removal process. Water is generated from the
electrochemical reaction between hydrogen and oxygen and is the major by-product of
the reaction in PEMFCs. It is usually found at the cathode side, although it is also present
in the anode. It is important to remove liquid water from the PEMFC to prevent flooding
of the electrodes and the blocking of pores in the gas diffusion layer (GDL). Conversely,
a sufficient amount of water/hydration is necessary to maintain high protonic
conductivity in the membrane.
However, the presence of some liquid water is unavoidable in most PEMFCs,
particularly at higher current densities, where gas-liquid two-phase flow occurs in the
porous electrodes and flow channels. The occurrence of two-phase flow leads to different
flow regimes and flow patterns, and is accompanied by a substantial pressure drop. The
two-phase flow regimes increase the complexity of water management. Liquid water
must be transported out of the catalyst layer, through the GDL, and into the flow
 2
channels, before being removed by the gas flow. Although practical experience and
empirically based design tools have allowed progress in water management techniques,
the processes are not well understood. Not much detailed information has been obtained
on the flow dynamics and associated flow regimes due to the difficulty of measuring and
directly observing water removal processes in a PEMFC.
The mechanisms of liquid water transport in the flow channel of a PEMFC include the
initial emergence of droplets in the GDL, their subsequent growth and detachment from
the GDL, and finally the interaction with the walls. These mechanisms lead to complex
two-phase flow throughout the channel. Figure 1.1 presents a schematic of a water
droplet emerging from the GDL and entering the cathode gas flow channel.

Figure 1.1. An emerged water droplet from the GDL entering the cathode gas flow
channel (reproduced from [89] with permission of Journal of Power Sources).

In classical two-phase flows in channels [1] it is often the case that both the liquid and
gas are introduced together into a flow channel or the gas phase is produced on the
channel wall due to boiling in a flow dominated by the liquid phase. Two-phase flows in
PEMFC flow channels are fundamentally different in several respects:
1. The gas phase is dominant (i.e. low saturation or high void fraction).
2. Water emerges from a porous layer (GDL) that interacts with the gas flow in the
channel.
 3
3. The channel walls have mixed surface wetting properties; three of the walls
comprised of the bipolar plate are typically hydrophilic, while the fourth wall
corresponding to the GDL is porous, rough and hydrophobic.
4. The length, time and flow scales and corresponding characteristic non-
dimensional numbers (Reynolds number, Capillary and Webber number) differ
significantly from those of “classical” two-phase flows.
These characteristics and the larger role of surface forces make the complexity of two-
phase flow in PEMFC flow channels challenging to analyze and predict. A further
complication is introduced by the variety of mechanisms by which water is transported
and produced in the fuel cell, including electro-osmotic drag, back diffusion, water
production from electrochemical reactions, etc. [2]. Finally, significant temperature
gradients can be present in a fuel cell and strongly impact the relative humidity and
saturation pressure, which affect the rate of phase change.
The size of the flow channel plays an important role in two-phase flow. It can
radically affect the flow regime and, hence, the fuel cell performance. A typical range for
hydraulic diameters of flow channels in PEMFCs is from 200 µm to 3 mm [3]. The flow
is laminar but unsteady; the liquid phase is dominated by surface forces instead of
volume forces and is strongly affected by the hydrophilic or hydrophobic properties of
the surface. The design of the cathode flow channel is the most important consideration
for water management [3]. Hence, improved understanding and prediction of two-phase
flow in the cathode flow microchannel of PEMFCs could pave the way for new or better
concepts in channel design as well as water management.

1.2. Scope and Organization of the Thesis

The objective of the work described in this thesis is to quantify the liquid water
droplet dynamics as well as the two-phase flow phenomena relevant to PEM fuels
cells. The approach combines experimental observations with numerical modelling.
Chapter 2 presents a literature review synthesizing experimental and theoretical work,
and open research question. In Chapter 3 we present flow visualization experiments in
laboratory flow channels relevant to fuel cells to observe and quantify the dynamics of
 4
liquid water droplets. Chapter 4 presents computational modelling of the two-phase flow
with a focus on the physical representation of the dynamic process. The correlation
between pressure signals and the two-phase flow regimes is investigated experimentally
in Chapter 5 as a preliminary step to diagnostics for operating fuel cells. Chapter 6
summarizes the contributions and outlines recommendations for future research.
 5

Chapter 2

2 Literature Review
In this section, reviews of salient features of two-phase flow in PEMFC microchannels
are addressed. These include the dynamic behavior of liquid water, the different flow
regimes, the pressure drop along the channel, the effects of channel geometry and surface
properties of the boundaries.

2.1. Dynamic Behavior of Liquid Water

Water management affects performance and durability of PEM fuel cells (PEMFCs),
and remains a pacing item in the development of commercial stacks. In spite of a number
of experimental and theoretical studies, the underlying liquid water transport, including
water droplet generation, growth, detachment and removal processes, are not well
understood. This stems from the modelling challenges associated with the complex
dynamics of the two-phase flow, and from the optically inaccessible and
electrochemically active nature of a fuel cell that make in-situ measurements difficult.
The experimental challenges and a review of the techniques used to date to visualize
liquid water transport are discussed by Bazylak [4]. One of the most widely used
approaches is optical visualization/photography in custom-designed transparent PEMFCs.
This approach is exemplified by the work of Tüber et al. [5] who, based on their
observations, concluded that hydrophilic diffusion layers resulted in increased current
densities and better fuel cell performance as a result of more uniform membrane
hydration. Yang et al. [6] visualized water droplet emergence from the GDL surface and
its subsequent behavior in the gas channel. They found that droplets appear only at
 6
certain preferential locations, and can grow to a size comparable to the channel
dimensions. A variety of phenomena in the gas channel were reported, including the
intermittent emergence of droplets, film formation and channel clogging. Chen et al. [7]
proposed a simplified force balance model based on droplet geometry approximations for
predicting the onset of water droplet instability on GDL surfaces, and supported the
model through their experimental observations. They found that droplet removal can be
enhanced by increasing the flow channel length or mean gas flow velocity, by decreasing
the channel height or contact angle hysteresis (the difference between the advancing and
receding contact angles), or by making the GDL more hydrophobic. Two different modes
of liquid water removal were identified by Zhang et al. [8], and based on a further
analysis using the force balance model, they proposed a relationship between the droplet
detachment diameter and air velocity. Independently of Chen et al. [7], Kumbur et al. [9]
conducted a similar theoretical and experimental study of the influence of controllable
engineering parameters, including surface PTFE coverage, channel geometry, droplet
chord dimensions, and operational air flow rate. The proposed water droplet instability
criterion was formulated in terms of the Reynolds number and droplet aspect ratio.
Theodorakakos et al. [10] constructed a platform consisting of a high aspect ratio channel
(2.7 mm × 7.0 mm) for visualizing the behavior of a single droplet placed on the surface
of different GDLs. They obtained side-view droplet detachment images and measured the
dynamic contact angles for input into their in-house volume-of-fluid (VOF) computer
simulations, and were able to correlate the critical droplet diameter at detachment with
the air velocity. The channel dimensions in these studies [9,10] were significantly larger
than PEMFC gas flow channels, which have hydraulic diameters ranging from 200 µm
to 3 mm [3], and though valuable, the insights from these studies are expected to have
limited relevance to PEMFCs because of the important impact of reduced dimensions on
the two-phase flow. Whereas previous experimental studies [5–10] investigated water
droplets or films initially resting on a GDL, the actual water transport process in a
PEMFC channel also involves the initial emergence of droplets into the GDL, their
subsequent growth and detachment from the GDL and, finally, their interaction with the
walls leading to a complex two-phase flow regime throughout the channel. The resulting
flow regimes may vary between surface tension dominated slug flow and inertia
 7
dominated annular flow (i.e. film flow as shown in Figure 2.1) depending on the load
[3]. The annular or dispersed droplet flow regime is desirable in fuel cell operation since
this provides a path for reactant flow in the presence of water. However, due to the
operational load of fuel cell stacks, the flow regime is at times likely to shift to slug flow
because at low power the shear force in the cathode channels is usually insufficient to
overcome the surface tension forces that hold water within the flow field channels and
GDLs.

2.2. Flow Regime

While microchannel two-phase flows have been the focus of many experimental and
theoretical studies targeted at PEMFCs (see review by Anderson et al. [11]), such flows
are also of relevance in applications ranging from heat sinks for electronic devices to
microreactors and microfluidic devices. Most of the work in these latter applications is
concerned with acqueous-acqueous systems, but a few studies have recently investigated
gas-liquid systems using flow visualization and particle image velocimetry. Günther et al
[12] investigated different gas–liquid rectangular microfluidic channels and networks
pertinent to chemical reactions and mapped bubbly, slug and annular flow patterns as a
function of the gas and liquid superficial velocities. A flow pattern map based on
dimensional analysis and visualization was also presented by Waelchli et al. [13] for
conditions related to microreactors. Based on a critical review of available observations,
it was noted that basing the similarity analysis on cross-sectional channel shape rather
than hydraulic diameters is key for reliable prediction of flow regime based on
experimental flow maps. This is consistent with the recent study of Kim et al. [14] who
identified three fundamentally different two-phase flow regimes (capillary bubbly,
segmented, annular) together with two transitory ones. They reported that the transition
boundaries depend on the geometry of the test microchannels and of the injection port.
They also noted differences in pressure drop for the capillary bubbly and segmented flow
regimes when using different microchannel materials. Note that these microfluidic
studies were all performed using hydrophilic surfaces as opposed to the hydrophobic
surfaces of the flow channel in PEM fuel cells. The flow regimes identified in these
studies, such as bubbly (or capillary bubbly), slug (or segmented), annular (or film) have
 8
different characteristics and are not all relevant to fuel cell channels. For instance, the
bubbly flow is unlikely occurs in PEM fuel cells since it is only present under very high
liquid velocity conditions (i.e. liquid phase dominated) with minimal gas flow.
As in the case of microfluidic applications, two-phase flow regimes relevant to
PEMFC cathodes can be broadly classified as slug, droplet and film flows, but the
detailed characteristics and transition boundaries differ significantly due to the properties
of the GDL, and particularly hydrophobicity and roughness. Each of the three regimes
can occur solely in a channel, or two or more might occur simultaneously at different
locations of the flow channel [11]. Hussaini and Wang [15] presented an in-situ study of
cathode flooding using a transparent PEMFC. Gas relative humidities of 26%, 42% and
66%, current densities of 0.2, 0.5 and 0.8 A cm−2 and flow stoichiometries ranging from 2
to 4 were used to represent typical operating conditions for automotive applications.
Based on observed flow patterns and the superficial phase velocities, a flow map as show
in Figure 2.1 was developed and utilized to determine operating conditions that prevent
channel flooding.

Figure 2.1. Typical flow patterns in gas channel of PEMFC (reproduced from
[15] with permission of Journal of Power Sources).

However, there are inconsistencies in the observed flow patterns from various
studies. Direct ex-situ visualization of droplet evolution using laboratory models of fuel
cell microchannels has allowed some detailed analysis not possible in-situ. Hidrovo et al.
 9
[16] investigated water slug detachment in two-phase hydrophobic microchannel flows.
Due to the aspect ratio and geometry of the microchannel, water was observed to form
pancake-like slugs rather than a spherical cap droplet. More recently, Lu et al. [17][18]
presented flow visualization and pressure drop measurements over a broad range of flow
regimes, flow parameters, channel surface wettability, geometry and orientation.
Colosqui et al. [19] used an experimental set-up conceptually similar to the present one
but with larger channel cross-sections in which gravity and surface tension forces are of
the same order. Results showed that flow channel geometry and interfacial forces are the
dominant factors in determining the size of slugs and the required pressure drop for their
removal, those residual water droplets can alter the wetting properties and act as
nucleating agents that impact the dynamics of slug formation and detachment. The
interaction between the air and water flows that occur at the gas–liquid interface of a
droplet was examined by Minor et al. [20]. Using micro-digital-particle-image-
velocimetry (micro-DPIV) and examining seeded droplets first placed on a GDL, they
analyzed the relationship between air velocity in the channel, secondary rotational flow
inside a droplet, droplet deformation and contact angle hysteresis.

2.3. Pressure Drop

In operating PEM fuel cells, a low pressure difference across the flow field is generally
desirable to lower auxiliary power demand. On the other hand, a larger pressure drop is
desirable and often necessary to remove water from the flow field channels. Differential
pressure signals have been widely adopted as one of the quickest approaches for
identifying two-phase flow regimes. The pressure drop signal correlates closely with the
two-phase flow regime especially when slug flow occurs. The presence of liquid water in
fuel cell channels hinders the gas flow, thus creating a higher gas pressure gradient
compared with single phase flow. In fuel cell operating conditions, a greater pressure
drop will increase the risk of instability, cost, and lower performance. It is evident that
the flow maldistribution suffers primarily from the increase in the pressure drop due to
uneven water distribution from channel to channel and thus is detrimental to fuel cell
performance and durability. Yousfi-Steiner at al. [21] reviewed the voltage degradation
issues associated with the water flow factors and characterization of water management.
 10
In PEMFCs, significant pressure drop is caused by the frictions of the cell, especially
the gas flow field inside the porous structure of the electrode and the GDL. The other
major contribution to pressure drop in PEMFCs is due to the hydrodynamics of two-
phase flow inside the channel. The pressure drop usually can be measured between inlet
and outlet channels of the gas manifold as well as at the anode or cathode. Quantifying
the contribution due to channel flooding is a key factor and is still a challenge. In the case
of two-phase flow in the channels, Rodatz et al. [22] showed that the existence of liquid
water can reduce the cross sectional area available for gas diffusion which in turn reduces
the gas permeability, leading to an increase of the pressure drop according to Darcy’s
law. Therefore, the pressure drop in the flow fields varies as a function of the flooding
level, and can indicate the presence of liquid water. It is evident that due to the higher
amount of liquid water production and accumulation in the cathode channel, the pressure
drop observed at the cathode side is higher than at the anode side. He et al. [23]
monitored the pressure drop at the cathode and showed a strong correlation between the
flooding level and the amplitude of losses due to mass transfer limitations associated with
cell flooding. The fluctuation in the pressure drop at the cathode along with the cell
voltage of an operating fuel cell was documented by Trabold [3] and is shown in Figure
2.2.

Figure 2.2. Correlation between fluctuations in cathode P signal and cell voltage
(reproduced from [3] with permission of Heat Transfer Engineering).
 11
It clearly demonstrates that the mean cathode pressure drop increases as the cell
voltage decreases. This happens because more water accumulates at low voltages,
resulting in a significant variation in the pressure drop signal. The pressure drop signal
can be further examined to deduce a ratio known as the two-phase friction multiplier
[24],  g , defined as the ratio of two-phase flow pressure drop to the single gas phase
2

pressure drop:
P2
g 
2
(2.1)
Pg

where P2 and Pg are the pressure drops with two-phase flow and with single-phase
flow in the channel, respectively. Wang et al. [25] have demonstrated that this ratio can
be used as a simple indicator for the liquid buildup in a PEMFC channel. Lu et al. [18]
conducted an ex situ investigation of flow maldistribution and the pressure drop effect.
This ratio was presented in the form of superficial air velocities relative to superficial
water velocities as shown in Figure 2.3.

Figure 2.3. Two-phase friction multiplier versus the superficial air velocity under
different superficial water velocities (reproduced from [24] with permission of
Chemical Engineering Progress).
 12
At low air velocities, the two-phase flow friction multiplier is greater than unity due to
liquid water buildup. Thus, the highest  g , indicating high water buildup result in slug
2

flow. When the air velocity is increased,  g will decrease and approach unity. The
2

relationships between pressure drop and flow regimes are highly correlated. As reported
by Lu et al. [18], two-phase flow at low superficial air velocities is dominated by slugs or
semi-slugs and lead to large fluctuations in the pressure drop and severe flow
maldistribution. At higher air velocities, a water film flow regime produced smaller but
more frequent fluctuations in the pressure drop, resulting from the water buildup at the
channel-exit manifold interface. A further increase in the air velocity shifted the flow
regime into mist flow where there is very little water buildup and the pressure drop is
very small.

2.4. Effect of Channel Geometry and Surface

Properties
A typical flow field of a PEM fuel cell consists of a series of channels and ribs with a
cross-sectional area of the channels of the order of a square millimeter. The channel
geometry has a prominent influence on the performance of the fuel cell due to its impact
on flow and water management. The main geometric parameters of the flow field include
length, depth, width and the rib width of the channel. Further geometric characteristics
depend on the production process such as the rib-wall radii, the channel-bottom-wall radii
and the wall angle. A cross sectional view of a typical flow channel is shown in Figure
2.4 [26].

Figure 2.4. Cross sectional view of typical flow channel. r, rib width; c, channel width;
d, channel depth; α, wall angle (reproduced from [26] with permission of International
Journal of Hydrogen Energy).
 13
Scholta et al. [27] investigated fuel cell performance in several different channel
geometries. They found that channel and rib widths in the range of 0.5 to 1.0 mm are
advantageous for fuel cell performance. In addition, it was determined that narrow
channel dimensions are preferred for high current densities, whereas wider dimensions
are better for low current densities. Shimpalee and Van Zee [28] numerically investigated
the influence of rib and channel dimensions for a fixed depth of 0.55 mm and concluded
that for a fuel cell with 200 cm2 active area a wider channel (1.0 mm vs. 0.7 mm) with a
smaller land width (0.7 mm vs. 1.0 mm) improves performance and flow distribution
uniformity. However, for a fuel cell with 100 cm2 active area, the performance effects of
these two parameters depend on the operating conditions.
Aktar et al. [26] studied the effects of channel shape and aspect ratio on pressure drop.
They considered four different rectangular geometries with different aspect ratios and a
triangular geometry. The optimum geometry was found to be a rectangular channel 1 mm
wide and 0.5 mm deep. This particular geometry exhibited the best water removal
capability at a reasonable pressure drop. Similar findings were also reported by Kumber
et al. [9]. The triangular flow channel geometry did not improve the water removal
characteristics or influence the pressure drop enough to move the droplet. Zhu et al. [29]
numerically investigated the effects of different channel geometries on water droplet
dynamics in a single channel. Simulations for microchannels with different cross-
sectional shapes, including rectangles with aspect ratios from 0.1 to 2, a trapezoid, an
upside-down trapezoid, a triangle, a rectangle with a curved bottom wall, and a semicircle
were compared. For the cases of rectangular channel geometries, the longest detachment
time and the largest detachment diameter was seen in the geometry with an aspect ratio
(depth/width) of 0.5. The longest removal time was seen in the geometry with an aspect
ratio of 0.25. However, the pressure drop for the geometry with an aspect ratio of 0.1 was
the highest. They concluded that there is no optimum design in channel geometry in
terms of finding a low pressure drop and efficient water removal.
Considering the effect of capillary forces, Metz et al. [30] presented a secondary
channel design on top of a triangular cathode gas microchannel for passive water removal
as shown in Figure 2.5. It was concluded that cathode walls with low contact angles as
well as opening angles larger than 20 are best suited to facilitate water removal in
 14
realistic situations. Their flow field design stabilized the cell at 95% of its initial
performance compared to 60% when using the standard design without a secondary
channel.

Figure 2.5. Schematic of channel design (left) and integration into PEMFC (right)
(reproduced from [30] with permission of Sensors and Actuators A).

Owejan et al. [31] used neutron imaging to compare the water accumulation in
different GDLs and microchannels. They reported that hydrophobic coating flow field
channels retain more water. However, these channels also contain a greater number of
smaller slugs in the channel area, improving the fuel cell performance at high current
densities. They also found that the triangular channel geometry retained less water than
rectangular channels of the same cross-sectional area, and the water is mostly trapped in
the two corners adjacent to the diffusion media.
More recently, an experimental investigation on the effect of channel geometry and
orientation as well as the wettability was analyzed by Cheah et al. [32]. It was shown that
larger water slugs formed in a hydrophilic channel in spite of reducing surface energy for
water removal but will hinder the mass transfer for the reactant gas to get underneath the
slug and diffuse into catalyst layer. The hydrophobic channel, on the contrary, produces
smaller water slugs and requires more energy for removal but the reactant transport is
better improved. They suggested considering all these factors in order to obtain an
optimal gas channel design.
 15
It is still not clear as to whether the flow channel wall should be more hydrophobic
or more hydrophilic in order to facilitate water removal. Whether or not the water
exhibits a wetting or a non-wetting behavior inside the channels is unlikely due to the
single parameter of hydrophobicity. Rather, the wetting characteristics of water inside the
channel will also depend on the surface material properties, the surface roughness and the
geometry of the microchannel.

2.5. Status of Prediction Method

Water transport is a predominant concern in the design of PEMFC flow
microchannels. The flow regimes differ from classical two-phase flows because of the
confined geometry, non-wetting and rough surfaces, and the dominant effect of surface
tension and surface forces [33,34]. Visualization has been employed extensively to study
these flows [35], but quantitative measurements remain scarce. The potential of Volume-
of-Fluid (VOF) based CFD simulations, Lattice Boltzmann methods (LBM) and Level set
methods (LSM) for fuel cell flows has been demonstrated [36–41], but these remain
limited by the inability to resolve the roughness of the GDL surface and by fundamental
issues regarding the physics of moving wetting lines [42,43].
As described in Section 1.1, the water transport mechanisms in PEMFCs include
electro-osmotic drag, capillary force, back diffusion, and multi-component two-phase
flow. The water droplets passing through the GDL and appearing in the gas flow
channels had been observed and shown to be the source of water. Therefore, accounting
for the corresponding droplet dynamics on the channel surface is a central issue in the
description and modelling of the two-phase flow phenomena involved in gas channel.
The following section will summarize different prediction methods focused on the
droplet dynamic issue.

2.5.1. Force balance method on droplet detachment

Liquid water removal from GDL has been studied thoroughly [7,9,44] and the
corresponding water droplet detachment models were reviewed by Schillberg and
Kandlikar [45]. To better describe the problem, several definitions were introduced to
describe droplet dynamics. The droplet height and static contact angle are schematically
 16
shown in Figure 2.6a. The contact line is the boundary where the three phases meet at
the interface of air, water and the GDL. When a droplet has no transverse forces acting on
it, the contact angle is constant around the entire contact length which is the so-called
static contact angle. When the droplet is experiencing shearing forces in the fuel cell
channel, the droplet generally tilts toward the direction of flow and the contact angle
varies along the line of contact as illustrated in Figure 2.6b.

(a) (b)

(c) (d)
Figure 2.6. Schematic of (a) droplet height and contact angle, (b) droplet subjected to a
shear flow with resulting deformation and dynamic contact angles (c) spherical droplet
geometry (d) control volume (reproduced from [45] with permission of ASME
Conference Proceedings).

The advancing (A) and receding (R) contact angles represent contact angles in the
downstream and upstream directions, respectively. These angles are the so-called
dynamic contact angles and represent the effect of two forces: a drag force and a surface
adhesion force, which balance on the droplet surface and contact line. This variability in
the contact angle around the droplet is a way to measure its deformation and stability.
The drag force tends to move the droplet away from its location and is a sum of the shear
stress force and the pressure force, whereas the surface adhesion force comes from
surface tension which acts to hold the droplet in place, pushing against the flow. The
droplets depart from the surface as long as the drag force exceeds the surface adhesion
 17
force. The subject of displacing liquid droplets from solid surface is a fundamental
problem and has attracted considerable interests by many researchers [46–50]. However,
the large range of conditions experienced by the droplets makes their dynamics quite
complex.
A simplified model that is based on the force balances and droplet-geometry
approximations was presented by Chen et al. [7] to predict the onset of instability leading
to the removal of water droplets at the GDL and gas flow channel interface. A spherical
droplet under a control volume is schematically described in Figure 2.6c-d. The overall
pressure drop is then related to find the drag force and after balancing this force with the
surface tension force, the droplet detachment can be described by Eq. (18) in [7]. They
concluded that the droplet removal can be enhanced by (a) increasing flow channel length
or mean gas flow velocity, (b) decreasing channel height or contact angle hysteresis (the
difference between the advancing and receding contact angles) or (c) making the GDL
more hydrophobic. The effect of mean air flow velocity on the droplet was represented in
terms of an instability window by analyzing the contact angle hysteresis of different
shape of droplet, channel length and flow velocity in Chen et al. [7]. A larger unstable
area in the stability window is desirable for enhanced water removal, as this indicates a
higher probability of droplet detachment.
In comparing the experimental and numerical results, the force balance model
provides a reasonable agreement in describing the droplet dynamics. Based on the data
Chen et al. [7] obtained the following criterion for prevention of clogging of the channel
by water droplets:

 L   gasU  
    (2.2)
 2 B    12
where L is the length of the flow channel, µ gas is the viscosity of flowing gas, U is the
average velocity along the flow direction,  is the surface tension and 2B is the channel
height. The left hand side of Eq. (2.2) represents the product of the channel length-to-
height aspect ratio with the capillary number, and can be treated as an initial estimate of
the flow regime in a gas channel.
 18
Similar to Chen et al. [7], Kumbur et al. [9] conducted a combined theoretical and
experimental study of the influence of controllable engineering parameters on liquid
droplet deformation at the interface of the GDL and the gas flow channel. They proposed
additional parameters including Reynolds number and droplet aspect ratio which affect
the water droplet instability. They used an ex-situ approach consisting of a rectangular 5
mm × 4 mm flow channel. Although the channel dimensions were substantially larger
than those of a typical PEMFC channel, the proposed analysis and results provide a good
representation of salient droplet behavior in the gas channel and indirectly support the
approach leading to Eq. (2.2). The overall trend is that decreasing the channel height
makes droplet detachment more likely; however the predicted onset of
instability/detachment is delayed. Recently an improved analysis that yields improved
agreement with experiments was proposed by Miller [51] based on the same concept as
Chen et al. [7] and Kumbur et al. [9] but with a more rigorous application of the force
balance.
In summary, the most common results in droplet dynamics are [45]:
(1) Increasing hydrophobicity of the GDL tends to decrease the droplet height at
detachment and reduces interaction with channel walls.
(2) Decreasing the channel height or increasing the channel length makes droplet
detachment more likely.
(3) Increasing the gas channel mean velocity will result in a decrease in the droplet
height at departure.
(4) Taller droplets tend to have larger contact angle hysteresis.

2.5.2. Volume-of-Fluid (VOF) method

The VOF method is a computationally intensive method developed in the early 1980s
[52]. With the availability of increasing more powerful computing resources, the method
has gained increasing acceptance to simulate time dependent flows in immiscible multi-
phase systems. In VOF methods, the location of the interfaces is determined by applying
a surface tracking technique to a fixed Eulerian mesh; the method can thus be readily
implemented into established CFD frameworks. A volume fraction indicator is used in
conjunction with a reconstruction technique to determine the location and shape of the
 19
interface. A key feature of the VOF method is the ability to capture the effects of
surface tension using a continuum surface force approach. Surface tension is an important
and often dominant force in microchannel flows.
Quan et al. [53] studied a serpentine single channel using the VOF method. The inlet
air flow velocity was fixed at 10 m/s together with an initial water distribution. This
initial water distribution was by starting the simulations in the presence of a single
droplet or a series of droplets freely suspended at the inlet or attached on the channel
wall. The simulations show that the bend/switchback area plays an important role in
determining water behavior inside a U-shaped microchannel due to the combined effects
of shear stress, wall adhesion, gravitational force, and surface tension. The results
reported in this study though providing some interesting qualitative insights, should be
treated with a good deal of caution because, a) the grid resolution was very coarse and
inadequate, and b) in an actual fuel cell water is transported through the GDL and
emerges into the cathode gas channel as opposed to being freely suspended initially
attached to the wall. Cai et al. [54] studied the mobility of water droplets and films inside
a straight channel and investigated the effects of hydrophilic and hydrophobic properties
of microchannels using the VOF method. Again, initial conditions were prescribed with
water film covering the wall. Results showed that the material properties of the channel
side-wall play an important role in water transport, with faster water removal on a
hydrophobic surface. A hydrophobic membrane electrode assembly (MEA) surface and
hydrophilic channel side walls could prevent water accumulation on the MEA surface.
Zhu et al. [55] introduced an improved, though still idealized, representation of the
process in a PEMFC, by modelling the emergence of water form a pore into the cathode
microchannel. The first sets of exploratory VOF simulations presented by the group were
two dimensional and were later extended to three dimensional [29,36,56]. It was found
that accounting for the initial connection of a droplet to a pore could yield significantly
different dynamics; additionally, the critical air velocity for droplet detachment are also
significantly higher for cases when a droplet is assumed initially stagnant and sitting on
the surface. Simulations were also presented using water emergence from several pores
illustrating the even more complex dynamics arising from droplet to droplet interactions,
and merging. Another recent contribution using VOF is due to Chen [57] who conducted
 20
parametric simulation for different dimensions of a gas and water channel using the
commercial package Flow-3D. He derived an analytical solution using a simplified
explicit model, which is based on a simple force balance between pressure drag and
surface tension. The results showed that by making the GDL surface more hydrophobic,
decreasing contact-angle hysteresis, and decreasing channel height one can reduce the
critical air-flow velocity.
To study the two-phase flow patterns in PEMFCs gas flow channel, Ding et al. [58]
employed a numerical modeling using VOF method. The GDL surface structure was
simplified by opening 320 pores on the surface with the same diameter of 400 µm
considering a typical 50% porosity of open area on GDL surface [59]. Their results show
that the flow pattern evolves from corner droplet flow, film flow, annular flow and finally
slug flow.

2.5.3. Level set method (LSM)

Level set methods treat the location of a propagating interface as the zero value of a
higher order continuous level set function. Given the initial value of the speed of the
propagation of each point on the interface, the interface location can be formulated as an
initial value problem. Level set methods are advantageous in that they do not require any
adaptive meshing of the interface as it propagates, and they handle sharp corners and
merging fronts rather well. In addition, there is no need for reconstruction of the interface
as is done with VOF methods, and level set methods are easily extended to multiple
dimensions [60]. Some of the first work on using level set methods for two phase flow
simulations was done by defining the level set function to be the distance of a given point
from the air/water interface. In air, the distance is defined to be negative, while in water
the distance is defined to be positive. A projection operator is used to reduce the Navier-
Stokes equation, and an equation for the level set function is computed instead of solving
equations for density and viscosity. By solving an equation for the smooth level set
function, level set methods avoid the numerical difficulty of solving equations for density
and viscosity at the air/water interface where there are discontinuities. Instead, the
density, viscosity, and curvature (to compute the surface tension) are computed from the
 21
level set function [61]. To obtain higher resolution of water droplets, adaptive meshing
techniques were later applied to the level set method [62].
Several studies have been done using level set methods to simulate two-phase flow in
PEM fuel cells. Mukherjee and Kandlikar [40] analyzed the growing water droplet inside
a model cathode channel using a level set method applied in two dimensions. The
channel had a cross sectional area of 1 mm2 and a length of 2 mm. They analyzed the
growth and departure of a water droplet inside a channel with air flowing through it. The
droplet departure phenomenon was found to be hindered by excessive water flow rates
and the increase in surface wettability. They also found that the droplet departure
diameter decreased as the contact angle at the channel wall was increased. More recently,
Choi and Son [41] presented the droplet dynamic in a model cathode channel with
multiple poles using a level set method applied in two dimensions. They found that as the
air speed in the microchannel increases, the detachment diameters of water droplets
decrease, the surface area coverage of water decreases, and the time it takes for water to
leave the microchannel decreases. This indicates that the air speed is critical in the
removal of water from the microchannel. Numerical simulations of two pores in the
microchannel indicate that the merger of two water droplets is more likely to occur when
the distance between the pores is decreased or when the pores are oriented in the same
direction as the airflow. Merging water droplets covered more surface area and took
longer to leave the microchannel, indicating that these effects inhibit the removal of water
from the microchannel. Finally, it was observed that droplets formed which touch the
sidewall and bottom wall take longer to leave the microchannel, and are in general larger
than droplets that only touch the bottom wall.

2.5.4. Lattice Boltzmann method (LBM)

LBM is a relatively new technique which has its roots in lattice gas cellular automata
and particle methods. The technique has gained increasing adoption for simulating a wide
range of flows, due in part to the relatively fast and simple numerical algorithms. Unlike
classical CFD methods, which solve the conservation equations (mass, momentum, and
energy) for macroscopic quantities; LBM models the fluid as consisting of particles
which propagate and collide over a discrete lattice mesh. Due to its particulate nature and
 22
local dynamics, LBM has several advantages over other conventional CFD methods,
especially in dealing with complex boundaries (such as porous media), accounting for
microscopic interactions and non-continuum effects (i.e. larger Knudsen number flows),
and parallelization of the algorithm [63]. Intrinsic features enable the LBM to model
phase segregation and interfacial dynamics of multi-phase flow, which are difficult to
handle in Navier-Stokes based CFD methods. More recently, the LBM has been extended
to the simulation of multi-phase fluid flows [64], though issues remain to be addressed to
make the technique fully functional in such flows. A few applications of LBM to study
transport phenomena in PEMFCs have been reported. Fei and Hong [65] simulated the
two-phase flow of methanol/CO2 in a microchannel in a direct methanol fuel cell, and
investigated CO2 bubbly flow phenomena in the microchannel under different operation
conditions. Park and Li [66] presented a 2D two-phase LBM simulation of liquid water
flowing through a model of the fibrous structure of carbon paper (GDL). The LBM study
of water droplet dynamics on a hydrophobic surface of a gas flow channel in PEMFC
was published by Hao and Cheng [67]. They used the LBM multiphase free energy
model to simulate the formation of a water droplet emerging through a micro-pore on the
hydrophobic GDL surface and its subsequent movement under the action of shear flow.
The computational model and domain were similar to those used in Zhu et al [56]; the
domain considered consisted of a 1200 µm long microchannel with a cross section of 600
µm  300 µm a water pore of 90 µm to simulate the water transport from the GDL. The
water and air velocity were fixed at 0.075 m/s and 3.86 m/s, respectively. One benefit
when using mesoscopic level simulations with LBM is the improved description of the
dynamic change of contact angle. As discussed later in Section 2.5.5, this is a key
requirement for physically realistic simulations. Similar results to the VOF results of Zhu
et al. [56] were reported, consistent with well-established trends, i.e. water droplet
removal was found to be facilitated by a high gas flow velocity on a more hydrophobic
GDL surface, and a highly hydrophobic surface induced lifting of water droplet from the
GDL surface. An analytical model based on force balance was also presented to predict
the droplet detachment size. However, characteristic roughness of an actual GDL surface
was not considered in this work.
 23
Several approaches have been proposed to simulate two-phase flow using LBM. For
example, Gunstensen et al. [68] developed a multi-component LBM on the basis of two-
component lattice gas model; Shan and Chen [64] presented a LBM model with mean-
field interactions for multi-phase and multi-component fluid flows; Swift et al. [69]
proposed a LBM model for multi-phase and multi-component flows using the concept of
free energy; He et al. [70] developed a model using an index function to track the
interface of multi-phase flow. Although the LBM methodology is a promising method for
multi-component/phase flows, one key issue is that all the schemes listed above are
limited to small density ratio (typically less than 20), and become numerically unstable
for higher density ratios [71]. This is a serious limitation since in most liquid-gas systems
of interest the density ratio is usually larger than 100, and in the specific case of PEMFCs
the density ratio (water to air) is of order 1000. To overcome this limitation, Inamuro et
al. [72] proposed a LBM for incompressible two-phase flows with large density
differences by using the projection method. In this method, two particle velocity
distribution functions are used. One is used for calculating the order parameter to track
the interface between two different fluids; the other is for calculating the predicted
velocity field without pressure gradient. The corrected velocity satisfying the continuity
equation can be obtained by solving a Poisson equation.
On the other hand, liquid droplet dynamics in a PEMFC gas channel is the result of
shear and drag forces due to the gas flow and the capillary forces associated with the
liquid-gas-solid interface. A finite steady-state equilibrium or static contact angle, also
known as partial wetting, can be reached due to the balance of surface tension forces.
Based on these considerations, Briant et al. [73] developed an approach based on the free-
energy LBM model introduced by Swift et al. [69] to simulate partial wetting and contact
line motion in single or two-component, two-phase fluids. Unfortunately, this method has
inherently the disadvantage of the original free energy LBM model of Swift et al. and can
only be used to simulate two-phase problems with a small density ratio. The maximum
density ratio in the simulations of droplets on partial wetting surfaces was reported to be
just around 2. In order to simulate a flow of two-phase fluids with a large density ratio on
a partial wetting wall, a new LBM scheme of the LBM is required, for example the
approach proposed by Yan and Zu [71]. The multiple-relaxation time LBM method has
 24
also been proposed to overcome to density ratio limitation, but to date has not been
demonstrated for high density channel flows, but the results for air-water flow in gas
diffusion layers are promising [74].

2.5.5. Modeling of the dynamic contact angle

Volume of fluid (VOF) CFD simulations performed to date by our group and others to
investigate the dynamics of liquid water droplets in hydrophobic micro-channels do not
rigorously account for the moving contact line mechanisms, which in any case are not
fully understood. Noting that surface tension is dominant (Capillary number Ca <<1) in
the flows of interest, the contact line singularity and the dynamic contact angle will be
handled by a suitable semi-empirical formulation in conjunction with a generalization to
3D of a recent dynamic contact angle algorithm [75].
Given the prominent role of surface tension in determining the force balance on water
droplets, a resolution of the dynamic effects on the contact line is critical to achieving
physically representative CFD simulations. The three-phase contact line presents
challenging theoretical and numerical problem in the context of continuum simulations in
which a no slip condition is usually applied at solid boundaries and walls. This gives rise
to a singularity (infinite tangential stress) when considering the moving contact line.
Another issue is the modeling/tracking of the behavior of the dynamic contact angle,
which is dependent on the flow conditions around the droplet as well as the contact line
motion.
One approach of dealing with the singularity is to impose a slip condition at the
boundary. Some of the earlier models proposed in this context are summarized in Table
2.1 based on the assessment of provided in Shikhmurzaev [76]. The main weakness in
most of the models is in dealing with the stress singularity. Blake and Shikhmurzaev [77]
presented a theoretical approach based on non-equilibrium thermodynamics to deal with
this issue and derived the following model for the dynamic contact angle:
2V (  2s*e  1se*uo )
cos S  cos D  (2.3)
(1  1se* )[(  2s*e  V 2 )1 / 2  V ]
where V is the dimensionless contact line velocity, and the non-dimensionalized local
densities  ies*  ( ies /  0s ) are defined for each phase in reference to their equilibrium
 25
values, with indices 1 and 2 referring to the free surface and the liquid-solid interface
respectively.

Table 2.1. Slip models examined by Shikhmurzaev [76]

Author(s) Model Basis Drawbacks
Huh & Mason [78] Physical model of the Does not display rolling
liquid-gas interface motion motion. Leads to integrable
singularity at the solid
boundary.
Durbin [79] Slip prescribed by Fair but not exact
bounding maximum shear agreement with experiments
stress density
Baiocci & Pukhnachev [80] Problems with one-sided Fair but not exact
constraints for agreement with experiments
Navier-Stokes equations
and the dynamic contact
angle

The dimensionless contact line velocity is given by [77]:


V U (2.4)
  1  4 
s
0

and the parameter u 0 is defined as:

sin  d   d cos d
u0  (2.5)
sin  d cos d   d
The other parameters used found in the model are phenomenological coefficients,
 h
α = Effect of surface tension gradient on the velocity distribution    
 

 
β = Effect of shear stress on the velocity distribution    
 h

 5Sc 2 h 
τ = Surface tension relaxation time   
  1e (1  1se* ) 
 26

To date no documented attempt has been made in implementing this model in CFD
codes, probably because lack of numerical robustness. Empirical or semi-empirical
models appear for now to be more practical. A good overview of more recent semi-
empirical models used to account for the dynamic contact angle is provided by Sikalo et
al. [81]. Most of the models assume Young’s equation is valid throughout the dynamic
process and the solid-liquid and solid-vapor surface tensions vary with flow field
dynamics, i.e.
 cos  equilibrium   solid liquid   solid vapour (2.6)

A key aspect in the definition of the Capillary number (Ca) used in the models is the use
of the contact line velocity as the velocity scale, i.e. Ca = (Vel)µ/ . Thus the formulation
proposed by Cox [82] reads:

1 3
 Q 
 
Ca   ln  1 
Q1
 2  g  D   g  e   O 1 1 

 f  D  f  e    
 ln  
(2.7)

Where ε is a dimensionless parameter based on the static contact angle mechanics. Q1 and
Q2 are parameters based on the outer flow field and the slip conditions on the wall
respectively and the functions f and g depend on the dynamic and equilibrium contact
angles.
Another popular contact angle formula, which again uses a capillary number based on the
contact line velocity, is the Hoffman-Voinov-Tanner law [81]:

 D 3   e 3  Ct Ca with Ct  72 (2.8)

This formula also involves a capillary number based on the contact line velocity:

Ca 
Vel  
 (2.9)
 27
More recently, a formulation with potentially broader applicability was proposed based
on the Hoffman functions [81]:

 D  f Hoff [Ca  f Hoff

1
( e )] (2.10)

   x   
0.706

f Hoff ( x)  arccos1  2 tanh5.16 0.99   (2.11)

   1  1.31x   

This formulation provides a good fit with experimental data and was implemented
recently with some success in both LBM [83] and VOF [51] simulations.
It is particularly important to understand the surface tension force in predicting two-
phase flow in a PEMFC microchannel. The main difference between level set methods
and the VOF method is that the VOF method uses a discontinuous function (zero in one
phase and one in the other) while level set methods represent the interface by a certain
contour of a smooth function. Each of these methods has one main advantage and
disadvantage. Due to the discontinuity at the air/water interface, the VOF method can
suffer from poor accuracy in determining the position of the interface and the calculation
of the mean curvature, which in turn determines the surface tension force. The main
advantage of the VOF method is that the mass of each fluid is exactly conserved. Level
set methods, on the other hand, use a smooth function, allowing higher accuracy in
resolving the interface. However, one major drawback of the originally proposed level set
method is that mass is not conserved, and significant mass losses may occur. A hybrid
numerical method with a high order of accuracy, and good mass conservation properties
would be ideal. In this aspect, Sussman and Ohta [84] successfully developed LSM in
conjunction with VOF for treating surface tension in incompressible two-phase flow.
It should be emphasized that VOF simulations conducted to date for geometries
relevant to PEMFC channels, do not incorporate the physics of the dynamic contact line
which is expected to alter significantly the detachment and subsequent evolution of water
droplets. The robustness of the VOF methodology is due in some measure to the relative
simplicity of the piecewise linear interface calculation (PLIC) procedure used to
reconstruct the liquid-gas interface topology. The quality of the approximation impacts
 28
the computation of the surface curvature and hence surface tension forces. In this
respect, the contact angle is a critical parameter in evaluating the surface adhesive force.
VOF simulations conducted to date for geometries relevant to PEMFC channels, do not
incorporate the physics of the dynamic contact line which is expected to alter
significantly the detachment and subsequent evolution of water droplets, but rather
prescribe a static contact angle based on Young’s equation. Dynamic contact angle
models using the semi-empirical relations discussed in Section 2.5.5 do improve the
physics significantly as shown in recent work by Fang et al. [34] for two-phase slug
flows, and recent results obtained by Miller [51] for a model cathode flow, which
highlights the large impact the contact line prescription has on droplet dynamics.
Dynamic implementation, whether theoretical or empirical (through a series of static
angle corrections at each time step in the simulation) impacts the numerical convergence,
in order to ensure practical simulation times, this will probably necessitate the
implementation of alternative algorithms.
 29

Chapter 3

3 Experimental Investigation of Water

Droplet Emergence
The focus of this chapter is the experimental investigation of the water droplet
emergence process in a model PEMFC microchannel. Strategies for effective water
management rely on an understanding of two-phase flow in gas channels in PEMFCs
under realistic operating conditions. Since the primary source of the liquid phase in the
gas channel is the water droplet transfer and formation through the GDL, the analysis of
water droplet development, evolution and instability behavior on the GDL surface is
essential. Building on our preliminary experiments [85] and previous computational work
[36], we present an experimental characterization of the emergence, detachment and
subsequent dynamic evolution of water droplets in laboratory model of a fuel cell cathode
channel. The controlled ex-situ experiments are analyzed using flow visualization and
image analysis and provide further insight into the droplet dynamics process and
quantitative data for the validation of numerical models.

3.1. Method and Apparatus

3.1.1. Microchannel design

A specifically designed microfluidics platform was fabricated for the experimental
study of droplet dynamics. The water pore size and shape as well as the dimensions and
 30
surface properties of the microchannel must be considered in the design in order to
model the phenomena of water emergence relevant to fuel cell operation. The
microfluidic chip produced for this specific purpose (see Appendix A for detailed
fabrication procedure) includes a single channel with a cross section of 250 µm × 250 µm
for air flow, as well as a single water injection pore with a cross section of 50 µm × 50
µm located in the middle of the air channel. An elastomer, polydimethylsiloxane (PDMS,
Dow Corning Sylgard 184), was selected as the casting material not only due to its
transparent and optical clear characteristic but also because its surface is hydrophobic and
provides a static contact angle of ~110°, a value similar to that of carbon paper without
PTFE coating.
For an actual GDL, pore sizes and locations are randomly distributed thus making it a
challenge to predict and control the locations of water droplets in the fuel cells. However,
the size of the “pores” is of the order of 50 ~ 100 µm [10]. Experimentally, the
emergence pore location can be prescribed on a GDL as was also shown by Kimball et al.
[86] by creating 250 µm pores at specific locations. There is also evidence that
preferential pathways do form spontaneously in standard treated GDLs [87], probably as
a result of residual water. In the present study, the 50 µm square size of water pore is
prescribed to represent such conditions.
A complete chip is shown in Figure 3.1a. Viewing from the side, Figure 3.1b, allows
observation of an emerging droplet as illustrated. The total length of the microchannel is
around 37 mm, whereas the field of view of the microscopic images only covers 3 mm
(Figure 3.1c). The dimensions of the channel cross section are the same as reported by
Zhu et al. [36], and are suitable for representing a typical gas flow channel used in micro
PEM fuel cells [88]; the pore size for water emergence is close to the mean values for
carbon paper reported in the literature [10][89]. Keeping in mind that the intended use of
the data for validation of various numerical methods [36][67][41], a square geometry was
selected as it is computationally more convenient to reproduce without requiring grid
skewing or unacceptable grid aspect ratios. The experimental set-up exhibits the most
salient physical mechanisms and maintains the appropriate range of values for the key
non-dimensional parameters.
 31
Since gravity is perpendicular to the direction of water emergence, it is necessary to
clarify the potential gravitational effects. The Bond number, denoted as the ratio of the
gravitational force to the surface tension force, is defined as

gL2
Bo  (3.1)

where  is the density of water droplet, g is the gravitational acceleration, L is the
characteristic length of the system (the measured maximum droplet height is around 200
µm), and  is the surface tension force between water and air. The Bond number in the
current system is estimated to be less than 0.0054 indicating that surface tension forces
dominate and the effects of gravity are negligible.

Figure 3.1. (a) PDMS chip for droplet manipulation. (b) Cross sectional view of chip.
(c) Field of view in microscope.

3.1.2. Measurement apparatus

A schematic diagram of the experimental apparatus is shown in Figure 3.2. This setup
allows injection of water from a small pore that emerges into an air stream while
simultaneously monitoring the pressure drop across the gas channel (P), the injected
water volume (Qw), and the air flow rate (Qa). A side view visualization of two-phase
flow is obtained as droplets move and interact with the gas stream inside the
microchannel. The data acquisition system was developed using LabView 8.5 (National
 32
Instruments) to integrate the water syringe pump (Harvard Apparatus, PHD 2000), air
mass flow meter (Omega, FMA-1600A), differential pressure sensor (Honeywell,
FP2000), and high speed camera (Vision Research, Phantom MIRO-4). Data streams are
initialized by the high speed camera to trigger the image sampling process. Eight-bit,
gray scale images are captured as well as time-resolved pressure drops, delivered water
volumes, air flow rates, and inlet upstream pressures and temperatures. These can be
measured simultaneously allowing correlation of droplet dynamics with flow regimes.
The collected data provide initial conditions, making the experiments suitable for
numerical validation. Digital image processing was performed using the toolbox in
Matlab 7.8. Images were recorded through an inverted microscope (Zeiss, Axiovert
200M) and magnified by a 5x objective lens (Zeiss, EC Plan Neofluar). Bright-field
illumination using a white light LED was applied throughout the image acquisition
process.

Figure 3.2. Schematic diagram of experimental apparatus.

The high speed camera consists of a CMOS (complementary metal oxide

semiconductor) image sensor with a maximum of 800 × 600 pixels. It has a capability of
12 bit pixel depth and a maximum full-resolution frame rate of 1000 Hz. However, in
 33
order to resolve the droplet dynamics, the images were adjusted to covering only
800×72 pixels. This boosts the camera frame rate up to 6006 Hz and provides around six
seconds of recording time with the 2GB built-in memory in the camera. Time resolved
images were fed into the contact angle measurement software (DataPhysics Instruments,
SCA202 v3.51.3 build 187). The tangent method for sessile drop measurement that fit the
three-phase point where the liquid touches the solid surface was used. A third degree
rational function was used to evaluate the advancing and receding contact angles as well
as to measure the droplet height and chord length. The accuracy of the contact angle
analysis depends upon the image processing technique and parameters used to extract the
drop edge from images. Droplet images are converted into 8 bit gray scale by
thresholding, resulting in images consisting of a dark foreground representing the drop,
and a white background. When moving from the background into the droplet, a sharp
decline in pixel intensity is observed, and the edge of the drop is expected to be in a
region along the curve (spline) corresponding to the highest gradient. Using a maximum
of 80 points to define the interface profile, a least square fit is used to determine the
tangent at the advancing and receding contact points. The raw image from the high speed
camera provides a resolution of 4.54 µm/pix. Due to the small droplet dimension (100-
350 µm), the raw images are imported into the contact angle evaluation program after a
three-fold magnification. Combining the uncertainties result in a maximum uncertainty of
±1.3 of contact angle evaluation.

3.1.3. Flow conditions

The flow channel described in Section 3.1.1 (250 µm × 250 µm cross section and 37
mm length) nominally simulate an operating fuel cell active area (Aact) of 0.0925 cm-2.
The volumetric flow at the cathode due to water production and anode-to-cathode
transport is given by:

1  2 iAact M H O
Qw  2
(3.2)
2F
Where α is the net drag coefficient, Aact is the active area, i is the current density, M H 2O

is the molecular weight of water, and F is Faraday’s constant. For operation at i = 2.0 A
 34
−2 -1 -2
cm with α = 1, the water flux is 33.6 µL min cm . The injection rates used in this
study range from 32.4 to 162.2 µL min-1 cm-2. The higher water flow rates were
specifically chosen to allow us to conduct the experiments in a reasonable time and
without incurring secondary effects such as evaporation. Appropriate normalization of
the time scale and the negligible impact of water injection rates on contact angles and
droplet dimensions discussed in the results below support this. The Reynolds numbers for
the air flow range from approximately 50 ~ 1200, representative of a range of operating
conditions for parallel and serpentine PEMFC configurations. The inlet flow conditions
used in the experiment are summarized in Table 3.1.

Table 3.1. Flow inlet conditions.

Test matrix for flow regime

Air flow rate, Qa (sccm) 12 ~ 300
Air superficial velocity, Va (m s-1) 3.2 ~ 80
-1
Water flow rate, Qw (µL min ) 3, 6, 15
Water superficial velocity, Vw (m s-1) 0.02, 0.04, 0.10
Reynolds number of air, Re 51 ~ 1268
Temperature, T (K) 296
Two specific cases for qualitative droplet dynamic analysis, Case 1: Va = 10 m s-1, Vw =
0.02 m s-1 and Case 2: Va = 10 m s-1, Vw = 0.04 m s-1.

3.2. Results and Discussions

The flow visualization experiments were analyzed both qualitatively, to develop a
flow regime map, and quantitatively, to determine the evolution of the dynamic contact
angle, the characteristic detachment frequency and droplet size.
 35

3.2.1. Flow regimes

In order to investigate the water droplet dynamics subjected to a shear flow in the
microchannel corresponding to a range of fuel cell operating conditions, three water
superficial velocity conditions were investigated. The superficial velocity is defined as
the bulk velocity of the water flowing through the cross sectional area of the channel.
Visualization of the water emergence process could be typically categorized into three
distinct flow regimes shown in Figure 3.3:

Figure 3.3. Typical flow regime in microchannel

1. Slug flow (Figure 3.3a) occurring at low air superficial velocities. The drag force is

minimal at these velocities, and the droplets keep growing until they contact the

channel wall and are then either slowly convected or remain on the wall until caught

up by another droplet with which they coalesce and then move along the channel

toward the outlet. Slug flow obstructs and limits air flow through the channel and thus
 36
increases the possibility of flooding and decreases the water removal effectiveness

in the channel.

2. Droplet flow (Figure 3.3b) occurring at increasing air velocity, with an individual

droplet emerging from the water pore and remaining pinned at the water pore due to

surface tension forces, until it grows to a critical size at which the air drag force

overcomes the surface tension force. In this regime, the effect of air velocity on the

changing shape of the emerging droplet and the time evolution of the emerging

process can be clearly identified.

3. Film flow (Figure 3.3c) that ensues with further increase in air flow velocities that drag

the droplets and induce the formation of an elongated and almost continuous wavy

water film. Large droplets are no longer observed. Instead, smaller protrusions appear

at the boundary of the film and are rapidly flattened due to a high air flow rate.

In summary, the water droplets can be shearing off the PDMS surface with a higher air
flow velocity, causing only a minimal obstruction. This occurs most often while droplets
are small compared with the channel geometry. Upon detachment, the droplet can interact
with the different channel walls. The detached water droplets can adhere to the surfaces,
coalesce and build up into water lumps, films or slugs, blocking the flow passage over
time. The different coalesced structures of water inside the flow channel will produce
different flow instabilities.
The different flow regimes are typically presented in two dimensional flow pattern
maps, with the two coordinates representing some appropriate hydrodynamic parameters.
The most popular coordinates are superficial phase velocities, but dimensionless Weber
numbers have also been used, as suggested by Akbar et al. [90]. Both representations are
shown in Figure 3.4. Dashed lines are used to delineate the three different regimes based
on the testing conditions, and provide the threshold air velocities for the transitions from
slug flow to droplet flow and then to film flow at a given water velocity. With increasing
 37
water velocity, the droplet flow regime region becomes smaller, and at some point the
slug and film regimes intersect, with transition occurring directly and bypassing the
droplet flow regime.

Figure 3.4. Flow map of water emergence phenomena in a model PEMFC cathode
gas microchannel.
 38

3.2.2. Droplet emergence frequency

The flow visualization presented above indicates that the periodicity of the droplet
emergence process is specific to the flow regime. Time series analysis was performed to
further characterize the droplet dynamics using digital imaging processing. Time series
raw images are first subtracted from a reference background image taken in the absence
of droplets. All images are cropped to the region of interest covering a single droplet in a
preselected registered image. This registered image is then compared with other images
by using a MATLAB image processing and signal analysis toolbox. Preselected
registered image A is compared with every other image B in the same recording series.
An algorithm to determine the two-dimensional correlation coefficient index, r, is applied
based on

 ( A mn  A )( Bmn  B )
r m n
(3.3)
 2  2
  Amn  A    Bmn  B  
 m n  m n 

where Amn and Bmn are m  n matrices, and A and B are the mean matrix elements. The
correlation coefficient values range from 0 to 1. When using the algorithm to compute
time resolved image correlations, a few images are randomly selected to serve as the
register image A, and are compared with image B in the image stack. The two
experimental flow conditions listed in Table 3.1 were selected to demonstrate the
concept.
A typical time domain signal for Case 2 is shown in Figure 3.5a and the corresponding
frequency spectrum obtained using a Fast Fourier Transform (FFT) analysis is shown in
Figure 3.5b, clearly identifying the periodicity. In case 2, the droplet emerges at a
dominant frequency of 13.2 Hz. A sample MATLAB source code for determining the
correlation coefficient index and emergence frequency is given in appendix B.
 39

Figure 3.5. (a) Time domain signal and (b) frequency distribution of droplet
emergence process.

The evolution of the dominant emergence frequency with flow conditions is presented
in Figure 3.6. The low emergence frequency at low air velocities allows droplets to grow
sufficiently large to contact the side walls. The frequency increases linearly with air
 40
velocity and reaches its highest value before transitions to film flow. No dominant
frequency is observed in the film flow regime, where the process is inherently a periodic
due changes in detachment frequency and necking, splitting and coalescence of the films.
The volume of a droplet can be further deduced from the water flow rate and the
emergence frequency. In the current studies, the droplet volumes reached up to the onset
of detachment range from 4.2 to 10 nL, 1.1 to 10 nL and 0.7 to 2.5 nL for water flow
rates of 0.02, 0.04 and 0.1 µL min-1, respectively. Interestingly, this suggests that this
chip platform is promising for manipulation of a nanoliter droplet.

Figure 3.6. Emergence frequency in droplet flow regime under different flow
conditions.
 41

3.2.3. Further image analysis of droplet emergence

Figure 3.7 displays the time series images captured at 6006 Hz illustrating the
evolution and cyclic nature of the water emergence into the gas stream. Only one cycle of
emergence is shown corresponding to Case 2. The emergence frequency is 13.2 Hz (with
a period of 75.75 ms), as presented earlier. The right boundary of the droplet is pinned at
the water outlet where the air, water and the PDMS surface meet. The dynamics of a
droplet is governed by the balance of pressure, shear and drag forces exerted on it.
Surface forces dominate the beginning of the growth and emergence cycle (0 to 10 ms)
shown in Figure 3.7a-d during which the droplet retains a symmetric shape.

Figure 3.7. Time resolved images of water emerging from a 50 µm square pore in
a 250 µm square gas microchannel with the flow condition of Case 2. a) t = 1 ms,
b) t = 3 ms, c) t = 5 ms, d) t = 10 ms, e) t = 15 ms, f) t = 20 ms, g) t = 25 ms, h) t
= 45 ms, i) t = 65 ms, j) t = 75 ms

Starting with Figure 3.7e, the larger droplet produces increasing blockage of the air flow
that results in increasing drag force inducing deformation with tilting toward the
downstream direction and the onset of contact angle hysteresis. However, the surface
adhesion forces still dominate drag and the droplet remains pinned. For these flow
conditions, the droplet does not grow sufficiently large to contact the top and side walls,
but eventually detaches (Figure 3.7j) and starts to move along the surface while another
emergence cycle begins.
 42

3.2.4. Dynamic contact angle

The high speed visualization images were further analyzed frame by frame to deduce
the evolution of the droplet chord and height and of the dynamic contact angles. The
advancing (a, upstream de-wetting side) and receding (r, downstream wetting side)
contact angles of a droplet vary from the static value as shown schematically in Figure
3.8a, and their evolution during the emergence process up to the point of detachment is
presented in Figure 3.8b for Case 2. Three different periods can be identified. Period I (0
~ 11 ms: surface tension dominated) corresponds to the beginning of emergence. In this
region, a and r increase at the same rate which corresponds to a symmetric change of
the droplet shape and thus no contact angle hysteresis is observed. In Period II (11 ms ~
17 ms: transition region), a keeps increasing but r stops increasing. This occurs as a
result of pressure forces acting on the upstream surface of the droplet and countering the
increase in r due to droplet growth. In Period III (17 ms ~ 75.75 ms) increasing pressure
and shear forces combined with growth induce further increase in the advancing contact
angle a up to 110º. During this process, significant distortion of the droplet occurs.
Conversely, the receding contact angle becomes smaller as the cycle progresses up to
detachment. The small variations in measurements around the trend lines are due to local
variations in the PDMS surface of each batch of fabricated chips. Additionally the
experimental surfaces are not perfectly homogeneous and occasional slip, jump and stick
motion of the contact line were observed in flow visualization using a high speed camera.
These phenomena are expected to be more pronounced in actual GDL, in which is likely
to play an even more prominent role.
The chord length of a growing droplet increases and reaches about 350 µm at
detachment which is larger than the channel dimension of 250 µm. This underscores the
significant deformation and elongation due to the hydrodynamic forces exerted by the air
stream. The assumption of a spherical droplet in the simplified models for predicting the
onset of droplet instability [7] are thus only very approximate. The data presented here
for the time evolution of the dynamic contact angles can be used to empirically prescribe
the dynamic contact angle in order to improve the physical realism of simulation.
 43

(a)

(b)
Figure 3.8. Dynamic contact angle evolution through an emergence cycle (13.2 Hz
under Case 2 flow conditions). Period I: surface tension force dominant. Period II:
transition. Period III: drag force dominant.

3.2.5. Droplet dynamic at the onset of detachment

When water is introduced into the air stream, the droplet initially emerges, grows and
remains pinned at the pore until it attains a critical size and detaches. The detachment
droplet size is a function of several parameters, including water and air flow rates surface
conditions, and pore size [55]. The droplet chord length, height, and the advancing and
receding contact angles can be deduced from the flow visualization at the onset of
 44
detachment. The dependence of the detachment droplet chord (C) and height (H) on air
velocity/Reynolds number is plotted in Figure 3.9 for various water injection velocities.

Figure 3.9. Effect of air flow velocity on characteristic droplet size (chord C and
height H) at detachment.

The air velocities can be interpreted as the critical operating velocities which must be
used in a PEMFC to ensure detachment and removal of liquid water to prevent flooding
and the subsequent blockage of the transport of reactants to the reaction sites [2]. The
windows in the figure delineate the size of droplet flow regime for three different water
injection velocities. The left and right hand sides of the window represent slug and film
flow regimes, respectively. The low velocity boundary of the windows at around Va = 8
m s-1 can be interpreted as the threshold air velocity that prevents the formation of slug
flow which is characterized by larger droplets and blockage. Slower water injection
velocities broaden the droplet flow regime window. The figure also shows a similar trend
to the numerical data presented by Zhu et al. [56] with the rate of size decrease tapering
at higher velocities.
The effect of the air flow on contact angle hysteresis ( = a - r) is plotted in Figure
3.10. This angle could be interpreted as the measure of the ability of the droplet to resist
the drag force and the capability of the flat PDMS surface to remove water under
 45
controlled experimental conditions. The data for different water velocities collapses to
the same curve and decreases linearly from 80° to 40°. As the air flow velocity increases,
the flow regime changes from a droplet pattern to a film pattern with a smaller advancing
contact angle and a reduction in the contact angle hysteresis.

Figure 3.10. Effect of air flow on contact angle hysteresis.

The analysis of Kumbur et al. [9] indicated that the contact angle hysteresis depends
on the channel air flow rate (or Reynolds number), droplet size (H and C) and surface
properties. The plot of the contact angle hysteresis versus the H/C ratio at different
Reynolds numbers is shown in Figure 3.11a. The data differs from that predicted by
Kumbur et al. [9]. However, in their analysis, in addition to different flow conditions, the
droplet was assumed to be hemispherical and symmetric and the effect of pore pinning
was not accounted for. The importance of accounting for pore connectivity in
determining the critical detachment size was highlighted by Zhu et al. [36]. The data
obtained here reflects the distribution on a flat PDMS surface under flow conditions
described in the test matrix of experiments.
 46

(a) (b)

Figure 3.11. Contact angle interpretations and effect of airflow on droplet aspect ratio
at the onset of detachment.

For a given droplet size (or aspect ratio, H/C), a higher contact angle hysteresis will
increase the capability of the droplet to resist the drag force. The droplet will thus have a
broader stability zone as shown by plotting the data appropriately in Figure 3.11a. In fuel
cell applications, a broader instability zone is desirable to enhance water removal
capabilities. Figure 3.11b shows the critical droplet size in terms of the aspect ratio versus
the Reynolds number. A higher Reynolds numbers at a given droplet size increases
instability and potential for removal of the droplet. The Reynolds number required to
sweep out a droplet is found to increase with decreasing droplet aspect ratio. Spreading of
the droplet, as discussed in [9], results in a decrease of the aspect ratio and requires
higher air velocities/ Reynolds number for detachment. As shown in the flow
visualization, a decrease in the aspect ratio occurs when the flow regime switches to a
film flow pattern. This flow regime makes water removal more difficult in fuel cell
applications. Another parameter that impacts the droplet detachment process is pore
diameter. Varying pore size was outside the scope of the present experiments, but
numerical simulations [55] indicate that although the aspect ratio decreases with
decreasing pore diameter, the critical Reynolds number for droplet removal decreases as
 47
a result of weaker surface tension due to the reduced connectivity with the pore. The
role of pore connectivity is clearly important and requires more in-depth analysis.

3.3. Summary
The dynamics of water droplets emerging from a pore in the presence of a cross flow
of air was investigated experimentally using a modeled PEMFC cathode gas channel.
Quantitative analysis of high speed flow visualizations was performed to deduce
characteristic process frequencies, contact angle hysteresis and critical droplet sizes. The
study shows

1. Three flow regimes: slug, droplet and film flow patterns, are identified under different

air and water velocities. At low air velocities, slug flow blocks the air flow through the

channel. At higher air velocities, a periodic pattern of droplet emergence, growth and

detachment appears. Further increase in air velocity induce wavy water film pattern. A

flow map of the flow regimes as a function of superficial air and water velocities was

presented.

2. For emerged droplets, significantly higher critical air velocities are observed compared

to results in the literature which only considered processes starting from droplet

initially static on a flat surface. This highlights the important impact of pore

connectivity.

3. The dynamic contact angle at the onset of detachment decreases as the air velocity

increases; as a result, the flow regime shifts from slug to droplet and film flow. This

angle is representative of droplet stability, i.e. its ability to resist the drag force on a

given surface.
 48
4. A decrease in the droplet aspect ratio and contact angle hysteresis is observed in the

film flow regime, and is found to reduce water removal capacity.

While the experiments presented here were obtained using a laboratory model of a
PEMFC cathode, they isolate and represent some of the salient flow features of operating
fuel cells, and the data and documented boundary conditions should prove useful for the
validation of simulation methods. Two key characteristic of GDLs that will be considered
in future work are roughness and inhomogeneities. These are expected to impact the
effective contact angle, induce additional pinning and alter the droplet growth and
dynamics. Finally, the critical impact of pore connectivity and pinning, which was also
qualitatively reproduced in earlier numerical work [36], needs to be incorporated into
practical force-balance models that can be used in design.
 49

Chapter 4

4 Numerical Simulation Using VOF

Method
The simulation of two-phase flow presents many challenges, particularly for
conditions corresponding to an operating fuel cell. To complement the experiments,
numerical simulations of the droplet emergence using VOF based free surface capture
method were investigated. VOF simulations have been presented for a similar problem as
noted in Chapter 1. However to date they relied on a static prescription of the contact
angle. The work in this chapter focuses primarily on the implementation of a dynamic
contact angle model. Two different approaches were implemented for the treatment of
dynamic contact angle of a moving contact line using commercial Computational Fluid
Dynamics (CFD) software, CFD-ACE+ and ANSYS FLUENT. Both are well established
and widely used for research and development and incorporate state of the art numerical
algorithm including a well-documented VOF module. The first approach used for the
dynamic contact angle model is through the velocity dependent contact angle function
derived from the experiments, and was implemented in the computational framework of
free surface module using CFD-ACE+ version 2010. The second approach is by
implementing the Hoffman function (cf. Eq. (2.11)) within the user defined function
framework of the multiphase flow module using ANSYS FLUENT version 14.0.
 50

4.1. Volume of Fluid Method

The VOF method is a computationally intensive approach first developed in the early
1980s by Hirt and Nichols [52] and refined continuously by many others. The VOF
model can model two or more immiscible fluids by solving a single set of momentum
equations and by tracking the volume fraction of each of the fluids through the solution
of an additional volume fraction transport equation. In each control volume, the volume
fractions of all phases sum to unity. The fields for all variables and properties are shared
by the phases and represent volume-averaged values, as long as the volume fraction of
each of the phases is known at each location. Thus the variables and properties in any
given cell are either purely representative of one of the phases, or representative of a
mixture of the phases, depending upon the volume fraction values [91]. The volume
fraction of fluid q in the cell is denoted as αq, where:
 q  0 : The cell is empty of the q th fluid.
 q  1 : The cell is full of the q th fluid. (4.1)
0   q  1 : The cell has an interface containing the q th fluid.

Based on the local value of αq, the appropriate properties and variables are assigned to
each control volume within the domain. The tracking of the interface between the phases
is accomplished by the solution of a continuity equation for the volume fraction of one
(or more) of the phases. For the qth phase, this equation has the following form:

t

 q  q      q  q Vq  0 (4.2)

in the absence of source terms and mass transfer between phases. The volume fraction
equation is not being solved for the primary phase. Instead, it is computed through
explicit time discretization based on the following constraint.
n


q 1
q 1 (4.3)

Using the explicit approach, finite-difference interpolation schemes are applied to the
volume fractions that were computed at the previous time step.
 qn 1   qn
t
 
V   U nf qn, f  0 (4.4)
f
 51
where
n+1 : index for new (current) time step
n : index for previous time step
αq,f : face value of the qth volume fraction, computed from second order upwind scheme
V : volume of cell
Uf : volume flux through the face based on normal velocity.
The properties appearing in the transport equations are determined by the component
phases in each control volume. In general, for n-phase system, the volume-fraction-
averaged composition variable takes the following form:
n
B    q Bq (4.5)
q 1

where B is a fluid property such as density or viscosity. This allows the calculation and
updates of the entire set of fluid properties in any given cell, based on the phases that are
present in that cell. A single momentum equation is solved throughout the domain, and
the resulting velocity field is shared among the phases. The momentum equation, shown
below, is dependent on the volume fractions of all phases through the properties ρ and µ
[91].

t
   
V    V V  P       V  V
 
T
   g  F


(4.6)

where P is the static pressure, F is body force term, ρ and µ are the volume averaged
density and dynamic viscosity which can be computed using Eq.(4.5).
In the present study of microchannel, due to the relatively small mass of water
droplets, body forces are negligible compare to surface tension forces. For a flow
channel, with droplet dimensions in the order of less than 0.5 mm and a water velocity
less than 0.04 m/s, the Weber number is:
U 2 L Inertia
We 
 Surface Tension
(4.7)

We 
998 0.04  0.0005  0.011
2

0.0728
 52
From experimental observation, the typical dimension of a water droplet is less than
0.5 mm even under extreme operating conditions. Hence, the surface tension force is
dominant and approximately 100 times the strength of the inertia force. Further
decreasing of water injection velocity will result in smaller droplets and an even smaller
Webber number.
Surface tension is accounted for by using the continuum surface force (CSF) model
[92] to represent the pressure jump induced by the surface tension within the transition
region. This pressure jump depends on the surface tension coefficient and the surface
curvature as measured by two radii in orthogonal directions, R1 and R2:
1 1 
p2  p2      (4.8)
 R1 R2 
where p1 and p2 are the pressures in the two fluids on either side of the interface. The
surface curvature is computed from local gradients of the surface normal at the interface
[91].
n   q (4.9)

where n is the surface normal, defined as the gradient of the volume fraction of the qth
phase. The curvature, k, is expressed by the divergence of the unit normal, n̂ .
k    nˆ (4.10)
where
n
nˆ  (4.11)
n

A continuous transition is implemented in the momentum equation as a volume force

source term, Fvol, with the following form when only two phases are considered:
ki  i
Fvol   ij (4.12)
1
2
i   j 
where ρ is the volume-averaged density computed using Eq. (4.5). This implies that the
surface tension source term for a cell is proportional to the average density in the cell. In
conjunction with the surface tension model while using wall contact angle as a boundary
condition, the surface normal is expressed in terms of the prescribed contact angle, w.
nˆ  nˆ w cos  w  tˆw sin  w (4.13)
 53
where n̂w and tˆw are the unit vectors normal and tangential to the wall, respectively. It
is worth noting that the combination of this contact angle with the calculated surface
normal one cell away from the wall determine the local curvature of the surface, and this
curvature is used to adjust the body force term in the surface tension calculation.
Therefore, suitable mesh refinement is necessary at the near wall location.
The geometric interface reconstruction is obtained from a piecewise linear interface
calculation (PLIC) method [93]. This assumes that the interface between two fluids has a
linear slope within each cell having a volume fraction greater than zero and less than one;
this linear shape is then used for the calculation of the advection of fluid through the cell
faces. A faithful determination of the interface normal vector is a critical step in the PLIC
method, since it determines the slope of the interface as well as the fluid volume
matching the given volume fraction for the cell.
In performing a time-dependent VOF calculation using an explicit scheme, the time
step used for the volume fraction calculation is not necessarily the same as the time step
used for the rest of the transport equations. A dimensionless Courant number, Co, is
introduced which compares the time step in a calculation to the characteristic time of
transit of a fluid element across a control volume.
u fluid t
Co  (4.14)
xcell
The Courant number is a measure of the transport distance crossed over a time step as a
fraction of the local length scale of the cell.

For instance, a Courant number 0.2 in a VOF computation would enforce a time step that
allows the fluid or the interface to cross up to 20 percent of the width of a cell during
each time-step, regardless of the actual velocity or cell dimension [94]. In order to
minimize calculation time while ensuring stability, the Courant number must formally
satisfy only Co < 1. However, in practice, values as low as 0.1 are often required. The
length of the next time step is calculated based on a fixed local Courant number in the
cells. The smallest time step is used as the characteristic time, which represents the
maximum velocity found in a given group of cells. Therefore, the time step size for a
 54
group of cells will be restricted to ensure that the free surface crosses less than a cell
during that time step.
The VOF free surface module in CFD-ACE+ provides several additional features for
solver control to remove tiny unphysical isolated droplets of liquid in gas region or tiny
isolated gas bubbles in liquid region [94]. These droplets and bubbles are collectively
called flotsam and jetsam, and are due to inadequate convergence, improper Currant
number or excessive skewness in the grid. Though, this can be minimized with the PLIC
reconstruction, it cannot be avoided entirely. Another feature of CFD-ACE+ is the
dynamic contact angle numerical implementation that can be achieved as an option by
simply introducing an user define function rather than having to implement this through
an user define subroutines in FLUENT.

4.2. Implementation of Dynamic Contact Angle

In the VOF method, a description of the contact angle is essential. The contact angle at
each location represents an interface or a boundary between a solid and a fluid zone. A
fixed value static contact angle (SCA, s) provides the simplest option. The value is often
pre-specified and remains unchanged regardless of changes in the flow or other
conditions. This is a suitable option if a quasi-steady or equilibrium solution is being
considered. In the absence of knowledge about the variation of the contact angle with the
flow or other conditions, applying the static value is the most expedient, but not
necessarily accurate approach. In the present study, however, the droplet emergence
phenomenon in the cathode gas channel is rarely quasi-steady. Instead, droplets emerge
randomly from the rough surface of the gas diffusion layer and interact dynamically with
the gas stream. Hence, treatments of dynamic contact angle using two different modeling
tools, CFD-ACE+ and FLUENT, are discussed in the following sections. The first model
is based in an empirical determination of the dynamic contact angle determined from the
experiments in the preceding Chapter. The second model is based in the theoretical
framework and Hoffman function discussed in Section 2.5.5.
 55

4.2.1. CFD-ACE+
The dynamical dependence of the contact angle on the local tangential speed of the
liquid-gas interface is important and is often specified based on empirical observations.
The dependence of the contact angle on the local gas flow velocity enables the simulation
of dynamic contact angle (DCA) effects, which are manifested by a difference in the
contact angle for the receding and advancing contact lines. Generally, the contact angle
varies continuously as the motion of the contact line changes from a receding to an
advancing pattern. This variation can be approximated by a contact angle function. The
variation in the contact angle with the velocity can be appreciable, and the required
functional dependence is obtained herein by curve fitting the experimental data. In this
section, two different methods are proposed to derive the velocity dependent contact
angle function; these are then implemented in the free surface module of CFD-ACE+.
The module provides options in the solver to activate the input of empirical dynamic
contact angle relations. The general form of such relations is:

  
 d  f  s , v t , x, s v t ,  lc ,  uc  (4.15)

 v t , if v t  F  0 
  
s vt  

 (4.16)
 v t , if v t  F  0

where
d: dynamic contact angle
s: static contact angle
x : location on boundary surface

v t : tangential velocity at solid boundary

v t : magnitude of v t

 
s v t : sign magnitude of v t

lc: lower cut off angle

uc: upper cut off angle
F : gradient of volume fraction.
 56
Several function forms are available in the module, such as single ramp, double ramp,
quadratic, cubic, hyperbolic tangent function, etc. The selection of the function and
associated parameters has to be made such that the contact angle function remains
bounded and smooth since the velocity might depart from the anticipated range especially
in the convergence process. Therefore, a lower contact angle cut-off of zero and upper
cut-off of 180 degrees is prescribed [94].

4.2.1.1. Method 1: derivation from droplet emergence experiment

The measurements of DCA and corresponding experimental conditions were described
in Chapter 3. The derivation of the velocity dependent contact angle function is based on
the water injection velocity (Vw) at 0.04 m/s and air flow velocity (Va) at 10m/s (as listed
as the Case 2 condition in Table 3.1). To obtain the velocity dependent contact angle
function during the period of one emergence cycle, i.e. from the onset of droplet
emergence into the air stream to the time just before it detaches from the water injection
pore. A snapshot of the droplet illustrates the associated parameters in Figure 4.1.

Figure 4.1. Image of water droplet subjected to the air flow stream. Points A and B are
the receding and advancing points in a 2-D plane of view, whereas r and a designate
the receding and advancing contact angle, respectively.

When a droplet emerges into the air stream, it will pin temporarily at receding point A,
however, the advancing point B will keep moving toward the channel outlet along the
surface until the droplet reaches a certain size, after which, it detaches and starts moving
along the surface. Details of the DCA evolution through an emergence cycle were
presented in Section 3.2.4. The relative motion of advancing point B is characterized as a
moving contact line. It should be noted that the experimental analysis are based on 2D
images, whereas the contact line of a droplet moves in a plane orthogonal to the image.
 57
Between the advancing leading edge and the receding trailing edge, it is assumed that
the contact line is connected in the same manner as the contact line velocity. Coordinates
of point B relative to point A as well as the contact angle distribution in one emerging
cycle in the time-series images are analyzed and presented in Figure 4.2.
To specify the dependence of the contact angle on the contact line velocity of a liquid-
gas interface along the local surface, the position data were fitted to a second order
polynomial function; the first derivative of this function was used to obtain the contact
line velocity function.

Figure 4.2. Position of advancing point (XB), velocity function and DCA distribution
for droplet emergence cycle Case 2.

The contact line advancing velocity, Vc, varies from about 0.001 m/s to 0.008 m/s
(where the droplet starts to detach). Plots of the positions and advancing velocity versus
time and the corresponding advancing contact angle are shown in Figure 4.3. The final
form of the velocity dependent contact angle function is deduced as a single ramp
function:

 D  25552 Vc  25 (4.17)

 58
It should be noted that the data and connection are somewhat ambiguous during the
initial emergence when the contact line is influenced by two surfaces: the pore surface
and the channel surface which are normal to each other. This occurs when the droplet just
enters the air channel and has to overcome the air stream pressure in the channel. The
droplet has to move from one quasi equilibrium position to the other and must overcome
the energy barrier between the interface of droplet and surface. An external disturbance
or a sudden change of equilibrium condition can affect the contact angle on a rough
surface [95].

Figure 4.3. Velocity dependent contact angle function from droplet emergence
experiment.

4.2.1.2. Method 2: derivation from capillary rise experiment

A cross sectional dimension of 250 µm square of PDMS microchannel was fed with
water from one side, with the other side open to atmosphere. The PDMS chip was placed
on the same experimental platform as described in Chapter 3. Water was pumped from
one side to the other, and the evolution of the shape of the water meniscus and of the
advancing contact angle was observed over a range of pumping velocities. Similarly, the
receding contact angle was observed and measured by sucking the water back into the
syringe pump as illustrated schematically in Figure 4.4.
 59

Figure 4.4. Schematic diagram of capillary rise experiment.

The flow images were analyzed using contact angle measurement. For each run of this
experiment, there is only one contact angle and only one velocity. This experiment
provides a more accurate measurement for advancing and receding angles. The pumping
flow rate provides the advancing data with positive velocity, and the withdrawing flow
rate provides the receding data with negative velocity. The resulting velocity dependent
contact angle function is presented in Figure 4.5. The data obtained with the first method
are also plotted for comparison.

Figure 4.5. Velocity dependent contact angle function.

The two different curve fitted functions, single and double ramp, were implemented as
the input parameters for the VOF free surface module of CFD-ACE+ using the dynamic
contact angle option. The final forms of the velocity dependent contact angle functions
for the single and double ramp cases, are:
 60
single ramp :  D  23263  Vc  76 (4.18)

  26250  Vc  67 for Vc  0
double ramp :  D (4.19)
 D  18744  Vc  67 for Vc  0

The double ramp function allows representation of the change in slope for the advancing
and receding contact angle. Since prior to detachment, the velocity is close to zero and
therefore the contact angle corresponds closely to the intercept angle (67) in the curve,
the specification of the intercept angle might affect the emergence phenomena slightly in
the simulations.

4.2.2. ANSYS FLUENT

Using user defined function (UDF) in FLUENT, the Hoffman function (cf. Eq.(2.11))
is implemented in two different steps in subroutines written in C code [51]. The first step
is to determine the unit normal for the phase boundary. The second step is to provide the
main program with a local contact angle by determining the local capillary number and
evaluating the Hoffman function. Details of the source code are listed in Appendix C.
In the step of determining the unit normal, the gradients and other phase information is
released from memory after the calculation of the phase composition. A dummy source
term (DEFINE_SOURCE) is reserved to exchange the information for evaluating the
VOF equations together with the associated data. The solver will call this dummy source
term code every time the VOF equations are solved. An additional function
(DEFINE_ADJUST) is also required to setup the system variables and allocate a memory
buffer to store the related information. When the dummy source term code is called, it
collects the phase gradients (C_VOF_G) for the boundary, and stores them in a user
defined function (C_UDMI). The code then returns a source term of zero to the
equations. This is a key step to enable the VOF gradient data to be accessible by the main
contact angle code specifically when the boundary conditions are imposed.
The second step involves the evaluation of local contact angle through the UDF
(DEFINE_PROFILE) to determine the Hoffman function which has the form of:

 D  f Hoff [Ca  f Hoff

1
( e )] (4.20)
 61
The inverse Hoffman function needs to be evaluated based on the static contact angle
through a zero finding function. The main code obtains the phase gradients (C_VOF_G)
from the UDF which was saved in the first step. The phase gradients are then brought
into the code and normalized. The contact line velocity can then be determined by taking
the dot product of flow field velocity (C_U, C_V, C_W) for the mixture with the
normalized unit vector along the direction of the interface normal vector (NV_DOT). The
capillary number is then determined by this velocity normal and the Hoffman function is
evaluated. Finally, the dynamic contact angle is returned to the main program at each
location and the next iteration starts.

4.3. Simulation Domain and Mesh, Boundary and

Initial Conditions

4.3.1. CFD-ACE+
A schematic of the three dimensional computational domain and the corresponding
mesh is shown in Figure 4.6.

Figure 4.6. Three-dimensional domain and mesh for the numerical simulations of
droplet emergence using CFD-ACE+.

The domain consists of a 3-D geometry with square air and water channels of width 250
µm and 50 µm, respectively. These dimensions are identical to those in the experimental
 62
setup up described in Chapter 3, except for the channel length used in numerical
simulation is 1,000 µm. While this is shorter than in the experiments to reduce
computational costs, the simulations still cover the flow region of interest. Uniform mesh
interval 12.5 µm (32,128 cells) and 6.25 µm (257,024 cells) were used. The uniform
velocity profiles are specified at the channel air inlet and water inlet. A convective
outflow condition is used at the channel outlet. Air flow enters the channel from the right-
hand side at a velocity U = 10 m/s, pressure of 101.3 kPa and temperature of 298K. The
water is injected into the channel with a velocity of V = 0.04 m/s through the pore located
on the midline of bottom surface (cf. Case 2 condition in Table 3.1). A no-slip boundary
condition is applied on all of the surfaces. Surface tension and contact angle are specified
on the wall as the boundary condition. The simulations are initialized with a uniform air
velocity field and no liquid water present in the channel. For the base case study, a static
contact angle (SCA), s = 110, is specified on all the surfaces; this coincides with the
experimental conditions in Chapter 3 which a microchannel using PDMS
(Polydimethylsiloxane) material. For the dynamic contact angle (DCA) simulation, a
velocity dependent contact angle function (discussed in Section 4.2.1) is specified. To
minimize the computation time without interfering with stability, the size of the next time
step Δt was computed before every new time step based on fixed Courant number (Co =
0.2) and the variable local velocity in the cells. Thus the largest velocity in the domain
determines the maximum time step. The resulting time step was of the order of 10-8 s ~
10-7 s. For the geometric reconstruction of the interfaces between water droplets and air, a
PLIC (piecewise linear interface calculation) method was used in conjunction with
Gauss’s theorem.

4.3.2. ANSYS FLUENT

Two different simulation conditions are discussed herein. The first addresses model
validation for the dynamic contact angle formulation using droplet impact on a horizontal
surface. The second one describes the condition for droplet emergence simulation used
for dynamic contact angle modeling.
 63
4.3.2.1. Droplet impact
The numerical code for dynamic contact angle is first validated using a well-
established two phase flow study before application to the droplet emergence simulation.
The droplet impact on a horizontal surface selected as a “benchmark” case is well
discussed in the literature and exhibits spreading and recoiling phenomenon that
correspond to advancing and receding contact lines [43]. Though a three dimensional
model is the best to predict the entire droplet impact process, it was shown that a two
dimensional axisymmetric simulation can represent the process faithfully [43].

Figure 4.7. Illustration of numerical grids used for the droplet impact computations.
The region of an adaptive refinement is presented.
 64
Thus, a two dimensional axisymmetric simulation is used in this study where the
changes in azimuthal direction are not considered for a droplet normal impact on a
horizontal wax surface. The computational domain consists a two dimensional plan of 12
× 12 mm large square in x-y plane as shown in Figure 4.7 where a single free falling
droplet along the x-direction impacts a plane in the y-direction. The largest uniform mesh
interval is 40 µm away from the boundary. Further mesh refinement using adaptive grid
is applied to the near wall boundary as well as the center line (cf. black region in Figure
4.7) where the major liquid droplet movement occurs. The mesh refinement is critical in
order to capture the particular details of lamella motion near the wall boundary. The
finest cell near the wall is 2.5 × 2.5 µm. The total number of cells used in this case is
597,447.
Water is used as the test liquid and encounters a normal impact on a horizontal waxed
surface. The static contact angle for water on the wax is in the range of 95 - 105,
therefore 105 is used for the SCA case study. Flow conditions and fluid properties are
summarized in Table 4.1. The water droplet size D is 2.7 mm and the impact Webber
number We is 90, equivalent to an impact velocity of 1.57 m/s along the x direction in
Figure 4.7. The droplet is initialized at the instant of impact using the FLUENT UDF
code for initialization which is listed in Appendix D.

Table 4.1. Properties of liquid and flow conditions.

Liquid  µ ρ D V
(N/m) (Ns/m2) (kg/m3) (mm) We (m/s) Re s

Water 0.0728 0.001 998.3 2.7 90 1.57 4223 105

Where  is surface tension, µ is dynamic viscosity, ρ is density, D is droplet diameter

before impact, V is impact velocity, Re is impact Reynolds number and s is the static
contact angle.
A variable time step is applied by controlling the Courant number, Co = 0.25, as a result
the time step is of the order of 10-6 s for SCA and 10-8 s DCA modeling. The introduction
 65
of the DCA model makes the numerical convergence process much slower and
necessitates variable time steps that are typically two orders of magnitude lower.

4.3.2.2. Droplet emergence

Figure 4.8 presents the three dimensional computational domain and corresponding
mesh for droplet emergence modeling. The width of the square air and water channels is
250 µm and 50 µm, respectively, and is described in Section 4.3.1.

Air
water

Figure 4.8. Computation domain and mesh illustration for the droplet emergence
study, presenting the region of an adaptive refinement.

These dimensions are identical to the experimental ones described in Chapter 3, except
for the channel length of 1,500 µm used in numerical simulation. Flow conditions of air
and water inlet velocities are 10 m/s and 0.04 m/s, respectively. These conditions are the
same in the case 2 condition in Table 3.1. A uniform mesh interval of 10 µm is used
 66
throughout the domain, except for the near-pore subdomain which is further refined.
To enhance the accuracy of the simulations and fully resolve the dynamic droplet
emergence and detachment, adaptive gridding was used to refine the mesh in the near
pore region and corresponding boundaries (cf. zoom-in region in Figure 4.8). The size of
the adaptive refinement region is adjusted to cover one emerging droplet in order to
capture the contact line movement. The finest cell near the wall is 2.5 × 2.5 × 2.5 µm.
The total number of cells used in this case is 896,700. A factor of 0.5 for under-relaxation
is applied for the pressure based solver control. The adaptive refinement and relaxation
factor play a key role for the droplet emergence simulation. Variable time step is applied
by setting the Courant number, Co = 0.25, as a result the time step is of the order of 10-7
~ 10-6 s for SCA modeling and 10-9 ~ 10-8 s for DCA modeling.
 67

4.4. Results and Discussions

4.4.1. CFD-ACE+

4.4.1.1. Static contact angle (SCA): base case

The base case using SCA is first presented at the instant in time when the droplet
detaches (Figure 4.9).

Top view

Side view

Figure 4.9. (a) Top view. (b) Side view of the instant at droplet detachment using SCA.

Examination of the animation sequence shows that each individual droplet emerging
from the water pore remains pinned around the pore due to surface tension forces until it
grows into a critical size before detaches from the pore. The effect of air velocity on the
changing shape of the emerging droplet and the time evolution of the emerging process
can be observed. The periodic emerging frequency is around 142.6 Hz (T = 7.0137 ms).
Images at different instants are analyzed to evaluate the advancing and receding contact
angles and further compared with the experimental results presented in Chapter 3. The
droplet shape at detachment differs significantly from the shape observed in the
 68
experiments. The long water trailing observed in the experiment does not present in the
simulations in which instead the droplet maintains a fuller height. Deviation of dynamic
contact angle evolution in an emergence cycle using normalized time scale (τ = t / T) is
presented in Figure 4.10.

Figure 4.10. Comparison of contact angle evolution of VOF simulation (SCA, mesh
12.5um) and experiments.

While the evolution of the advancing angle is reproduced reasonably by the simulation,
the static angle formulation results in there significant discrepancy between the
simulation and experimental receding angles. Additionally, a smoother evolution rather
than the stick and jump motion that is observed in the experiment. This highlights the
effect of surface roughness or inhomogeneities which are not taken into account in this
study. It’s worth noting that the droplet emergence is around 10.8 times faster than that of
the experimental results (142.6 Hz versus 13.2 Hz), and the simulation droplets are
smaller than observed in the experiment immediately after detachment.
 69
4.4.1.2. Dynamic contact angle (DCA): Method 1, Eq. (4.17)
DCA simulation results using Eq. (4.17) derived from the droplet emergence experiment
are presented in Figure 4.11.

(a)

(b)
Figure 4.11. (a) Side view of the instant at droplet detachment using DCA Eq.
(4.17). (b) Comparison of contact angle evolution of VOF simulation (DCA,
method 1, mesh 12.5 m) and experiment.
 70
The simulated droplet exhibits water trailing but its shape isn’t otherwise faithful to the
experiment. The contact angle evolution shows a monotonic decrease before crossing
over the experimental curve and over predicting the contact angle. The emergence
frequency is around 137 Hz. It is shown that using method 1 for deriving the velocity
dependent contact angle function Eq. (4.17) is clearly deficient for describing the contact
line motion. In this method, the data is deduced from point A and point B (cf. Figure 4.1
and Figure 4.2) and thus characterize the whole droplet behavior based on only two
points. Whereas, in realty, there is a different contact angle vs. velocity relation at each
point of the curve defined by the intersection of the droplet surface and the bottom wall
on which it moves. The contact angle at point A and B is not a function of their velocity
only. There are different contact angles along the footprint and they all influence each
other by deforming the droplet surface. Using a two dimensional image and axisymmetric
assumption is not suitable to delineate the velocity dependent contact angle function since
the effect of through plan motion is not well described. In addition, there is a window of
about 17 ms (cf. Figure 3.8 Period I and II) during which data points include a surface
tension dominant and transition region when the droplet is growing but held in place by
the strong surface tension; this phenomenon skews the data set and the velocity. Finally,
Eq. (4.17) doesn’t include a negative velocity dependent contact angle relation.

4.4.1.3. Dynamic contact angle (DCA): Method 2, Eq. (4.19)

Figure 4.12 presents the results using DCA Eq. (4.19) with a uniform mesh of 12.5 m.
The droplet emergence frequency is around 114 Hz, which is slightly lower than the SCA
results. An improvement of contact angle evolution is observed. Though the receding
angle is still not well captured, the evolution of the advancing angle is in better agreement
with the experiments. The surface tension dominates period (cf. Figure 3.8) ends around 
= 0.15, however no transition period is observed. Although some of the discrepancies can
be attributed to the dynamic contact angle model, detailed examination of these results
also indicated that a finer mesh is required by the VOF technique to resolve the droplet
dynamics, especially close to the solid surface.
 71

(a)

Figure 4.12. (a) Side view of the instant at droplet detachment using DCA Eq. (4.19)
and [3]. (b) Comparison of contact angle evolution of VOF simulation (DCA, method
2, mesh 12.5 m) and experiment.

Figure 4.13 presents the results using a finer 6.25 m mesh. This mesh resolution yields a
significant improvement in the droplet emergence frequency, but the value of about 70
Hz, is still five times higher than the experimental results. Although further mesh
refinement might be warranted, the discrepancy is primarily attributed to the effect of
surface inhomogeneities and roughness effects that are not accounted for in the
simulations.
 72

Figure 4.13. (a) Side view of the instant at droplet detachment using DCA Eq. (4.19)
and [3]. (b) Comparison of contact angle evolution of VOF simulation (DCA, method
2, mesh 6.25 m) and experiment.
 73

4.4.2. ANSYS FLUENT

4.4.2.1. Droplet impact on a horizontal surface

Base case study using SCA under the flow condition as listed in Table 4.1 is discussed
herein for the benchmark problem of a free falling droplet impacting a horizontal surface.
The results are compared against the experimental and numerical results presented by
Šikalo and Ganić [96] and are shown in Figure 4.14.

Figure 4.14. Comparison of time sequence of water droplet impact onto wax surface
(We = 90), experiment (left) and numerical (right) (reproduced from [96] with
permission of Experimental Thermal and Fluid Science).

As pointed out by Šikalo and Ganić [96], the liquid film (lamella) thickness is less than
6µm at the beginning of the impact observed from the experiment. This observation
emphasizes the importance of using grid refinement at the wall boundary. It was noted
that due to the depth field in the experimental images, the experimental images present
some blurriness especially during the later spreading process which exhibits wrinkling on
the edge due to circumferential instability. Consequently, only the sharpest image
boundary should be considered when comparing with the numerical results. The image at
t = 11.0 ms shows an uprising central jet in the axisymmetric simulation which azimuthal
variations present in the latter phase of the experiments.
The time series images during the spreading and recoil phases using SCA model are
shown in Figure 4.15 and Figure 4.16, respectively.
 74

t = 0.06 ms

t = 0.15 ms

t = 0.30 ms

t = 0.41 ms

t = 0.52 ms

t = 0.60 ms

t = 0.85 ms

t = 1.01 ms

t = 1.25 ms

t = 1.96 ms

t = 2.50 ms

t = 3.03 ms

t = 3.86 ms

Figure 4.15. Time series images during the spreading phase for SCA modeling of
droplet impact on wax surface.
 75

t = 4.65 ms

t = 5.02 ms

t = 5.72 ms

t = 6.01 ms

t = 7.03 ms

t = 7.80 ms

t = 8.01 ms

t = 9.50 ms

t = 10.0 ms

t = 11.0 ms (not in scale)

Figure 4.16. Time series images during the recoiling phase for SCA modeling of
droplet impact on wax surface.

Additionally, the time series images during the spreading and recoiling phases using the
DCA model are shown in Figure 4.17 and Figure 4.18, respectively.
 76

t = 0.07 ms

t = 0.14 ms

t = 0.30 ms

t = 0.43 ms

t = 0.54 ms

t = 0.60 ms

t = 0.86 ms

t = 1.01 ms

t = 1.25 ms

t = 1.96 ms

t = 2.54 ms

t = 3.06 ms

t = 3.80 ms

Figure 4.17. Time series images during the spreading phase for DCA modeling of
droplet impact on wax surface.
 77

t = 4.64 ms

t = 5.08 ms

t = 5.70 ms

t = 6.03 ms

t = 7.07 ms

t = 7.88 ms

t = 8.10 ms

t = 9.60 ms

t = 10.0 ms

t = 11.0 ms

t = 12.0 ms

t = 13.0 ms (not in scale)

Figure 4.18. Time series images during the recoiling phase for DCA modeling of
droplet impact on wax surface.
 78
In each case, an air micro-bubble is shown near the impact point at the origin of the
coordinate system due to the deformation of the droplet surface at the instant of contact.
The entrapment of air bubble is numerically and experimentally observed over a wide
range of impact conditions by Mehdi-Nejad et al. [97] and Thoroddsen et al. [98],
respectively. The non-dimensional spreading diameter (d/D) and apex height (y/D) of a
droplet, shown schematically in Figure 4.19, are quantified for further analysis as a
functions of dimensionless time (tV/D) from impact.

Figure 4.19. Schematics of spreading diameter and apex height of drop impacts.

The numerical are with the work of Šikalo et. al [81] for the spreading diameter are
shown in Figure 4.20 for different contact angle approach.

Figure 4.20. Numerical simulation of the temporal evolution of the spread

diameter in comparison with the results of Sikalo et. al [81].
 79
The present results agree with the experimental data very well in the early spreading
phase when the inertial force is dominant in the early phase and the effect of wettability is
relatively minor. During the later advancing phase of the spreading while close to the
maximum spreading distance and the drop begins to recoil, the results of using DCA
algorithm provide better evidence the ability to track the wettability effect together with
the drop deformation. The maximum spreading diameter using DCA approach is about
1.2% larger than the numerical results in Šikalo et. al [81]. The simulated and
experimental apex heights are also compared in Figure 4.21.

Figure 4.21. Numerical simulation of the temporal evolution of the apex height in
comparison with the results of Šikalo et. al [99].

The dimensionless apex height decreases from 1, the highest value at the instant of
impact, to around 0.05 in the spreading phase. There is no significant difference in the
spreading phase between DCA and SCA simulations, and the results coincide with the
experimental results of Šikalo et. al [99]. This is expected since the inertial force is
dominant in the spreading phase until the droplet reaches the maximum spreading
distance; by the time the recoil phase takes over as a viscous effect starts to play a major
role. During the recoil phase, the DCA approach follows the experimental data faithfully.
 80
4.4.2.2. Droplet emergence in a fuel cell channel
The base case for droplet emergence study is using SCA approach. Figure 4.22 shows the
time evolution of the droplet in one emergence cycle.

Figure 4.22. One cycle of droplet emergence time resolved images of numerical results
using static contact angle (SCA, s = 110) approach in FLUENT. Air inlet velocity,
Va = 10m/s; water inlet velocity, Vw = 0.04 m/s.

The edges inside the channel delineate the near-wall grid refinement region as shown in
Figure 4.8. The droplet is essentially hemispherical before 5 ms, and then tilts slightly
toward the channel outlet after 7 ms. The droplet detaches from the water pore at around
25.4 ms which corresponds to a shedding frequency of 39 Hz.
 81

Figure 4.23. One cycle of droplet emergence time resolved images of numerical results
using dynamic contact angle (DCA) approach in FLUENT. Air inlet velocity, Va = 10
m/s; water inlet velocity, Vw = 0.04 m/s.

Similarly, the time evolution using dynamic contact angle (DCA) simulation is presented
in Figure 4.23. In comparison with the results in SCA, the droplet interface in DCA is
less spherical with a longer chord length and a shorter height. Droplets emerge faster
compared to the SCA case, and detach from the water pore at around 22.8 ms which is
equivalent to a shedding frequency of 44 Hz. The emergence frequency is smaller in SCA
simulations thus allowing the droplet to grow slightly higher in comparison with DCA
simulations.
 82
To further investigate the treatment of the dynamic contact line, the advancing and
receding contact angles as well as the droplet chord length are presented in Figure 4.24
for both SCA and DCA using FLUENT.

Figure 4.24. Comparison of experimental result and the evolution of dynamic contact
angle using FLUENT.

The results using DCA approach provide a significant improvement in dynamic contact
angle as well as the chord length compare to the SCA results. Considerable deviation in
receding contact angle between SCA and experimental results is observed. Using DCA
approach, the results of receding angle as well as the chord length agree with the
experiments well overall except in the early emergence phase ( < 0.2). The higher
advancing and receding contact angles in early emergence phase cause the droplet to
form a nearly spherical shape.
 83
Furthermore, comparing the results of DCA simulations using CFD-ACE+ (cf.
Figure 4.13) and FLUENT, both numerical approaches track the dynamics of the droplet
reasonably well in a dimensionless time scale. Although the emergence frequency (70 Hz
in CFD-ACE+ vs. 44 Hz in FLUENT) are still higher than the value in experiment (13.2
Hz in Figure 3.5). Note that a uniform grid size 6.25 µm is used in CFD-ACE+
simulation, whereas the finest adaptive grid size for FLUENT is 2.5 µm. The CFD-ACE+
simulations are consistent with the experimental observations in the early emergence
phase (i.e. period I: surface tension force dominant in Figure 3.8). However, they deviate
significantly for the receding angle in the later part of the emergence phase (i.e. period
III: drag force dominant). The DCA approach in CFD-ACE+ using Eq. (4.19) derived
from the one-dimensional capillary rise experiment provides better capability to capture
the wettability in the surface tension dominated regime but is not satisfactory in the drag
force dominated regime where the effect of flow field is important in the vicinity of the
moving contact line. On the other hand, the DCA modeling using the Hoffman function
in FLUENT presents better results in both regime though deviations still exist in the
surface tension force dominated period which might be due to surface roughness and
inhomogeneities physically present in the experiment and not accounted for in the
numerical model.

4.5. Summary
The dynamics of single water droplet emerging from a pore in the presence of a cross
flow of air was investigated numerically in a modeled PEMFC cathode gas channel. The
importance of accounting for the pore connectivity during the initial emergence phase
was highlighted in the prior work of Zhu et al. [56]. In the present work, we explicitly
resolve the pore and also examine the critical role of dynamic contact angle modelling.
Different approaches of dynamic contact angle implementations were examined, and the
results were compared with the results using static contact angle model as well as the
experimental measurements. The emphasis here is placed upon the implementation of
two empirical velocity dependent contact angle functions as well as the theoretically
based Hoffmann equation in two different numerical software CFD-ACE+ and FLUENT.
There are two main conclusions of this study. First, grid refinement is critical in order to
 84
capture the salient motion of a moving droplet. With the extensive use of nonuniform,
adaptive refinement in the region of droplet movement the dynamic evolution of contact
angles are well modeled in one emergence cycle without sacrificing the computational
cost. Second, a velocity dependent contact angle function derived either from droplet
emergence experiments or from one-dimensional capillary rise experiments doesn’t
provide satisfactory tracking of the contact line motion of an emerging droplet in a single
channel because droplets encounter different flow regimes and emergence phase.
Significantly improved simulations are obtained by introducing the Hoffmann equation
which takes into account more fundamental aspects of the local contact line velocity as a
function of the flow field. The dynamic contact angle model was validated against
numerically and experimentally documented droplet impacts on a horizontal surface.
Spreading and recoiling parameters were compared and presented for validation. This
model was then used to investigate a flow condition from our own experiments in a
model fuel cell microchannel and the results show that the Hoffmann equation improves
the accuracy of predicting the advancing and receding contact angle of an emerging
droplet inside a channel.
 85

Chapter 5

5 Pressure Signature and Diagnostic Tool

Experimental data and signal analysis of pressure drop measurement from a two-
phase flow fixture are presented in this chapter. The pressure drop setup consists of a
single-channel fixture, flow delivery for air and water, data acquisition and a high speed
digital camera. The flow fixture was designed to allow replaceable GDL and bipolar
plates for testing of different channel geometries and materials. High accuracy dynamic
pressure sensors were installed at different locations of the fixture. Experiments were
conducted for flow rates of air and liquid water similar to those in an operating fuel cell.
Power spectrum analysis was performed on the pressure time series data and the results
were used to construct a flow regime map. When comparing with footage from the high
speed digital camera, it was found that the presence of water at the outlet was the major
cause of a low-frequency periodic oscillation of pressure signal. The present study
concludes with several suggestions for future work, including verification of data
repeatability, multiple water injection, modification of design of flow plate outlet to avoid
water accumulation, and real time monitoring using power spectrum of pressure signal.

5.1. Introduction
Understanding the behaviour of liquid water in PEMFC gas channels is a key to
effective water management. From an on-board diagnostic aspect, the pressure signature
of a unit cell shows significant advantages in providing immediate and real time linkage
between the air/water two-phase flow behaviours [3]. An increase of pressure drop would
 86
indicate water accumulation or water film build-up in the channel, while an abrupt
decrease would imply water removal at the outlet [18]. However, using pressure drop as a
diagnostic tool for the water behaviour, particularly in water removal, still has
deficiencies. For example, high frequency oscillations in pressure drop signal, which
arise from unstable two-phase flow transport in the channel, are usually observed. Such
oscillations might not indicate any water removal/build-up at the outlet yet they are
difficult to distinguish and separate from regular increase/decrease patterns of pressure
drop profile. Moreover, considering the requirement of relatively low pressure drop
(usually < 10 mbar) across fuel cell channels, signal noises from disturbances and
measurement errors may contribute considerably to the results as well as reducing the
accuracy and even changing the pattern of the pressure drop profile. Real-time
visualizations are used here to assist the interpretation of the pressure drop profiles, and
to identify the potential of using the pressure drop as a diagnostic tool for two-phase flow
in PEMFCs. The specific objectives are to perform pressure signal measurement in the
two-phase flow fixture with liquid water fluxes representative of real operating
conditions (i.e. current density) and identify and delineate the water management regimes
in the channel and GDLs.

5.2. Method and Apparatus

Figure 5.1 shows the experimental setup which comprises a two-phase flow fixture
(supplied by Ballard Power System), water pump (Harvard Apparatus, PHD2000), mass
flow controller of air (Brooks, MFC 5850i), differential pressure sensor (Honeywell,
FP2000), gauge pressure senor (Noshok, 100 series) and data acquisition system
(National Instruments, LabVIEW 8.6 and DAQ). The combination of small hydraulic
diameter of the flow channel with dry (air only) and two-phase flow regimes, the pressure
drops span a large range of values making precise measurements challenging. The mass
flow controller provides fine adjustment of air flow rate from 0 to 1,000 sccm in a 1%
full scale precision. The measurement range of the differential pressure sensor is up to 5
psi, with a 0.1% accuracy and 3 KHz frequency response.
 87

Figure 5.1. Experimental setup in pressure drop measurement.

The user control interface was developed using LabVIEW 8.6 to provide real-time
equipment control, monitoring and data logging. The two-phase flow fixture can be
configured with several options including a transparent plastic flow channel or a non-
transparent carbon flow plate (CFP), and an untreated carbon paper which is clamped
together using an air pressurized device. The clamping pressure (Pb) is maintained at
around 75 psi to provide reliable sealing without leakage. The chamber of the clamping
device provides 29 pressure ports (denoted as P01 to P29) and 28 water injection ports
(W01 to W28) which are 15 mm apart, respectively. Figure 5.2 illustrates the
corresponding locations of pressure and water ports.

Figure 5.2. Illustration of the pressure measurement and water injection ports.
 88
A Toray TGP-H-060 non-woven fibrous GDL material was used. A series of 1 mm
diameter holes were punched into the GDL at the corresponding locations of pressure and
water ports. In this arrangement, the water droplet is thus allowed to enter the flow
channel through the 1 mm GDL holes. A design sample of the flow channel plate is
presented in Figure 5.3 and is either made from transparent plastic or non-transparent
carbon.

Figure 5.3. Flow channel plates (Ballard Power System). Unit: mm

The schematics and dimensions of several flow channels are shown in Figure 5.4.
Dimensions are measured under a microscope to ensure accuracy in prescribing the
hydraulic diameter (Dh) required for calculating the theoretical pressure drop across the
channel. Note that the drawing of CFP-0D is not shown here due to confidentiality.

Figure 5.4. Schematics and dimensions of the flow channels.

The pressure drop (P) in this study is measured between P01 and P28. First of all, the
system is examined by measuring the pressure drop of a dry channel (i.e. only gas phase)
under different air flow rates. A sample of the real time signal of the flow rate and the
 89
pressure drop (P0128) across port P01 and P28 in flow plate Plastic-0A is shown in
Figure 5.5. The corresponding Reynolds number range from 110 to 1100 where a laminar
flow regime is maintained.

Figure 5.5. Real time signal of flow rate versus pressure drop in flow plate Plastic-0A.

Four different flow plates are investigated and the comparison of Darcy friction factor
between the experimental and theoretical value using Poiseuille equation is shown in
Figure 5.6. The Darcy friction factor, f, determined from the pressure drop measurement
is expressed as:
2P0128 Dh
f  (5.1)
 AVA 2 L0128
where L0128 is the channel length between pressure port P01 and P28. The curves are in
good agreement. This verification provides a calibration curve and a reference base for
wet channel measurements.
 90

Figure 5.6. Pressure drop calibration curve for various flow plates.

The matrix of flow conditions (each element represents the case number) shown in Table
5.1 comprises the combinations of air and water flow rates derived from typical fuel cell
operating conditions. Two plastic plates (0A and 02) and two carbon flow plates (010A
and 0D) are examined using this matrix.

Table 5.1. Test conditions matrix

 91
In each measurement, the pressure drop signal is initially recorded for three minutes in
the dry state to ensure the pressure drop P0128 falls within the calibration curve range
shown in Figure 5.6. Single water injection from water port W01 as well as dual injection
from W01 and W15 are both tested for certain flow plates. To maintain the amount of
water injected into the channel constant, the flow rate is divided by two for the dual
injection case. After each measurement and before starting a new test, the flow channel is
air purged for at least 30 minutes to ensure both channel and GDL are in a dry state.

5.3. Results and Discussions

5.3.1. Characterization of mean pressure drop

A sample of the real-time pressure drop signal of Plastic-0A flow plate using single
injection and under the flow condition of Case 11 is shown in Figure 5.7. While water is
introduced, a spike appears in the signal due to the water accumulation in the chamber
before being pushed into the channel; some oscillations are also evident. A time span of
more than 20 minutes is recorded and averaged to represent the mean pressure drop for a
given test condition. The mean pressure drop in this case is 10,506 Pa compare to a value
of 9,731 Pa in a dry state. That is, an increase of 775 Pa due to water coverage in the
channel is observed. The increase of pressure drop due to water accumulation is a key
challenge of water management and may lead to flow misdistribution, and reduction of
the fuel cell performance and durability.
The signal exhibits various oscillation patterns in different stage which indicate
dynamic characteristics. Though there seems to be a characteristic frequency present in
some intervals, a dominated value could not be identified. The dynamic features underly
the complexity of the flow and the challenge in achieving an analytic model of the
process.
 92

Figure 5.7. Pressure drop signal of Plastic-0A using single injection under the Case 11
test condition.

The ratio of two-phase flow to the single gas phase pressure drop can be used as a simple
indicator of the liquid buildup in fuel cell flow channel. This ratio, also known as the
normalized pressure drop or the two-phase friction multiplier [24], is defined by Eq.
(5.2):
P2
g 2  (5.2)
Pg

where P2 and Pg is the pressure drop with a two-phase flow and a single-phase gas
flow in the channel, respectively. A higher ratio represents a higher water buildup in the
channel. The mean pressure drop of Plastic-0A flow channel using single injection for
different flow conditions is listed in Table 5.2. The additional database of pressure drop
in different test cases is listed in Appendix E. It is shown that while varying the air flow
rate (Case 11 to 14), the difference of pressure drop (P2 - Pg) decreases; though the
friction multiplier (g2) increases. The highest friction multiplier which appears in Case
24, implies a higher tendency of slug flow is formed. Whereas, in highest air flow rate
 93
(Qa = 407 sccm), the friction multiplier tends to unity and implies a strong air flow
inducing a mist two-phase flow regime in the channel.

Table 5.2. Pressure drop measurement of Plastic-0A using single water injection.

Plastic-0A single injection, P0128

Dry Wet Wet-Dry Wet/Dry
Case Qw Qa ReA
(Pg) (P2 ) (P2  Pg ( g 2 )
#11 407 891 9,731 10,506 775 1.080
#12 306 670 7,276 7,924 648 1.089
0.075
#13 204 447 4,814 5,359 545 1.113
#14 61 134 1,488 1,794 306 1.206
#21 407 891 9,699 10,322 623 1.064
#22 306 670 7,238 8,083 845 1.117
0.056
#23 204 447 4,818 5,279 461 1.096
#24 61 134 1,431 1,752 321 1.224
#31 407 891 9,804 10,625 821 1.084
#32 306 670 7,307 7,858 551 1.075
0.037
#33 204 447 4,792 5,429 637 1.133
#34 61 134 1,422 1,679 257 1.181
#41 407 891 9,793 9,983 190 1.019
#42 306 670 7,325 7,555 230 1.031
0.011
#43 204 447 4,850 5,248 398 1.082
#44 61 134 1,435 1,698 263 1.183

The difference of pressure drop between dry and wet channel is expected to be controlled
at below 1 kPa (10 mbar) for a given flow channel in order to minimize the excess power
loss1. In this study, as shown in Table 12.2 in Appendix E, the difference of pressure drop
of Plate-02 is worth noticing since the majority of its value exceeds 1 kPa except in the
lowest water injection rate (i.e. Case 41 to Case 44). The flow channel geometry of
Plastic-02 might not be optimized when large amount of water exist in the gas channel at
high current density conditions. However, CFP-010A and CFP-0D flow channels present
a relatively low pressure drop deviation and are all below 1 kPa criterion. As discussed in
Figure 2.5, a secondary channel design on top of a triangular gas channel can provide
passive water removal due to capillary forces. This type of channel design provides an

1
Private communication, Ballard Power Systems
 94
effective way to let water path through the secondary channel in the upper portion
while maintaining a suitable air passage in the lower portion of the channel. Plots of the
difference of pressure drop in test conditions are shown in Figure 5.8.

(a)

(b)
Figure 5.8 The effects of air flow (ReA) and flow channel geometry on the
difference between wet and dry pressure drop under single water injection. (a)
Plastic-0A and Plastic-02, (b) CFP-010A and CFP-0D.
 95

(a)

(b)
Figure 5.9. The effects of air flow (ReA) and flow channel geometry on the friction
multiplier under single water injection. (a) Plastic-0A and Plastic-02, (b) CFP-010A
and CFP-0D.

To better describe the effect of air flow on the water coverage in the flow channels,
comparison of the friction multiplier using single water injection are plotted against the
 96
Reynolds number of air flow as shown in Figure 5.9. The friction multiplier decreases
as ReA increases. It generally varies from nearly 1.0 to around 1.7. However, at lowest
ReA, Plastic-02 channel exhibits a higher ratio in some test conditions which indicate a
higher tendency for slug flow formation. Moreover, carbon flow plate CFP-0D reveals a
lower friction multiplier ratio in most cases.

5.3.2. Dynamic characteristics

An oscillation of the pressure drop signal is observed in Plastic-0A flow channel plate
using single water injection under Case 34 flow conditions. A sample of real time
pressure drop signal is presented in Figure 5.10.

Figure 5.10. Real-time pressure drop signal of Plastic-0A channel using single water
injection under Case 34 flow conditions.

The oscillation frequency is identified as 2.8 Hz using DFT (Discrete Fourier Transform).
This frequency was found to be correlated to the water accumulation and drainage
process at the opening of channel outlet. The accumulated water block the opening for a
certain period of time which causes the pressure to increase until it levels off when water
has been purged out of the outlet. Some snap shot images shown in Figure 5.11 highlights
 97
this process. It is worth to noting that a proper design of the end of the channel linking
to the outlet as well as the shape and surface property of the outlet are critical in order to
enhance water expulsion.

Figure 5.11. Snap shot images of water accumulation and drainage process.

Though the frequency is related to the water expelling process, it does not present in all
the flow plates in this study. Therefore, we focus our attempt to devise a flow regime
diagnostic relation to plate CFP-010A due to its lower pressure drop difference as
discussed earlier.

Figure 5.12. Typical real time pressure drop signal of CFP-010A channel using single
water injection under Case 14 flow condition.

Figure 5.12 shows the pressure drop signal using single water injection in CFP-010A
under Case 14 flow condition. Under this condition, the air flow rate is at the lowest
value whereas the water flow rate is highest. The analysis reveals lower frequency
oscillations which implies slug regime in the channel. This is consistent with the fact that
slug flow is characterized by a slower water build up process. The overall picture of the
 98
frequency analysis is shown in Figure 5.13. The frequencies range mostly below 1 Hz
but do not present a clear trend for identifying the flow regimes.

Figure 5.13. Dominant frequency in two-phase flow regime of CFP-010A.

At lower water injection rate e.g. Qw = 0.011 and Qw = 0.037, increasing the air flow
rate reduces the frequency from about the order of 10 Hz to 1 Hz indicating a flow
regime changes. On the other hand, the trend is reversed while using higher water
injection rate e.g. Qw = 0.056 and Qw = 0.075 since the frequency increase with air flow
rate. To further investigate the evolution of flow regime and how the strength of the
pressure drop signal is distributed in the frequency domain respective to the strengths of
other background signals associated with the flow regime, the power spectral density
(PSD) of the pressure drop signal is deduced by squaring its amplitude. A plot of power
spectrum of CFP-010A using different flow conditions to illustrate the extreme
conditions is presented in Figure 5.14. That is, at highest water injection, Case 11 and
Case 14 represent highest air flow versus lowest air flow, respectively. At lowest water
injection, Case 41 and Case 44 represent the effect of highest and lowest of air flow.
 99

Figure 5.14. Power spectrum of Case 11, 14, 41 and 44 of CFP 010A.

The amplitude of the spectrum reveals the change between flow regimes. For slug slow
regime, the amplitude of the signal surpasses that for the other flow regime even though
the frequency is low. This implies a strong pulse-like signal relative to the background
which might due to the water accumulation and purging process.

Figure 5.15. Flow regime identification using DFT (Discrete Fourier Transform)
power of the pressure drop signal.
 100
Figure 5.15 synthesizes the relationship between the spectra power amplitude and the
flow regimes under the effects of water and air flow rates. The amplitude corresponds to
the highest value in each test condition. This successfully identifies the boundaries
between different flow regimes under various test conditions. These boundaries might
change slightly if test conditions varied further. This is a distinct, order of magnitude
difference between flow regimes. The power spectrum analysis of the pressure drop is
proposed as a useful and practical approach to identify and characterize dominant flow
patterns and regime where in operating fuel cells where flow visualization is not possible
due to the opaque flow channel. This approach also has the potential to be extended to
real time monitoring of flow regimes in the channels and stack.

5.4. Summary
The ability to monitor or diagnose flow regime change and the onset of flooding in
operating fuel cells would open the door for active control of flow conditions to achieve
more effective water management strategies. The experiments and analysis presented in
the Chapter show that pressure signal analysis cab be used to provide real time
monitoring of air/water two-phase flow behaviors. The experimental work was conducted
to characterize the two-phase flow in a single flow channel using an ex-situ two-phase
flow fixture. Several flow channels design identical to those used in the commercial fuel
cell stacks were tested under different operating conditions. Database of pressure drop as
well as the characterization of flow regimes and the development of flow diagnostics tool
were explored. A flow channel with secondary channel design located on the top was
found to provide a better passage to lift the water due to the effect of capillary force. The
accumulation and purging process of water at the outlet is the primary cause for the low
frequency periodic oscillation of pressure drop due to water. Finally, using power
spectrum analysis on the pressure drop signal provides a suitable approach for
constructing the flow regime.
 101

Chapter 6

6 Conclusion and future works

6.1. Conclusion
An experimental and numerical investigation of the two-phase flow in a fuel cell gas
channel was presented in this thesis. First, observations of the dynamics of a single
droplet emerging from a pore into a model PEMFC flow channel were made using
microfluidics technology. The development of the transparent microchip platform
allowed quantitative flow visualization of droplet dynamic under a range of air flows and
water injection rates. Different flow patterns were analyzed and a flow map was
constructed, and the time evolution of the contact line during the droplet emergence and
detachment cycle was determined. These results set the foundation for validation of the
numerical model. Numerical simulations using the VOF method were then undertaken
with to resolving the time-dependent behaviour of a single droplet throughout and
emergence cycle. The focus of the modeling was treatments of the contact angle at the
three-phase interface on a hydrophobic surface. Both empirical and theoretical velocity
dependent contact angle functions were implemented in CFD software and their
performance examined through a set of simulations representing experimental conditions.
The last part of the thesis was focused on identifying relationships between pressure drop
signals and flow regimes as a basis for developing design and diagnostic tools for water
management. A two-phase flow fixture consisting of gas flow channel identical to those
used in a commercial fuel cell stack was used to further study the two-phase flow
behaviour via ex-situ measurements and observations. This two-phase flow fixture
 102
provides allowed examination of the effect of channel geometry and flow conditions
on the pressure signature. The results were analyzed and the potential for using pressure
drop signals as a diagnostics tool was discussed.
The main contributions of this work are:
1. Single droplet dynamics: The detachment and subsequent dynamic evolution of
single droplet in microchannel was analyzed experimentally. Three different flow
patterns were identified and a flow regime map was obtained. The time evolutions
of the advancing and receding contact angles in the droplet regime were also
determined.
2. Modeling of dynamic contact angle: VOF simulations with a dynamic contact
angle treatment were performed for the first time. A velocity dependent contact
angle function deduced from experiments as well as theoretical contact line
relation (Hoffmann function) were implemented into VOF simulations. The
simulations highlight the critical importance of a dynamic contact line treatment
and show that using the Hoffmann function yields more physical simulations.
3. Flow diagnostic tool: Pressure signature were acquired and analyzed using a
custom designed two-phase flow apparatus with flow field plates identical to
those in commercial stacks. The results provide evidence of the importance of
channel design with secondary channels. It is shown that the flow regimes can be
characterized using power spectrum density of the normalized pressure drop
signal. This provides a good basis for ex-situ diagnostics.

6.2. Future works

Future work includes further investigation on the experimental and numerical works
for droplet emergence as well as the flow diagnostics tool.

Experimental work on droplet emergence

In present study, only smooth surface channels were examined. To further study the
effect of surface properties, implementation of surface roughness is essential. Some
preliminary experiments were performed, but were not analyzed as systematic
experimental tests were not completed due to the challenges of microfabrication process.
 103
Experiments with multiple water injection would also be valuable to approach fuel
cell operating conditions when simultaneous emergence and interaction of multiple
droplets occurs.

Numerical work on droplet emergence

To complement the experimental work, modelling developments should be undertaken
to extend the dynamic angle model (Hoffmann function) to rough surfaces representative
of GDLs. Simulations with further grid refinement should also be undertaken with high
performance computing resources.

Flow diagnostics tool

The pressure signal-flow regime correlation identified in this study opens the
possibility for practical flow diagnostics. Several avenues could be pursues: (1)
repeatability verification of pressure drop; (2) the design of the flow plate outlet to
minimize the water accumulation and increase drainage capability; (3) real time
monitoring using power spectrum of pressure drop; (4) validation and further analysis
with the VOF; (5) multiple water injections.
 104

7 Bibliography
[1] N. Shao, A. Gavriilidis, P. Angeli, Flow regimes for adiabatic gas-liquid flow in
microchannels, Chem. Eng. Sci. 64 (2009) 2749–2761.

[2] T. Berning, N. Djilali, A 3D, multiphase, multicomponent model of the cathode

and anode of a PEM fuel cell, J. Electrochem. Soc. 150 (2003) A1589–A1598.

[3] T.A. Trabold, Minichannels in Polymer Electrolyte Membrane Fuel Cells, Heat
Transf. Eng. 26 (2005) 3–12.

[4] A. Bazylak, J. Heinrich, N. Djilali, D. Sinton, Liquid water transport between

graphite paper and a solid surface, J. Power Sources. 185 (2008) 1147–1153.

[5] K. Tüber, D. Pócza, C. Hebling, K. Tuber, D. Pocza, Visualization of water

buildup in the cathode of a transparent PEM fuel cell, J. Power Sources. 124 (2003)
403–414.

[6] X.G. Yang, F.Y. Zhang, A.L. Lubawy, C.Y. Wang, Visualization of liquid water
transport in a PEFC, Electrochem. Solid State Lett. 7 (2004) A408–A411.

[7] K.S. Chen, M.A. Hickner, D.R. Noble, Simplified models for predicting the onset
of liquid water droplet instability at the gas diffusion layer/gas flow channel
interface, Int. J. Energy Res. 29 (2005) 1113–1132.

[8] F.Y. Zhang, X.G. Yang, C.Y. Wang, Liquid water removal from a polymer
electrolyte fuel cell, J. Electrochem. Soc. 153 (2006) A225–A232.

[9] E.C. Kumbur, K. V Sharp, M.M. Mench, Liquid droplet behavior and instability in
a polymer electrolyte fuel cell flow channel, J. Power Sources. 161 (2006) 333–
345.

[10] A. Theodorakakos, T. Ous, A. Gavaises, J.M. Nouri, N. Nikolopoulos, H.

Yanagihara, Dynamics of water droplets detached from porous surfaces of
relevance to PEM fuel cells, J. Colloid Interface Sci. 300 (2006) 673–687.

[11] 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.

[12] A. Gunther, S.A. Khan, M. Thalmann, F. Trachsel, K.F. Jensen, Transport and
reaction in microscale segmented gas-liquid flow, Lab Chip. 4 (2004) 278–286.
 105
[13] S. Waelchli, P.R. von Rohr, Two-phase flow characteristics in gas-liquid
microreactors, Int. J. Multiph. Flow. 32 (2006) 791–806.

[14] N. Kim, E.T. Evans, D.S. Park, S.A. Soper, M.C. Murphy, D.E. Nikitopoulos,
Gas-liquid two-phase flows in rectangular polymer micro-channels, Exp. Fluids.
51 (2011) 373–393.

[15] I.S. Hussaini, C.-Y. Wang, Visualization and quantification of cathode channel
flooding in PEM fuel cells, J. Power Sources. 187 (2009) 444–451.

[16] C.H. Hidrovo, F.-M. Wang, J.E. Steinbrenner, E.S. Lee, S.S. Vigneron, C.-H.
Cheng, et al., Water Slug Detachment in Two-Phase Hydrophobic Microchannel
Flows, in: Proc. 3rd Int. Conf. Microchannels Minichannels, ASME, 2005: pp.
ICMM2005–75261.

[17] Z.J. Lu, C. Rath, G.S. Zhang, S.G. Kandlikar, Water management studies in PEM
fuel cells, part IV: Effects of channel surface wettability, geometry and orientation
on the two-phase flow in parallel gas channels, Int. J. Hydrogen Energy. 36 (2011)
9864–9875.

[18] Z. Lu, S.G. Kandlikar, C. Rath, M. Grimm, W. Domigan, A.D. White, 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, Int. J.
Hydrogen Energy. 34 (2009) 3445–3456.

[19] C.E. Colosqui, M.J. Cheah, I.G. Kevrekidis, J.B. Benziger, Droplet and slug
formation in polymer electrolyte membrane fuel cell flow channels: The role of
interfacial forces, J. Power Sources. 196 (2011) 10057–10068.

[20] G. Minor, N. Djilali, D. Sinton, P. Oshkai, Flow within a water droplet subjected
to an air stream in a hydrophobic microchannel, Fluid Dyn. Res. 41 (2009).

[21] N. Yousfi-Steiner, P. Mocoteguy, D. Candusso, D. Hissel, A. Hernandez, A.

Aslanides, A review on PEM voltage degradation associated with water
management: Impacts, influent factors and characterization, J. Power Sources. 183
(2008) 260–274.

[22] P. Rodatz, F. Büchi, C. Onder, L. Guzzella, Operational aspects of a large PEFC

stack under practical conditions, J. Power Sources. 128 (2004) 208–217.

[23] W. He, G. Lin, T. Van Nguyen, Diagnostic tool to detect electrode flooding in
proton-exchange-membrane fuel cells, Aiche J. 49 (2003) 3221–3228.

[24] R.W. Lockhart, R.C. Martinelli, Proposed correlation of data for isothermal two-
phase two-component flow in pipes, Chem. Eng. Prog. 45 (1949) 39–48.
 106
[25] Y. Wang, S. Basu, C.Y. Wang, Modeling two-phase flow in PEM fuel cell
channels, J. Power Sources. 179 (2008) 603–617.

[26] N. Akhtar, A. Qureshi, J. Scholta, C. Hartnig, M. Messerschmidt, W. Lehnert,

Investigation of water droplet kinetics and optimization of channel geometry for
PEM fuel cell cathodes, Int. J. Hydrogen Energy. 34 (2009) 3104–3111.

[27] J. Scholta, G. Escher, W. Zhang, L. Kuppers, L. Jorissen, W. Lehnert,

Investigation on the influence of channel geometries on PEMFC performance, J.
Power Sources. 155 (2006) 66–71.

[28] S. Shimpalee, J.W. Van Zee, Numerical studies on rib & channel dimension of
flow-field on PEMFC performance, Int. J. Hydrogen Energy. 32 (2007) 842–856.

[29] X. Zhu, Q. Liao, P.C. Sui, N. Djilali, Numerical investigation of water droplet
dynamics in a low-temperature fuel cell microchannel: Effect of channel geometry,
J. Power Sources. 195 (2010) 801–812.

[30] T. Metz, N. Paust, C. Mfiller, R. Zengerle, P. Koltay, Passive water removal in

fuel cells by capillary droplet actuation, Sensors and Actuators a-Physical. 143
(2008) 49–57.

[31] J.P. Owejan, T.A. Trabold, D.L. Jacobson, M. Arif, S.G. Kandlikar, Effects of
flow field and diffusion layer properties on water accumulation in a PEM fuel cell,
Int. J. Hydrogen Energy. 32 (2007) 4489–4502.

[32] M. Cheah, I.G. Kevrekidis, J.B. Benziger, Water slug formation and motion in gas
flow channels: the effects of geometry, surface wettability, and gravity, Langmuir.
29 (2013) 9918–9934.

[33] N. Djilali, P.C. Sui, Transport phenomena in fuel cells: from microscale to
macroscale, Int. J. Comut. Fluid Dyn. 22 (2008) 115–133.

[34] C. Fang, C. Hidrovo, F. Wang, J. Eaton, K. Goodson, 3-D numerical simulation of

contact angle hysteresis for microscale two phase flow, Int. J. Multiph. Flow. 34
(2008) 690–705.

[35] K. Tüber, D. Pócza, C. Hebling, Visualization of water buildup in the cathode of a

transparent PEM fuel cell, J. Power Sources. 124 (2003) 403–414.

[36] X. Zhu, P.C. Sui, N. Djilali, Three-dimensional numerical simulations of water

droplet dynamics in a PEMFC gas channel, J. Power Sources. 181 (2008) 101–115.

[37] A.D. Le, B. Zhou, A general model of proton exchange membrane fuel cell, J.
Power Sources. 182 (2008) 197–222.
 107
[38] Y. Tabe, Y. Lee, T. Chikahisa, M. Kozakai, Numerical simulation of liquid
water and gas flow in a channel and a simplified gas diffusion layer model of
polymer electrolyte membrane fuel cells using the lattice Boltzmann method, J.
Power Sources. 193 (2009) 24–31.

[39] 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. 195 (2010) 7278–7288.

[40] S.G. Mukherjee, A; Kandlikar, A Numerical Analysis of Growing Water Droplet

inside an Air Supply Channel of a PEM Fuel Cell, in: 2006.

[41] J. Choi, G. Son, Numerical study of droplet dynamics in a PEMFC gas channel
with multiple pores, J. Mech. Sci. Technol. 23 (2009) 1765–1772.

[42] T.D. Blake, The physics of moving wetting lines., J. Colloid Interface Sci. 299
(2006) 1–13.

[43] S. Sikalo, H.-D.D. Wilhelm, I. V. Roisman, S. Jakirlic, C. Tropea, S. Šikalo, et al.,

Dynamic contact angle of spreading droplets: Experiments and simulations, Phys.
Fluids. 17 (2005) 062103.

[44] F.Y. Zhang, X.G. Yang, C.Y. Wang, Liquid Water Removal from a Polymer
Electrolyte Fuel Cell, J. Electrochem. Soc. 153 (2006) A225.

[45] C.H. Schillberg, S.G. Kandlikar, A Review of Models for Water Droplet
Detachment From the Gas Diffusion Layer-Gas Flow Channel Interface in
PEMFCs, ASME Conf. Proc. 2007 (2007) 299–310.

[46] E.B. Dussan V., R.T.-P. Chow, On the ability of drops or bubbles to stick to non-
horizontal surfaces of solids, J. Fluid Mech. 137 (1983) 1–29.

[47] E.B. Dussan V., On the ability of drops or bubbles to stick to non-horizontal
surfaces of solids. Part 2. Small drops or bubbles having contact angles of arbitrary
size, J. Fluid Mech. 151 (1985) 1–20.

[48] E.B. Dussan V., On the ability of drops to stick to surfaces of solids. Part 3. The
influences of the motion of the surrounding fluid on dislodging drops, J. Fluid
Mech. 174 (1987) 381–397.

[49] P. Dimitrakopoulos, J.J.L. Higdon, Displacement of fluid droplets from solid

surfaces in low-Reynolds-number shear flows, J. Fluid Mech. 377 (1997) 351–378.

[50] P. Dimitrakopoulos, J.J.L. Higdon, On the displacement of three-dimensional fluid

droplets from solid surfaces in low-Reynolds-number shear flows, J. Fluid Mech.
377 (1998) 189–222.
 108
[51] C. Miller, Liquid Water Dynamics in a Model Polymer Electrolyte Fuel Cell
Flow Channel, Dept. Mech. Eng. MASc (2010).

[52] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free
boundaries, J. Comput. Phys. 39 (1981) 201–225.

[53] P. Quan, B. Zhou, A. Sobiesiak, Z.S. Liu, Water behavior in serpentine micro-
channel for proton exchange membrane fuel cell cathode, J. Power Sources. 152
(2005) 131–145.

[54] Y.H. Cai, J. Hu, H.P. Ma, B.L. Yi, H.M. Zhang, Effects of
hydrophilic/hydrophobic properties on the water behavior in the micro-channels of
a proton exchange membrane fuel cell, J. Power Sources. 161 (2006) 843–848.

[55] X. Zhu, P.C. Sui, N. Djilali, Dynamic behaviour of liquid water emerging from a
GDL pore into a PEMFC gas flow channel, J. Power Sources. 172 (2007) 287–295.

[56] X. Zhu, P.C. Sui, N. Djilali, Numerical simulation of emergence of a water droplet
from a pore into a microchannel gas stream, Microfluid. Nanofluidics. 4 (2008)
543–555.

[57] K.S. Chen, Modeling Water-Droplet Detachment From GDL/Channel Interface in

PEM Fuel Cells, ASME Conf. Proc. 2008 (2008) 797–803.

[58] Y. Ding, H.T. Bi, D.P. Wilkinson, Three dimensional numerical simulation of gas-
liquid two-phase flow patterns in a polymer-electrolyte membrane fuel cells gas
flow channel, J. Power Sources. 196 (2011) 6284–6292.

[59] B. Sundén, Transport phenomena in fuel cells, WIT Press, Boston, 2005.

[60] J.A. Sethian, Level set methods and fast marching methods : evolving interfaces in
computational geometry, fluid mechanics, computer vision, and materials sciences,
Cambridge Univ. Press, Cambridge [u.a.], 2006.

[61] M. Sussman, P. Smereka, S. Osher, A Level Set Approach for Computing

Solutions to Incompressible Two-Phase Flow, J. Comput. Phys. 114 (1994) 146–
159.

[62] M. Sussman, A.S. Almgren, J.B. Bell, P. Colella, L.H. Howell, M.L. Welcome, An
Adaptive Level Set Approach for Incompressible Two-Phase Flows, J. Comput.
Phys. 148 (1999) 81–124.

[63] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev.
Fluid Mech. 30 (1998) 329–364.
 109
[64] X.W. Shan, H.D. Chen, Lattice Boltzmann model for simulating flows with
multiple phases and components, Phys. Rev. E. 47 (1993) 1815–1819.

[65] K. Fei, C.W. Hong, All-angle removal of CO2 bubbles from the anode
microchannels of a micro fuel cell by lattice-Boltzmann simulation, Microfluid.
Nanofluidics. 3 (2007) 77–88.

[66] J. Park, X. Li, Multi-phase micro-scale flow simulation in the electrodes of a PEM
fuel cell by lattice Boltzmann method, J. Power Sources. 178 (2008) 248–257.

[67] L. Hao, P. Cheng, Lattice Boltzmann simulations of liquid droplet dynamic

behavior on a hydrophobic surface of a gas flow channel, J. Power Sources. 190
(2009) 435–446.

[68] A.K. Gunstensen, D.H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model
of immiscible fluids, Phys. Rev. A. 43 (1991) 4320–4327.

[69] M.R. Swift, E. Orlandini, W.R. Osborn, J.M. Yeomans, Lattice Boltzmann
simulations of liquid-gas and binary fluid systems, Phys. Rev. E. 54 (1996) 5041–
5052.

[70] X. He, G.D. Doolen, Thermodynamic Foundations of Kinetic Theory and Lattice
Boltzmann Models for Multiphase Flows, J. Stat. Phys. 107 (2002) 309–328.

[71] Y.Y. Yan, Y.Q. Zu, A lattice Boltzmann method for incompressible two-phase
flows on partial wetting surface with large density ratio, J. Comput. Phys. 227
(2007) 763–775.

[72] T. Inamuro, T. Ogata, S. Tajima, N. Konishi, A lattice Boltzmann method for

incompressible two-phase flows with large density differences, J. Comput. Phys.
198 (2004) 628–644.

[73] A.J. Briant, A.J. Wagner, J.M. Yeomans, Lattice Boltzmann simulations of contact
line motion. I. Liquid-gas systems, Phys. Rev. E. 69 (2004) 31602.

[74] X.D. Niu, T. Munekata, S.A. Hyodo, K. Suga, An investigation of water-gas

transport processes in the gas-diffusion-layer of a PEM fuel cell by a multiphase
multiple-relaxation-time lattice Boltzmann model, J. Power Sources. 172 (2007)
542–552.

[75] I. V Roisman, L. Opfer, C. Tropea, M. Raessi, J. Mostaghimi, S. Chandra, Drop

impact onto a dry surface: Role of the dynamic contact angle, Colloids Surfaces A
Physicochem. Eng. Asp. 322 (2008) 183–191.

[76] Y.D. Shikhmurzaev, The moving contact line on a smooth solid surface, Int. J.
Multiph. Flow. 19 (1993) 589–610.
 110
[77] T.D. Blake, Y.D. Shikhmurzaev, Dynamic wetting by liquids of different
viscosity, J. Colloid Interface Sci. 253 (2002) 196–202.

[78] C. Huh, S.G. Mason, The steady movement of a liquid meniscus in a capillary tube,
J. Fluid Mech. Digit. Arch. 81 (1977) 401–419.

[79] P.A. Durbin, Considerations on the moving contact-line singularity, with

application to frictional drag on a slender drop, J. Fluid Mech. Digit. Arch. 197
(1988) 157–169.

[80] C. Baiocci, V. V Pukhnachev, Problems with one-sided constraints for Navier-

Stokes equations and the dynamic contact angle, J. Appl. Mech. Tech. Phys. 31
(1990) 185–197.

[81] S. Šikalo, H.-D.D. Wilhelm, I. V. Roisman, S. Jakirlić, C. Tropea, S. Sikalo, et al.,

Dynamic contact angle of spreading droplets: Experiments and simulations, Phys.
Fluids. 17 (2005) 062103.

[82] R.G. Cox, Inertial and viscous effects on dynamic contact angles, J. Fluid Mech.
357 (1998) 249–278.

[83] S. Mukherjee, J. Abraham, Investigations of drop impact on dry walls with a

lattice-Boltzmann model, J. Colloid Interface Sci. 312 (2007) 341–354.

[84] M. Sussman, M. Ohta, A stable and efficient method for treating surface tension in
incompressible two-phase flow, Siam J. Sci. Comput. 31 (2009) 2447–2471.

[85] T.C. Wu, N. Djilali, Qualitative and quantitative analysis of moving droplets in a
modelled PEMFC cathode gas microchannel, in: 20th Int. Symp. Transp. Phenom.,
Victoria, BC, Canada, 2009.

[86] E. Kimball, T. Whitaker, Y.G. Kevrekidis, J.B. Benziger, Drops , Slugs , and
Flooding in Polymer Electrolyte Membrane Fuel Cells, Aiche J. 54 (2008) 1313–
1332.

[87] A. Bazylak, D. Sinton, N. Djilali, Dynamic water transport and droplet emergence
in PEMFC gas diffusion layers, J. Power Sources. 176 (2008) 240–246.

[88] N.T. Nguyen, S.H. Chan, Micromachined polymer electrolyte membrane and
direct methanol fuel cells - a review, J. Micromechanics Microengineering. 16
(2006) R1–R12.

[89] S. Litster, D. Sinton, N. Djilali, Ex situ visualization of liquid water transport in

PEM fuel cell gas diffusion layers, J. Power Sources. 154 (2006) 95–105.
 111
[90] M.K. Akbar, D.A. Plummer, S.M. Ghiaasiaan, On gas-liquid two-phase flow
regimes in microchannels, Int. J. Multiph. Flow. 29 (2003) 855–865.

[91] ANSYS FLUENT V14.0 Theory Guide, 15317 (2011).

[92] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface
tension, J. Comput. Phys. 100 (1992) 335–354.

[93] D.L. Youngs, Time-dependent multi-material flow with large fluid distortion,
Numer. Methods Fluid Dyn. 24 (1982) 273–285.

[94] CFD-ACE+ v2010.0 Modules Manual, Part 2, (2010).

[95] B. He, N.A. Patankar, J. Lee, Multiple equilibrium droplet shapes and design
criterion for rough hydrophobic surfaces, Langmuir. 19 (2003) 4999–5003.

[96] S. Sikalo, E.N. Ganic, Phenomena of droplet-surface interactions, Exp. Therm.

Fluid Sci. 31 (2006) 97–110.

[97] V. Mehdi-Nejad, J. Mostaghimi, S. Chandra, Air bubble entrapment under an

impacting droplet, Phys. Fluids. 15 (2003) 173.

[98] S.T. Thoroddsen, T.G. Etoh, K. Takehara, N. Ootsuka, Y. Hatsuki, The air bubble
entrapped under a drop impacting on a solid surface, J. Fluid Mech. 545 (2005)
203.

[99] S. Sikalo, M. Marengo, C. Tropea, E.N. Ganic, Analysis of impact of droplets on

horizontal surfaces, Exp. Therm. 25 (2002) 503–510.

[100] A. Mata, A.J. Fleischman, S. Roy, Fabrication of multi-layer SU-8 microstructures,

J. Micromechanics Microengineering. 16 (2006) 276–284.

[101] Y.N. Xia, G.M. Whitesides, Soft lithography, Annu. Rev. Mater. Sci. 28 (1998)
153–184.
 112

8 Appendix A – Microfluidic chip

fabrication
Chips incorporating the microchannels were fabricated on a 3″ silicon wafer by a
sequential multi-layer SU-8 (MicroChem Corp.) fabrication process [100]. Two different
high resolution laser photo-plotted masks were prepared for the water and air delivery
networks. A 25 µm layer of SU-8 25 was first spin coated, baked, exposed and developed
to form the water delivery path. This process was repeated for a second 125 µm SU-8 100
layer to form the air delivery path. The finished multi-layer SU-8 structure served as the
molding master for soft lithography techniques [101]. Two different mixture ratios of
PDMS, 15:1 (with excess vinyl groups) and 5:1 (with excess Si-H groups), were cast
separately on the multi-layer molding master and thermally cured on a hotplate at 80 ºC
for 40 minutes. The cured PDMS polymer was peeled off from the master mold, punched
for fluidic access for water and air, diced, ethanol cleaned and finally air dried. The
molding master was treated with a few drops of Repel-Silane ES (General Electric
PlusOne) liquid for 5 minutes before each use for casting to prevent adhesion of the
PDMS on the SU-8. Two symmetric elastomers with different mixture ratios of PDMS
were carefully aligned to each other by predefined alignment marks under a stereo
microscope providing ±5 µm precision. The aligned chip was thermal bonded on a
hotplate at 80 ºC for 2 hours or overnight to increase the strong interfacial bonding
between the two pieces. The secure bonding between the two pieces is important to
prevent any leakage in the air stream at pressures up to 3 atmospheres. After bonding, a
complete chip consists of a 50 µm square water pore and a 250 µm square air channel as
shown in Figure 3.1a.
 113

9 Appendix B – MATLAB source code

Correlation coefficient and droplet frequency
114
115
 116

10 Appendix C – FLUENT UDF source

code for DCA
117
118
119
 120

11 Appendix D – FLUENT UDF source

code for droplet impact
 121

12 Appendix E – Pressure drop

measurement data

Table 12.1. Pressure drop measurement of Plastic-0A using dual water injection.

Plastic-0A dual injection, P0128

Dry Wet Wet-Dry Wet/Dry
Case Qw Qa ReA
(Pg) (P2 ) (P2  Pg ( g 2 )
#11 407 891 9,686 10,491 805 1.083
#12 306 670 7,253 7,947 694 1.096
0.075
#13 204 447 4,818 5,308 490 1.102
#14 61 134 1,447 1,714 267 1.185
#21 407 891 9,707 10,262 555 1.057
#22 306 670 7,268 7,809 541 1.074
0.056
#23 204 447 4,818 5,201 383 1.079
#24 61 134 1,443 1,702 259 1.179
#31 407 891 9,715 10,116 401 1.041
#32 306 670 7,272 7,636 364 1.050
0.037
#33 204 447 4,827 5,134 307 1.064
#34 61 134 1,439 1,641 202 1.140
#41 407 891 9,790 10,084 294 1.030
#42 306 670 7,281 7,419 138 1.019
0.011
#43 204 447 4,824 4,987 163 1.034
#44 61 134 1,416 1,826 410 1.290
 122

Table 12.2. Pressure drop measurement of Plastic-02 using single water injection.

Plastic-02 single injection, P0128

Dry Wet Wet-Dry Wet/Dry
Case Qw Qa ReA
(Pg) (P2 ) (P2  Pg ( g 2 )
#11 407 953 23,517 25,091 1,574 1.067
#12 306 716 17,902 20,534 2,632 1.096
0.075
#13 204 478 11,999 14,265 2,266 1.102
#14 61 143 3,550 6,140 2,590 1.185
#21 407 953 23,450 25,514 2,064 1.057
#22 306 716 17,866 19,817 1,951 1.074
0.056
#23 204 478 11,904 13,621 1,717 1.079
#24 61 143 3,562 5,202 1,640 1.179
#31 407 953 23,513 24,894 1,381 1.041
#32 306 716 17,807 18,968 1,161 1.050
0.037
#33 204 478 11,971 13,326 1,355 1.064
#34 61 143 3,602 4,704 1,102 1.140
#41 407 953 23,687 24,235 548 1.030
#42 306 716 17,995 18,601 606 1.019
0.011
#43 204 478 12,053 12,466 413 1.034
#44 61 143 3,590 4,214 624 1.290
 123

Table 12.3 Pressure drop measurement of CFP-010A using single water injection.

CFP-010A single injection, P0128

Dry Wet Wet-Dry Wet/Dry
Case Qw Qa ReA
(Pg) (P2 ) (P2  Pg ( g 2 )
#11 407 617 2,986 3,380 394 1.132
#12 306 464 2,246 2,493 247 1.110
0.075
#13 204 309 1,485 1,698 213 1.144
#14 61 92 447 568 121 1.271
#21 407 617 2,995 3,231 236 1.079
#22 306 464 2,259 2,475 216 1.096
0.056
#23 204 309 1,463 1,640 177 1.121
#24 61 92 444 598 154 1.348
#31 407 617 2,968 3,152 184 1.062
#32 306 464 2,219 2,375 156 1.070
0.037
#33 204 309 1,470 1,603 133 1.091
#34 61 92 444 526 82 1.185
#41 407 617 2,971 3,068 97 1.033
#42 306 464 2,223 2,296 72 1.033
0.011
#43 204 309 1,473 1,542 68 1.046
#44 61 92 448 483 35 1.078
 124

Table 12.4 Pressure drop measurement of CFP-0D using single water injection.

CFP-0D single injection, P0128

Dry Wet Wet-Dry Wet/Dry
Case Qw Qa ReA
(Pg) (P2 ) (P2  Pg ( g 2 )
#11 407 684 6,543 6,850 307 1.047
#12 306 514 4,911 5,192 281 1.057
0.075
#13 204 343 3,266 3,507 241 1.074
#14 61 103 981 1,186 205 1.209
#21 407 684 6,586 6,862 276 1.042
#22 306 514 4,939 5,174 235 1.048
0.056
#23 204 343 3,275 3,481 206 1.063
#24 61 103 981 1,139 158 1.161
#31 407 684 6,527 6,756 229 1.035
#32 306 514 4,891 5,092 201 1.041
0.037
#33 204 343 3,261 3,428 167 1.051
#34 61 103 981 1,106 125 1.128
#41 407 684 6,528 6,672 145 1.022
#42 306 514 4,884 5,010 126 1.026
0.011
#43 204 343 3,250 3,361 112 1.034
#44 61 103 976 1,072 96 1.098