UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Numerical simulation of interfacial flows in micropores Ahmadlouydarab, Majid 2013

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

Item Metadata

Download

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

Full Text

Numerical Simulation of Interfacial Flows in Micropores by Majid Ahmadlouydarab  B.Sc., Sahand University of Technology, 2001 M.Sc., Sharif University of Technology, 2004  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Chemical and Biological Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) June 2013 c Majid Ahmadlouydarab 2013  Abstract Recent technological developments in microfluidics and fuel cells have given special significance to interfacial dynamics in small pores. Using a diffuseinterface model and a finite-element code, I have simulated three associated problems numerically: gas-liquid flow regimes in micropores; relative permeability for two-phase flow through a model porous medium; and dynamics of sessile drops under the simultaneous action of a wettability gradient and an external flow. For two-phase flows in corrugated microchannels driven by a pressure drop, a number of flow regimes were observed: gas flow, blockage, liquid flow, bubble-slug flow, droplet flow, annular flow and annular-droplet flow. Some of the regimes are known from prior studies in macroscopic pipes, but the others are new and specific to the micropores. Then a map of flow regimes has been constructed in the plane of the liquid saturation and the imposed pressure drop. The transitions among certain flow regimes show significant hysteresis, largely owing to the pinning of the interface at sharp corners in the flow conduit. As an extension of the above study, I computed the relative permeability of a model porous media made of corrugated tubes, using an averaging scheme over a pore-size-distribution of a real porous medium. I discovered that the flow rates vary nonlinearly with the pressure gradient, and that the extended Darcy’s law does not hold in general. In the third project, I found that for each prescribed wetting gradient, there is a narrow range for the cross flow within which a stationary drop can be achieved. The drop motion exhibits strong hysteresis, i.e. sensitivity to initial conditions and forcing history. Two drops merge or separate depending on the competition between wettability and external flow. In general, the wettability gradient favors catch-up and coalescence while the external flow favors separation. These numerical simulations have demonstrated that novel interfacial dynamics can be produced in micropores where capillary forces and contact line dynamics play more important roles than in larger spatial dimensions. The numerical results may serve as guidelines to future experiments and technological development in microfluidics and lab-on-chip devices. ii  Preface This PhD thesis entitled “Numerical Simulation of Interfacial Flows in Micropores” presents the main features of the research that I led and carried out during my PhD study under supervision of Professor J. J. Feng. In this preface, the contributions and collaborations to the papers published or submitted for publication from current thesis are briefly explained. • A version of chapter 3 has been published. M. Ahmadlouydarab, Zhong-Sheng (Simon) Liu, and J. J. Feng (2011), Interfacial flows in corrugated microchannels: flow regimes, transitions and hysteresis. International Journal of Multiphase Flow. 37: 1266-1276. Under supervision of J. J. Feng, and collaboration with Zhong-Sheng (Simon) Liu, I did a comprehensive study of two phase flow regimes in corrugated microchannels and drafted the paper. J. J. Feng put his ideas and helped me to prepare the final version of paper. • A version of chapter 4 has been published. M. Ahmadlouydarab, Zhong-Sheng (Simon) Liu, and J. J. Feng (2012), Relative permeability for two-phase flow through corrugated tubes as model porous media. International Journal of Multiphase Flow. 47: 85-93. Under supervision of J. J. Feng, and collaboration with Zhong-Sheng (Simon) Liu, following first study I studied effects of important parameters on relative permeability of two phase flow in a corrugated tube bundle model and drafted the paper. J. J. Feng helped me to prepare the final version of paper. • A version of chapter 5 has been submitted for publication. M. Ahmadlouydarab and J. J. Feng (2013), Motion and coalescence of sessile drops driven by substrate wetting gradient and external flow. Through a systematic research, I studied the motion and coalescence of a single or a pair of sessile drops on the solid substrate in the presence of wetting gradient and external flow. I conducted this study under supervision of J.J. Feng. I prepared the draft of the paper, but final version of the paper was prepared with help of J. J. Feng. iii  Table of Contents Abstract  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  Preface  Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xv Dedication  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi  1 Introduction . . . . . . . . . . . . . . . . . . . . 1.1 Motivation . . . . . . . . . . . . . . . . . . . 1.2 Objectives and contributions of this research 1.3 Thesis outline . . . . . . . . . . . . . . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  1 1 6 7  2 Methodology Of Research . . . . . . . . 2.1 The diffuse interface model . . . . . . . 2.2 Cahn-Hilliard model parameters . . . . 2.3 Discretization of computational domain 2.4 Numerical algorithm and validation . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  9 9 11 13 13  3 Flow regimes in corrugated microchannels . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Flow regimes . . . . . . . . . . . . . . . . . . . . . 3.3.2 Hysteresis and dependence on initial configuration 3.3.3 Geometric effect: rounded corners . . . . . . . . . 3.3.4 Wettability of the solid . . . . . . . . . . . . . . .  . . . . . . . .  . . . . . . . .  16 16 17 19 19 23 26 28  . . . . .  . . . . .  . . . . .  iv  Table of Contents 3.4 3.5  Comparison with previous experiments . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4 Relative permeability of two phase flow in model porous media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Problem setup and averaging scheme . . . . . . . . . . . . . 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Capillary number effect . . . . . . . . . . . . . . . . . 4.3.2 Saturation effect . . . . . . . . . . . . . . . . . . . . . 4.3.3 Wettability effect . . . . . . . . . . . . . . . . . . . . 4.3.4 Viscosity ratio effect . . . . . . . . . . . . . . . . . . . 4.3.5 Effects of pore geometry and initial configuration . . 4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Motion and coalescence of sessile drops . . . . . . . . . . . . 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Problem setup and methodology . . . . . . . . . . . . . . . . 5.3 Results: the motion of a single sessile drop . . . . . . . . . . 5.3.1 Single drop driven by wettability gradient . . . . . . 5.3.2 Single drop driven by external flow . . . . . . . . . . 5.3.3 Competition between external flow and wetting gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Results: drop-drop coalescence . . . . . . . . . . . . . . . . . 5.4.1 Coalescence driven by wettability gradient . . . . . . 5.4.2 Coalescence driven by external flow . . . . . . . . . . 5.4.3 Competition between external flow and wetting gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . .  29 34 37 37 38 41 42 44 48 51 53 56 58 58 59 61 61 65 67 70 70 71 74 77  6 Conclusion and recommendations . . . . . . . . . . . . . . . 80 6.1 Summary of key findings . . . . . . . . . . . . . . . . . . . . 81 6.1.1 Interfacial flow regimes in corrugated microchannels . 81 6.1.2 Relative permeability of interfacial flows in a model porous media . . . . . . . . . . . . . . . . . . . . . . . 81 6.1.3 Sessile drops motion and coalescence . . . . . . . . . 82 6.2 Significance and limitations . . . . . . . . . . . . . . . . . . . 82 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  85  v  List of Figures 1.1  1.2 2.1  2.2  2.3  3.1  3.2  Schematics of the flow regime maps for microchannels in size of ∼ 1mm oriented (a) horizontally, (b) Vertically [35]. VL and VG are superficial velocities of liquid and gas phases respectively. . . . Water drops formed on the GDL pores of proton exchange fuel cell [56]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A static contact line viewed in (a) the sharp-interface model and (b) the diffuse-interface model. n is the outward normal to the wall. The phase field variable φ represents the two fluid bulks by ±1 and the fluid-fluid interface by φ = 0. Other symbols are defined in the text. Reprinted from Ref. [80] by permission. . . . . . . . . . . . Convergence of the moving contact line to the sharp-interface limit with decreasing Cn, shown for steady interfacial shapes in a Couette flow at three capillary number Ca values. Walls contact angle θs = 90. The m value corresponds to Λ = ld /D = 0.01. x is flow direction and walls are at y = 0 and y = 1. Adapted from Yue et al. [84] with permission. . . . . . . . . . . . . . . . . . . . . . (a) An unstructured triangular mesh generated by GRUMMP with interfacial refinement. The parameters are GF = 3, outer boundary mesh size h2 = 0.5, interior mesh size h3 = h2/2, and interfacial mesh size h1 = h2/64. (b) Magnified view of a portion of the interfacial region. The bold curve indicates the interface φ = 0, which is centered in a belt-like region of refined triangles. Reprinted from Ref. [85] by permission. . . . . . . . . . . . . . . . . . . . . Schematics of the periodic flow geometry. The top of the domain is the axis of symmetry, and periodic boundary conditions are imposed between the left and right sides. Initially the liquid rests in the groove representing a pore chamber. . . . . . . . . . . . . . . Schematic illustration of the condition for contact line pinning. The solid curve is the gas-liquid interface, and the dashed curves indicate its limiting positions before the contact line depins from the corner.  2 5  9  12  14  18  20 vi  List of Figures 3.3  3.4  3.5 3.6  3.7 3.8 3.9  3.10  3.11  3.12  3.13 3.14  Flow regime map for the baseline model. Note that the columns of symbols are spaced for clarity and not according to the exact S values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The gas flow regime achieved for a low saturation (S = 0.38) under a gentle driving force F = 0.25. The meniscus is circular at the start of the simulation, with no flow (t = 0). It deforms moderately in the steady state (t = 43.6). The upper edge of the plots corresponds to the axis of the pore. Time is scaled by µl Rt /σ. . . . . . . . . . Development of complete blockage in the microchannel at F = 0.001 and S = 0.75. . . . . . . . . . . . . . . . . . . . . . . . . Development of the liquid flow regime with F = 0.75 and S = 0.38. Because of the axisymmetry of the geometry, as the liquid moves from the chamber into the pore throat, its projected area on the meridian plane being plotted increases even though its volume is conserved. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Development of the bubble-slug regime at F = 1.75 and S = 0.38. . Development of the annular flow regime at F = 3 and S = 0.38. . Variants of the annular flow regime: (a) droplet flow at F = 3 and S = 0.31; (b) annular-droplet flow at F = 3 and S = 0.75. Both are time-periodic configurations, and the liquid film on the ridge in (a) has stopped moving downstream. . . . . . . . . . . . . . . . . Hysteresis or path dependence of the flow regimes at saturation S = 0.31 (a) and S = 0.71 (b). Other parameters are the same as have produced the flow regime map in Fig. 3.3. . . . . . . . . . . Dependence of flow regimes on the initial configuration. (a) Droplet flow develops from a liquid bridge initially wetting the upstream wall of the pore chamber. (b) Bubble-slug flow develops from a liquid bridge initially in the middle of the pore chamber. (c) Annulardroplet flow develops from a liquid bridge initially attached to the downstream wall of the pore chamber. F = 3, S = 0.38, and the other conditions are the same as in the baseline model. . . . . . . Flow regimes in the baseline model with rounded corners, compared with those taken from the flow map of Fig. 3.3 with sharp corners. S = 0.38 and all other conditions are the same as in Fig. 3.3. . . . Development of the shell flow regime with rounded corners at F = 0.5 and S = 0.38. . . . . . . . . . . . . . . . . . . . . . . . . . . Development of the droplet flow regime with rounded corners at F = 1 and S = 0.38. . . . . . . . . . . . . . . . . . . . . . . . .  21  22 22  23 23 24  24  25  26  28 29 29  vii  List of Figures 3.15 Hysteresis in micropores with rounded corners. At S = 0.38, forward and backward paths are traced by increasing F gradually and then decreasing it, starting from (a) the gas flow regime and (b) the droplet flow regime. The regimes achieved from the baseline model with rounded corners (Fig. 3.12) are shown for comparison. Other parameters are the same as in the baseline model. . . . . . . 3.16 Effect of wall wettability on the interfacial evolution. (a) Annular flow in hydrophilic pore, θs = 70◦ ; (b) gas flow in neutrally wetting pore, θs = 90◦ ; (c) droplet flow in hydrophobic pore, θs = 135◦ . All other parameters are the same: F = 2, S = 0.6, with geometrical parameters α = 0.5, β = 2, γ = 2. . . . . . . . . . . . . . . . . . 3.17 Flow regime map in terms of the superficial velocities: comparison of our numerical prediction with the experiment of [35]. Symbols are our results and solid lines indicate regime boundaries of [35] separating the regimes indicated at the bottom. . . . . . . . . . . 3.18 Gas holdup ratio as a function of the homogeneous gas fraction: comparison between our numerical prediction and the experiment of [72]. The computations are for Re = 0.5 mm and the experimental data are for a circular tube of inner diameter 1.09 mm. . . . . . .  4.1  4.2  4.3  Schematics of the corrugated axisymmetric tube as a model for a pore. The computational domain is half of the meridian plane. The top of the domain is the axis of symmetry, and periodic boundary conditions are imposed between the left and right sides. The geometry is specified by three length ratios: α = Lc /Lt , β = Rc /Rt and γ = Lt /Rt . Initially the liquid forms a disc that blocks the cross-section of the groove representing the pore chamber. . . . . . Pore size distribution (PSD) of a model porous medium. The distribution function ψ(R) is converted from the data of [44] on pore volume distribution. . . . . . . . . . . . . . . . . . . . . . . . . Temporal development of three flow regimes for the baseline parameter values: length ratios α = 1, β = 2, γ = 2; liquid-gas viscosity ratio M = 18, liquid saturation S = 40% and wetting angle θ = 135◦ . Time is scaled by µn Re /σ. (a) Blockage at FR = 0.38, when the imposed pressure fails to dislodge the liquid bridge pinned at the corner of the pore chamber; (b) liquid flow regime at FR = 0.56, with the gas entirely trapped in the pore chamber; (c) drop flow regime at FR = 1.4, with large liquid drops being surrounded and conveyed by a continuous gas phase. . . . .  30  31  33  34  39  40  42  viii  List of Figures 4.4  4.5  4.6  4.7  4.8  4.9  4.10  4.11  4.12  4.13  Dependence of the liquid and gas capillary numbers on the driving force F . All other parameters are at the baseline values given in Fig. 4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Variation of the relative permeability of (a) the wetting (gas) phase and (b) the non-wetting (liquid) phase with the respective capillary number. All other parameters are at the baseline values given for Fig. 4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Relative permeabilities as functions of the liquid (non-wetting) saturation S. F = 1.4 and all other parameters are at their baseline values given for Fig. 4.3. . . . . . . . . . . . . . . . . . . . . . . Comparison of the relative permeabilities Kri (S) with previous studies. (a) Krw of the wetting phase; (b) Krn of the non-wetting phase. Open symbols denote computational results while closed ones experimental results. For our result, the body force F = 1.4 and the other parameters are the same as in Fig. 4.3. . . . . . . . Relative permeabilities as functions of the liquid wetting angle θ on the solid surface. F = 1.4, and the other parameters are at the baseline values of Fig. 4.3. The insets show dominant flow patterns among the pores responsible for the overall Kri . . . . . . . . . . . Effect of the contact angle θ on the relative permeabilities Krw (a) and Krn (b): comparison with lattice-Boltzmann computations in the literature. . . . . . . . . . . . . . . . . . . . . . . . . . . . Relative permeabilities as functions of the liquid-gas viscosity ratio. F = 1.4 and all other parameters are at the baseline values as given for Fig. 4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Flow patterns in a pore of radius Rt = 18.25 µm for increasing M values. Shell flow prevails for M = 0.2 while bubble-slug flow obtains for larger M . All other parameters are at the baseline values as given for Fig. 4.3. . . . . . . . . . . . . . . . . . . . . . Comparison of (a) Krw and (b) Krn as functions of the viscosity ratio M with previous studies. Open symbols denote computations and filled ones experiments. . . . . . . . . . . . . . . . . . . . . Effect of pore geometry and initial interfacial configuration. (a) Pore with sharp corners and the liquid bridge initially at the upstream end of the pore chamber. (b) Pore with rounded corners and the liquid bridge initially at the upstream end of the pore chamber. (c) Pore with rounded corners and the liquid bridge initially at the downstream end of the pore chamber. The parameters are the same for all three cases: effective pore radius Re = 20 µm, α = 1, β = 3, γ = 2, S = 40.5%, FR = 1.4, θ = 135◦ and M = 18. . . . . . . . .  43  45  46  47  49  50  52  53  54  55 ix  List of Figures 5.1  5.2  5.3  5.4  5.5  5.6  5.7  5.8  5.9  Schematic of the initial configuration for 2D planar computations. The drops have an effective diameter De , and the geometry is specified by the length ratios α = L/De , β = H/De , and for two drops, γ = Dc /De . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Experimental conditions of Moumen et al. [54] study used to benchmark our simulation. (a) Contact angle profile along the solid substrate. (b) Profile of the dimensionless wettability gradient G. . . . Comparison of the droplet center velocity, non-dimensionalized into a capillary number, between our 2D computations and the wedgeflow model predictions of [68]. For each Λ value, a slip length ls = 2.5ld = 2.5ΛDe is used in the theoretical model. . . . . . . . . Comparison of the drop center velocity among our 3D computation, experimental measurement and theoretical predictions of the wedge-flow model [54]. Our computation is based on Cn = 2.03 × 10−2 and Λ = 6.76 × 10−3 , and one of the wedge-flow predictions is at the corresponding slip length. A much smaller ls /De = O(10−7 ) is needed to predict the experimental data. . . . Instantaneous velocity of the sessile drop on a substrate with a constant wettability gradient for four values of the viscosity ratio M = 1, 10, 50 and 100. G = 0.144, θm = 92.3◦ . . . . . . . . . . . The average velocity of the drop measured between xc = 2.08 to xc = 4.17 as a function of the dimensionless wettability gradient G. Our computation is compared with the wedge-flow model of [54]. . Dependence of the steady-state drop velocity, in terms of the capillary number Ca, on the driving force F for two contact angles. The insets indicate the degree of drop deformation at the largest F value in each case. . . . . . . . . . . . . . . . . . . . . . . . . . Trajectory of drops starting from three initial shapes with instantaneous contact angles of 45◦ , 90◦ and 135◦ . The drop is initially centered at x0 = 3.33. The constant wetting gradient, with G = 0.245 and θm = 130◦, drives the drop toward the right, while the external flow at F = 6.12 × 10−2 goes to the left. The insets illustrate the initial and final shapes of the drop. . . . . . . . . . . . . . . . . . Final position for a single drop that has come to rest under the opposing actions of F and G. The y-axis shows the displacement of the drop center scaled by the effective drop diameter. Initially the drop is centered at x0 = 3.33. . . . . . . . . . . . . . . . . .  60  62  62  63  64  66  67  69  70  x  List of Figures 5.10 Coalescence of two identical drops on a substrate with a constant  5.11  5.12  5.13  5.14  5.15 5.16  wetting gradient G = 0.244 and maximum contact angle θm = 164◦. The drops are initially centered at x = 3.30 and 4.44 (γ = 1.14), and the two solid curves indicate the trajectories of the front contact line of the trailing drop and the rear contact line of the leading drop. After merging at t = 30, the center of the new drop moves according to the dotted line. . . . . . . . . . . . . . . . . . . . . Coalescence of two identical drops on a substrate with uniform contact angle θ = 105◦ driven by an external flow at F = 0.309 toward the right. The drops are initially centered at x = 1.46 and 2.48 (γ = 1.02), and the curves indicate drop positions in the same way as in Fig. 5.10. Coalescence starts at t = 164. To magnify the view of the curves near the point of coalescence, we have omitted the initial period until t = 150. . . . . . . . . . . . . . . . . . . Comparison of (a) the x component of the velocity and (b) the pressure distribution on two drops at a dimensionless time t = 123. The velocity and pressure values are taken along the contour of φ = 0.9, just outside the interface, and are plotted against the polar angle θ measured from the midpoint between the two contact lines. Velocity and pressure are made dimensionless by σ/µd and σ/De , respectively. . . . . . . . . . . . . . . . . . . . . . . . . . Drop coalescence under co-current F and G: the catch-up time t and catch-up distance d as functions of F for a fixed G = 0.183 and θm = 136◦ . Initially the drops are centered at x0 = 3.30 and 4.44, with initial separation γ = 1.14. The catch-up distance d is defined as the distance travelled by the center of the trailing drop up to the time of coalescence. . . . . . . . . . . . . . . . . . . . Drop coalescence under opposing effects of the external flow (toward the left) and the wetting gradient (toward the right). G = 0.244, θm = 164◦, and F = 0.01. The initial configuration is identical to that of Fig. 5.10, and the curves have the same meanings. Coalescence occurs at t = 46. . . . . . . . . . . . . . . . . . . . . . . . The catch-up time t and catch-up distance d both increase with F , other conditions being the same as in Fig. 5.14. . . . . . . . . . . Drop separation under opposing effects of the external flow (toward the left) and the wetting gradient (toward the right). G = 0.244, θm = 164◦. (a) At F = 0.051, the right drop continues to the right but the left one moves to the left. (b) At F = 0.154, both drops move to the left while separating from each other. . . . . . . . .  71  72  73  75  76 77  78  xi  Nomenclature Symbol A Ad Ae B Ca Cai Cam Can Caw Cn D De fm fw F F FR g G hG hL H HG K Kri Krn Krw ld L m M  Units (SI) m2 m2 m2 N.m−3 − − − − − − m m N.m−2 N.m−1 N − − m2 .s−1 − − − m − m−2 − − − m m m2 .P a−1 .s−1 −  Description Total cross sectional area of porous medium Area of drop Effective cross sectional area of pore Body force Capillary number Capillary number of ith phase External capillary number Capillary number of non-wetting phase Capillary number of wetting phase Cahn number Characteristic length Effective diameter of drop Fluid-fluid mixing energy Wall energy Free energy of the fluid-wall system Effective bond number Local effective bond number Gravitational acceleration Dimensionless wetting gradient Gas holdup Liquid holdup Height of domain Homogeneous gas fraction Intrinsic permeability Relative permeability of ith phase Relative permeability of non-wetting phase Relative permeability of wetting phase Diffusion length Length of domain Cahn-Hilliard mobility Viscosity ratio xii  Nomenclature n p Q Qi QG QL r R Re Rc Re Rt s S t t∗ u ui um ur ux U VG VL W We x x0 xc xr y  − pascal N.m−2 m3 .s−1 m3 .s−1 m3 .s−1 m m − m m m − − − s m.s−1 m.s−1 m.s−1 m.s−1 m.s−1 m.s−1 m.s−1 m.s−1 m m−1 m − − − m  Unite normal vector Pressure Bulk chemical potential Flow rate through a tube of radius R Gas flow rate Liquid flow rate r component of Cylindrical coordinates Radius of tube Reynolds number Pore chamber radius Pore effective radius Pore throat radius positive number Volume fraction of phases time time Velocity magnitude Velocity of fluid ith Characteristic velocity of the external flow r component of velocity field x component of velocity field Velocity field vector Superficial velocity of gas Superficial velocity of liquid Width of domain Wetting gradient x component of cartesian coordinates Instantaneous position (trajectory) of drop center Initial position of drop center Final position of drop center y component of cartesian coordinates  Units (SI) − m degree degree degree  Description Geometrical parameters Capillary width Contact angle Maximum contact angle Static contact angle  Greek symbols Symbol α, β, γ θ θm θs  xiii  Nomenclature Λ λ µ µg µi µl µm µn µw ρ ρg ρi ρl σ σgs σls σw1 σw2 φ ϕ ψ Ω  − N P a.s P a.s P a.s P a.s P a.s P a.s P a.s kg.m−3 kg.m−3 kg.m−3 kg.m−3 N.m−1 N.m−1 N.m−1 N.m−1 N.m−1 − − − m3  Dimensionless diffusion length Mixing energy density Effective viscosity of two phase flow Viscosity of Fluid2 Viscosity of fluid ith Viscosity of Fluid1 Viscosity of medium Viscosity of non-wetting phase Viscosity of wetting phase Effective density of two phase flow Density of fluid2 Density of fluid ith Density of fluid1 Fluid-fluid interfacial tension Gas-solid interfacial tension Liquid-solid interfacial tension Fluid1-solid interfacial tension Fluid2-solid interfacial tension Phase field variable Porosity Distribution function Fluid domain  xiv  Acknowledgements I completed this PhD thesis with the help and support of some outstanding people. First, I would like to express my deeply-felt thanks to my supervisor, Professor James J. Feng, for his excellent and thoughtful support and encouragement during my PhD studies. He showed me how to be an excellent scientific researcher. I learnt from him to be more critical and to dig the problems thoroughly, and not to be disappointed when the hard times of research come. Also, I would like to thank Prof. Zhong-Sheng (Simon) Liu and Prof. Xiaotao (Tony) Bi, my PhD committee advisors, for their great suggestions during my research. Great thanks go to Professor Pengtao Yue and Professor Peng Gao who helped me to learn our in-house codes Amphi and Grummp. I am indebted to my lovely wife, Maryam Hanifehnezhad, for her invaluable support and continual encouragement and love which undoubtedly was one of the most important reasons of every success that I had in my graduate life at UBC. Finally a special and deeply-felt thanks go to my little son, Ali, for his exclusive spiritual support.  xv  Dedication To my father, Moharram, my lovely wife, Maryam Hanifehnezhad, and my little son, Ali, my step-mother, Rogaiieh Bikdel, and in memory of my mother, Zahra Nabavi.  xvi  Chapter 1  Introduction 1.1  Motivation  Interfacial flows through porous media are scientifically and technologically important. Traditionally, these flows in pores of microscopic dimensions are studied in the context of displacing oil by water in enhanced oil recovery [47, 51]. More recently, such flows have taken on new significance in applications such as microfluidics [3, 17, 40] and proton exchange membrane fuel cells (PEMFC) [57, 67, 87]. The small spatial dimensions have several implications. For one, gravity and inertia are often negligible since the Bond number and Reynolds number are much below unity. Meanwhile, capillarity plays a much more important role. This is manifested not only by a curvature-related pressure difference between the two phases, but also by detailed features of the interfacial dynamics, including interfacial deformation and rupture and the motion of three-phase contact lines on the solid wall. This thesis focuses on three aspects of the problem: • Interfacial flow regimes, transitions and hysteresis; • Transport in terms of a relative permeability; • Dynamics of sessile drops driven by external flow and wetting gradients. The background for each will be briefly reviewed below. The gas-liquid flow regimes through channels and microchannels have been summarized into maps in terms of the superficial gas and liquid velocities [8, 50, 70, 71]. In macroscopic pipes and tubes, on the order of 1 cm in diameter, three categories of flow regimes have been reported: separated (e.g., annular and stratified flows), intermittent (e.g., slug and churn flows), and distributed (e.g., bubbly and droplet flows). Microchannels, with a diameter on the order of 1 mm or less, have received much attention recently [35] (Fig. 1.1). The shrinking length scale causes some subtle changes in the flow regimes. For example, in horizontal microchannels, the stratified 1  1.1. Motivation  (a)  (b) Figure 1.1: Schematics of the flow regime maps for microchannels in size of ∼ 1mm  oriented (a) horizontally, (b) Vertically [35]. VL and VG are superficial velocities of liquid and gas phases respectively.  regime is never seen because the importance of gravity is much reduced by the small dimensions [14, 27]. Mishima et al. [52] also reported suppression of the churn flow regime in small channels. This probably reflects the declining role of inertia, as churn flows prevail for the highest flow rates for both the gas and the liquid. Aside from these, the familiar bubbly, slug, churn and annular flow regimes are seen, and the regime map is qualitatively the same as in larger tubes [17, 35]. So far, research on this subject has relied almost exclusively on experimental observations. As the length scale 2  1.1. Motivation shrinks, spatial and temporal resolution becomes increasingly difficult. In this context, it will be very useful to carry out fully-resolved hydrodynamic computation of the fluid dynamics and interfacial evolution on the scale of the pores. Little numerical work of this kind has been done. In terms of two-phase transport through porous media, relative permeability has served as a key parameter. It arises in the context of modelling two-phase transport by using the so-called extended Darcy’s law. Conventional studies on immiscible two-phase flow in porous medium are based on “separate flow models”, where relative permeability of each phase is calculated from an extension of Darcy’s law for single phase flow: ui =  KKri (∇pi − ρi g) µi  (1.1)  where ui , pi , and µi are the velocity, pressure, and viscosity of fluid i respectively, ρi g the body force. K is the intrinsic permeability determined by pore structure of the porous medium alone, and K ri is the relative permeability of phase i, which is often assumed to depend only on saturations of the phases. The extended Darcy’s law assumes that each component is driven by its own pressure gradient and body force and there is no interaction across the fluid-fluid interface. This has been challenged by earlier authors, and shown to be flawed under many practical conditions [5, 23, 42, 48]. Surveying the literature on this topic, therefore, one finds wide discrepancy among measured Kr data. There is not only large quantitative scatter but sometimes qualitative contradictions in the literature [18]. For instance, the assumption of separate, uncoupled single-phase flow implies K rn + Krw 1 [18]. This seems to be confirmed by some measurements [e.g. 15] but violated by others [e.g. 4, 6]. Notably, Ehrlich [23] and Yiotis et al. [79] reported a K rn that varies with S non-monotonically, with a maximum at an intermediate saturation that is greater than 1. There is also a long-running controversy over whether the interfacial tension and viscosity ratio between the two fluids affect the relative permeability [2, 4, 18]. A more general review of the experimental data can be found in [18] and [6]. Our interpretation of this discrepancy is that K r depends on many parameters besides the saturation S of the wetting phase, including the capillary number Ca [31], solid wettability (or contact angle θ) [48], viscosity ratio M [79], initial morphology of the interfaces (or setup of flow: displacement of one fluid by another or pumping both fluids) [2, 48], flow direction (co-current or countercurrent flows) [42], pore geometry and flow history [29]. Dependence of Kr on parameters other than saturation is because of 3  1.1. Motivation interaction among multiphase flow components even at very low Re and Ca. The idea, often called “viscous coupling”, was proposed initially by Rose [62, 63] in the 1970’s. Since then, researchers have used various experimental and computational tools to study how viscous coupling depends on the geometric, physical and operating parameters. Computation of relative permeability has largely relied on simplified geometries exemplified by the tube bundle model [23] and pore-network model [9]. Such calculations assume quasi-static interfaces and postulate local solutions based on Poiseuille flow, and as a result do not account for interfacial dynamics. More recently, pore-scale computation by solving the NavierStokes equations has become possible. The lattice-Boltzmann method has met with considerable success [e.g. 33, 34, 48, 79]. However, resolution of the interface remains a numerical bottleneck because the regular lattice precludes adaptive refinement at the interface. Much remains to be done to clarify the physics underlying the apparently confusing K r data. The last part of my research has been motivated by efforts to optimize the air-water transport through the gas diffusion layer (GDL) in protonexchange-membrane fuel cells (PEMFC). This is a porous medium through which oxygen flows to the reaction site while the water, a product of the reaction, is discharged [57, 75]. Water vapor condenses in the hydrophobic pores, and small droplets are believed to coalesce and then be driven out partly by wetting gradients (Fig. 1.2). Meanwhile, the incoming air creates a counterflow, which tends to drive the water droplets in the opposite direction [32]. How do sessile drops move and coalesce on a substrate bearing a wetting gradient while being subject to an external flow? This is the question that we set out to investigate. A survey of the literature shows that considerable work has been done on drop motion driven by wetting gradients. Early theoretical models [10, 30] balance the driving force—the capillary force due to the asymmetry in contact angles—against the viscous friction to predict the velocity of drop motion. Although details differ, these models share the following general features. Consider a sessile drop on a substrate having a one-dimensional wetting gradient W e = d(cos θ)/dx, θ being the contact angle. The footprint of the drop is assumed to be a circle of radius R. The driving force turns out to be proportional to W eR 2 , while the viscous friction is proportional to the uR. Therefore, the velocity u comes out to be proportional to W eR. The models differ mostly in treatment of the contact line singularity and determination of the viscous friction. One can estimate the friction locally from Tanner’s law [10], or calculate it using Cox’s wedge solution or a lubrication solution [68], and in the end, the only difference amounts to a different 4  1.1. Motivation  Figure 1.2: Water drops formed on the GDL pores of proton exchange fuel cell [56].  prefactor. Recently, Xu and Qian [77] have reported detailed phase-field solutions of the problem. Qualitatively the same results are obtained, even though the contact line is now treated by Cahn-Hilliard diffusion. Experimentally, the idea of moving droplets by wetting gradients was demonstrated by Chaudhury and Whiteside [11]. Ito et al. [38] reported that u increases with W e on a number of substrates, but the dependence is nonlinear. More recent data of Daniel et al. [16] show linear dependence. Moumen et al. [54] measured drop motion on a nonlinear wetting gradient, and concluded that the experimental data confirm theoretical predictions based on balancing the driving force and the friction force. The last two studies have also considered the effects of contact-angle hysteresis. Based on this brief summary of theoretical, numerical and experimental work, the motion of a single droplet on a substrate having a wetting gradient is well understood. In contrast, little can be found in the literature on drop-drop interaction on wetting gradients. The most relevant studies dealt with coalescence between a drop moving on a wetting gradient and a stationary one resting on an adjoining substrate of uniform wettability [46, 76]. Besides, drop motion driven by an external flow has been studied for some time, e.g. Ref. [19, 20, 22, 43, 53, 65, 88], but only on a substrate with uniform wettability. Besides, to our knowledge, there have been no studies of drop coalescence driven by external flow. As it turns out, the dynamics of sessile drops under two opposing driving forces, one due to the wettability gradient on the substrate, and the other due to the external flow, presents many  5  1.2. Objectives and contributions of this research interesting features that cannot be anticipated from the effect of each factor alone.  1.2  Objectives and contributions of this research  A review of the literature shows at best a fragmental knowledge of two-phase flows in micropores. From the perspective of fundamental research, much of the physics rides on the detailed interfacial morphology, which is unfortunately not easy to observe and record in experiments. On the practical side, empirical and sometimes ill-grounded rules have been used in place of sound design principles. Thus, my thesis has identified the following as its general objectives: • Identify the flow regimes in general two-phase flows in micropores; • Understand the transitions between the flow regimes; • Establish a connection between the relative permeability and the flow regimes and morphology; • Confirm or refute the “extended Darcy’s law” as the basis for engineering design in porous media; • Establish the fundamental physics of drop motion and coalescence under simultaneous influence of wettability gradients and an external flow. As will be seen, these have been accomplished to a good degree in the research described in Chapters 3–5. The main contributions of the thesis can be summarized into four general items: • We have constructed a flow-regime map covering a reasonable range of parameters, with the transition between regimes, including hysteresis, identified and associated with the interfacial morphology. • We have determined that the pressure-flow relationship through a porous medium is generally nonlinear for two-phase flows, and the extended Darcy’s law does not hold in general. Mechanisms of phase interaction have been delineated based on pore-scale hydrodynamics. • We have discovered novel dynamics of sessile drops driven by substrate wetting gradient and an external flow. Hysteresis plays an important 6  1.3. Thesis outline part, and drops may separate or merge depend on the relative strength of the two forcing agents. • Our new understanding of the close association between macroscopic measurable quantities and pore-scale interfacial dynamics argues for a more rational design strategy based on flow-dynamic principles and high-resolution data, as opposed to empirical rules lacking general applicability.  1.3  Thesis outline  Chapter 2 presents details of the research methodology. We present a formulation of the Cahn-Hilliard diffuse-interface model, which will be used for all the simulations to be presented in the thesis. This model has been implemented in our in-house software toolkit AMPHI. This toolkit has two major ingredients: A flow solver which is based on Petrov-Galerkin finite-element algorithm and an adaptive meshing scheme. Additionally, we explain the protocol on how to extract physical results from the numerical data. Chapters 3–5 deal with the three research projects in turn. Chapter 3 presents detailed results of the simulation of gas-liquid two-phase flows in microchannels periodically patterned with grooves and ridges. Depending on the pressure gradient, capillary force and the level of liquid saturation, a number of flow regimes may appear: gas flow, blockage, liquid flow, bubbleslug flow, droplet flow, annular flow and annular-droplet flow. Gas flow and blockage flow regimes are new. We construct a flow regime map that summarizes these regimes. In addition, we report significant hysteresis on boundaries between regimes on the map; this is largely owing to the pinning of the interface at sharp corners in the flow conduit. Under the same operating conditions, different flow regimes can be realized from different initial conditions. Chapter 4 presents our computation of the relative permeability for gasliquid two-phase flows through a model porous medium made of corrugated tubes. This is computed by averaging over a pore-size-distribution of a real porous medium. A constant body force is applied on both fluids to simulate a pressure-driven creeping flow. Our results show that the flow rates vary nonlinearly with the pressure gradient, and the extended Darcy’s law does not hold in general. The interaction between the two phases, known as viscous coupling, is a prominent feature of the process. As a result, the relative permeability depends not only on saturation, but also on the capillary number, viscosity ratio, wettability of the solid wall, pore 7  1.3. Thesis outline geometry, and the initial configuration. The effects of these factors are explored systematically and compared with previous studies. Chapter 5 presents simulations of drop motion and coalescence on a substrate subject to a wetting gradient and an external flow. For a single drop or a pair of identical drops, we consider three scenarios of driving their motion: by a wetting gradient on the substrate, by an external flow due to a pressure gradient parallel to the substrate, and by a combination of the two. Results show that both mechanisms for driving the drop depend intimately on the initial configuration and there is considerable hysteresis in the drop dynamics. We examine how co-current and counter-current forcing induces drop coalescence or separation under different conditions. The final chapter, Chapter 6, summarizes the key results of the thesis, outlines the limitations of the current work, and makes recommendations for future work.  8  Chapter 2  Methodology Of Research 2.1  The diffuse interface model  The prominent features of the fluid dynamic problems include the deformation and movement of the interface, the three-phase contact line, and the surface wettability. These will be handled in a diffuse-interface framework using the Cahn-Hilliard model based on free energy formulation [39, 86, 90]. The three-phase contact line is one of the most significant unresolved problems in multiphase flows computations. Consider the steady-state displacement of fluid 2 by fluid 1 on a solid wall, Fig. 2.1, where the solid wall is stationary. A straightforward formulation of the problem runs into a singularity at the contact line, where based on no-slip velocity conditions it requires pinning of the interface. The diffuse-interface framework has proved to be convenient for interfacial problems, especially in regularizing the singularity of the three-phase contact line. The fluid interface is treated as a thin but diffuse layer where the two components mix to some extent. Thus, even with the no-slip condi-  Figure 2.1: A static contact line viewed in (a) the sharp-interface model and (b) the diffuse-interface model. n is the outward normal to the wall. The phase field variable φ represents the two fluid bulks by ±1 and the fluid-fluid interface by φ = 0. Other symbols are defined in the text. Reprinted from Ref. [80] by permission.  9  2.1. The diffuse interface model tion on solid walls, a contact line may move thanks to mass flux across the interface through diffusion. We introduce a scaled “concentration” φ such that in the two-fluid bulks φ = ±1 and interface is given by φ = 0. φ varies steeply but continuously across the interface. Denoting the fluid domain by Ω and the solid surface by ∂Ω, the free energy of the fluid-wall system F can be written as F=  Ω  fm (φ, ∇φ)dΩ +  fw (φ)dA,  (2.1)  ∂Ω  where Ω is the fluid domain and ∂Ω is the solid surface. f m is the fluid-fluid mixing energy and can be written in Landau-Ginzburg form [74]: fm (φ, ∇φ) =  λ λ |∇φ|2 + 2 (φ2 − 1)2 , 2 4  (2.2)  and fw is the wall energy [39, 73]: fw (φ) = −σcosθs  φ(3 − φ2 ) σw1 + σw2 + , 4 2  (2.3)  where λ is the mixing energy density, and is the capillary width representative of the thickness of interface. The gradient energy term λ2 |∇φ|2 and the bulk energy term 4λ2 (φ2 − 1)2 represent the philic and phobic tendencies between the fluid species, their competition ends to the equilibrium φ profile. As fm is the diffuse interface counterpart of the interfacial tension, so we can relate the conventional interfacial tension σ to the parameters in the mixing energy. For instance, from an equilibrium hyperbolic-tangent φ profile that is the 1D energy minimizer, we obtain [85]: √ 2 2λ σ= (2.4) 3 The wall energy is designed so that away from the contact line, f w gives the fluid-solid interfacial tensions σ w1 and σw2 for the two fluids, which determine the static contact angle θs through Young’s equation σw2 − σw1 = σcosθs . In equilibrium, minimizing mixing energy F shows that φ contours are parallel lines intersecting the wall at the angle θ s . Along the normal to the interface, φ retains its characteristic hyperbolic-tangent profile, undisturbed by fw . A variational procedure of δF δφ leads to the bulk chemical 2  potential Q = λ −∇2 φ + φ(φ 2−1) [39, 69], and the gradient of Q controls the diffusion and advection of φ through Cahn-Hilliard equation ∂φ + U · ∇φ = m∇2 Q, ∂t∗  (2.5) 10  2.2. Cahn-Hilliard model parameters where m is the Cahn-Hilliard mobility, and U is velocity field. Now the moving contact line can be formulated using the Cahn-Hilliard equation (2.5) and Navier-Stokes equations ((2.6) and (2.7)). ∇·U  ρ(  = 0,  (2.6)  ∂U + U · ∇U ) = −∇p + ∇ · [µ(∇U + ∇U T )] + Q∇φ + B, (2.7) ∂t∗  where t∗ is time, p is pressure, and Q∇φ is the diffuse-interface representation of the interfacial tension. For the flow in microscopic pores of interest here, and droplets of submillimeter diameters the Bond number is much bellow unity. Thus we have neglected gravity [44, 87]. B is a constant body force that is used to impose a pressure gradient on the fluids. The effective viscosity µ and effective density ρ are defined as an average between that of the fluids weighted by the volume fractions (1 + φ)/2 and (1 − φ)/2: µ=  1+φ 1−φ µl + µg , 2 2  (2.8)  1+φ 1−φ ρl + ρg , (2.9) 2 2 Now the gas-liquid interface is no longer a boundary that requires boundary conditions. ρ=  2.2  Cahn-Hilliard model parameters  The diffuse-interface model introduces two “mesoscopic” parameters, the Cahn number Cn = /D and the parameter Λ = l d /D. The former is the ratio between the interfacial thickness and the macroscopic characteristic length, while the latter between the diffusion length l d = (µl µg )1/4 m1/2 and D [80, 84]. These must be chosen judiciously. Cn should be small enough so the sharp-interface limit is approached. The diffusion length, or the CahnHilliard mobility m, is directly related to the motion of contact line, which, as mentioned earlier, happens via interfacial diffusion. A sharp-interface limit with zero diffusion would produce a pinned contact line. Instead, we should seek a limit with a finite diffusive flux across the interface. Such a limit is obtained by keeping the Cahn-Hilliard mobility m fixed while reducing value. Yue et al. [84] has presented a scaling arguments in order to rationalize independency of bulk chemical potential Q from . They showed that the 11  2.2. Cahn-Hilliard model parameters  Figure 2.2: Convergence of the moving contact line to the sharp-interface limit with decreasing Cn, shown for steady interfacial shapes in a Couette flow at three capillary number Ca values. Walls contact angle θs = 90. The m value corresponds to Λ = ld /D = 0.01. x is flow direction and walls are at y = 0 and y = 1. Adapted from Yue et al. [84] with permission.  chemical potential Q varies not over but a larger diffusion length l d = (µl µg )1/4 m1/2 , and that the variation of Q across the contact line scales with (µl µg )1/2 u/ld . Thus, as → 0 at a fixed m, both ld and the variation of Q stay constant, resulting in a finite diffusive flux that, of necessity, produces the slip velocity u. To support this sharp-interface limit conceptualization, Yue et al. [84] used numerical evidence. They computed the steady-state solutions for an interface sheared between parallel plates or driven by a pressure gradient in a capillary tube. Over all the capillary numbers tested, the interface approaches a definite limit as → 0 while all other parameters, especially the mobility m, are fixed (Fig. 2.2). Beside this the other important consideration is the implications of Cahn-Hilliard diffusion on mass conservation when using a phase-field model for simulating two-phase flows. Although globally the phase-field variable φ is conserved, a drop may shrink spontaneously while φ shifts from its expected values in the bulk phases. The shift of the phase-field variable φ and shrinkage of drops are phenomena inherent to the Cahn-Hilliard dynamics. The magnitude of φ shifting as well as drop shrinkage are proportional to the Cahn number Cn = /D. In particular, the mass loss becomes negligible if a small enough capillary width is used [83]. 12  2.3. Discretization of computational domain Based on the procedure proposed by Yue and Feng [80] for determining the Cahn-Hilliard model parameters to simulate moving contact lines, first we choose the smallest manageable . Then λ is calculated using equation (2.4). Based on the criteria Cn ≤ 4Λ [80], we choose Λ to reach the sharp-interface limit which is in the range 10 −5 ≤ Λ ≤ 10−2 as recommended in Ref. [85]. In the simulations presented in this thesis, I have used Cn ≤ 0.02 and Λ ≤ 0.01, which according to our previous studies are sufficient for the sharp-interface limit [28, 84].  2.3  Discretization of computational domain  To get accurate results of interfacial coalescence, rupturing, and deformation, one needs not only a proper physical model for the interface, but also very fine spatial resolution of the interface. However, using a very fine mesh throughout the whole computational domain is expensive, sometimes prohibitively so depending on the flow geometry and scale of the problem. Adaptive mesh generation/refinement is a powerful technique to enable accurate simulations are manageable costs. In this technique a high quality dense mesh is produced around the interface while the mesh density is low far away from the interface. In our calculations such adaptive meshing is achieved by using a general-purpose mesh generator GRUMMP, which stands for Generation and Refinement of Unstructured Mixed-Element Meshes in Parallel [26]. GRUMMP produces triangular elements in 2D and tetrahedral elements in 3D and controls the spatial variation of grid size using a length scale LS, which specifies the intended grid size at each location in the domain. This length scale is related to three parameters h 1 , h2 , and h3 which control the density of mesh in the interface and inside the bulks respectively. A grading factor GF is used to produce generally satisfactory transitions among different regions of the mesh. Fig. 2.3 shows an example of the mesh inside a square containing an ellipse. The basic schemes in GRUMMP follow the work of Ruppert [64], with several significant improvements in the areas of cell size and grading control by Oliver-Gooch and Boivin [61]. In the problems simulated here, we have done mesh-refinement tests as well to ensure that the numerical results have converged with the grid size.  2.4  Numerical algorithm and validation  We solve the governing equations with the boundary conditions, which will be laid out for each problem in the ensuing chapters, using an in-house finite13  2.4. Numerical algorithm and validation  Figure 2.3: (a) An unstructured triangular mesh generated by GRUMMP with interfacial refinement. The parameters are GF = 3, outer boundary mesh size h2 = 0.5, interior mesh size h3 = h2/2, and interfacial mesh size h1 = h2/64. (b) Magnified view of a portion of the interfacial region. The bold curve indicates the interface φ = 0, which is centered in a belt-like region of refined triangles. Reprinted from Ref. [85] by permission.  element package AMPHI [86, 90]. Our numerical algorithm has two major ingredients: a finite-element flow solver and an adaptive meshing scheme. The discretization of the governing equations follows weak formulation of the standard Galerkin method [36]. However, the Cahn-Hilliard equation requires special attention. With C 0 elements, which are smooth within each element and continuous across their boundaries, one cannot represent spatial derivatives of higher order than 2. Thus the fourth-order CahnHilliard equation (2.5) has been decomposed into following two second-order equations using variable Ψ [85]. ∂φ mλ + U.∇φ = 2 ∆(Ψ + sφ), ∂T  (2.10)  Ψ = − 2 ∆φ + (φ2 − 1 − s)φ,  (2.11)  where s is a positive number that enhances the convergence of the iterative solution of the final linear system. On an unstructured triangular mesh, piecewise quadratic (P2) elements for velocity and piecewise linear (P1) elements for pressure have been used. For temporal discretization, the Crank Nicholson scheme and the three-point backward difference scheme have been used [36], and these second-order 14  2.4. Numerical algorithm and validation implicit schemes give nearly identical results. After spatial discretization of the equations, all equations have been listed in ref. [85], the linear system within the Newton loop is solved using preconditioned Krylov methods such as the generalized minimum residual (GMRES) method and the biconjugate gradient stabilized (BCGSTAB) method. Yue et al. [86] have presented detailed descriptions of the numerical method and performed extensive validations by calculating benchmark interfacial problems, including drop deformation in shear flow, contraction of a torus and drop retraction. In the last problem, these authors considered a spheroidal drop that retracts toward a sphere under interfacial tension in a quiescent matrix fluid. This process has received much attention as an experimental procedure for measuring the interfacial tension between the two phases. The simulations of Yue et al. [86] using AMPHI was in excellent agreement with experiments and with predictions of the Maffettone and Minale model [49]. Another notable benchmark problem is the rise or sedimentation of a drop in an immiscible ambient fluid. There is a wealth of experimental data on this process. Using AMPHI, Yue et al. [86] simulated the rise of Newtonian drops in a Newtonian matrix where inertia plays a major role. The terminal velocity and the drop shape were in very good agreement with those reported in the literature. More computations on validation of our numerical toolkit can be found in Ref. [80, 82, 89, 90]. In view of the extensive validation on AMPHI that has been done before, I will present limited validation in the following chapters, and instead focus on the novel physics of the new problems.  15  Chapter 3  Flow regimes in corrugated microchannels 3.1  Introduction  As discussed in chapter 1 the gas-liquid flow regimes have been summarized into maps in terms of the superficial gas and liquid velocities [8, 50, 70, 71]. The shrinking length scale from macro to micro causes some subtle changes in the flow regimes. For example, in horizontal microchannels, the stratified regime is never seen because the importance of gravity is much reduced by the small dimensions [14, 27]. Mishima et al. [52] also reported suppression of the churn flow regime in small channels. This probably reflects the declining role of inertia, as churn flows prevail for the highest flow rates for both the gas and the liquid. Aside from these, the familiar bubbly, slug, churn and annular flow regimes are seen, and the regime map is qualitatively the same as in larger tubes [17, 35]. The flow-regime maps in microscopic channels and tubes are important and useful. But these are empirical observations only, and by themselves do not offer a clear understanding of the underlying fluid dynamics. How does the interfacial dynamics on the pore scale determine the overall flow regimes? What are the fundamental physical principles that operate over different geometries? Such insights are difficult to obtain from microscale experiments because of limitations on spatial and temporal resolution. Thus, we turn to numerical simulations. Our work differs from previous researches in several aspects. First and foremost, we strive for a fully-resolved hydrodynamic computation of the fluid dynamics and interfacial evolution on the scale of the pores. The aim is to establish an understanding of the mechanisms producing the flow regimes from first principles. This contrasts with the existing literature that relies almost exclusively on experimental observations. Second, we adopt a geometric and mechanical setup of the problem relevant to water transport in the gas-diffusion medium of PEM (proton ex-  16  3.2. Problem setup change membrane) fuel cells [21, 56], which differs from previous studies. For one, we adopt a periodically corrugated axisymmetric flow conduit to simulate the geometry of complex pores with pore chambers and pore throats. In addition, we consider pore sizes on the order of 10 microns [66], much smaller than the tubes and pipes of existing studies. As a result, interfacial tension is expected to dominate gravity and inertia. Moreover, the flow is driven by an external pressure gradient applied to both components. This way, the gas and liquid flow rates are not control parameters but outcomes of the fluid dynamics. Finally, we give special attention to the wetting properties of the solid walls. In gas-diffusion medium, hydrophobic coating is frequently applied to facilitate water removal [45]. In the two-phase flow literature, a few authors have experimented with modifying the wettability of the smooth surface of tubes [7, 13].  3.2  Problem setup  As mentioned in Chapter 2 for the flow in microscopic pores of interest here, estimation from typical pore sizes and operating conditions in the gas diffusion medium of PEM fuel cells gives Reynolds and Bond numbers on the order of 10−2 [44, 87]. Thus, we neglect gravity and inertia and focus on creeping flows in which capillary forces play an important role. We consider flow through the periodic geometry of Fig. 3.1 driven by a constant pressure gradient. It is convenient to replace the prescribed pressure drop over one period of the channel by a constant body force B acting on both components along the flow direction such that the pressure will be periodic as well as the velocity field, and the magnitude of the effective body force B will be a key control parameter for the flow problem. The initial condition typically has the liquid and gas at rest in the domain, with a certain initial configuration for the gas-liquid interface. For example, the liquid may initially occupy the pore chamber as shown in Fig. 3.1. This initial configuration specifies the liquid saturation S in the pore, defined as the volume fraction of the liquid relative to the entire pore volume. The saturation S is an important control parameter of the problem, and remains constant for the duration of each flow simulation. Under the constant driving force B, both components start to flow and eventually a steady, or time-periodic flow pattern establishes itself. Periodicity is imposed on U , p and φ between the left and right boundaries of the axisymmetric computational domain depicted in Fig. 3.1. The top of the domain is the axis of symmetry on which we impose 17  3.2. Problem setup  Figure 3.1: Schematics of the periodic flow geometry. The top of the domain is the axis of symmetry, and periodic boundary conditions are imposed between the left and right sides. Initially the liquid rests in the groove representing a pore chamber.  symmetry conditions ∂/∂r = 0 and Ur = 0. On the solid substrate, the following conditions are used: U  = 0,  (3.1)  n · ∇Q = 0,  (3.2)  λn · ∇φ + fw (φ) = 0,  (3.3)  where n is the unit normal vector pointing from the fluid into the wall. Of these, the first condition asserts zero slip and that the contact line motion is solely due to Cahn-Hilliard diffusion. The second denotes zero flux through the walls. The third is the diffuse-interface specification of the wetting angle. It is the natural boundary condition arising from the variation of the wall energy fw : φ(3 − φ2 ) σls + σgs + , (3.4) fw (φ) = −σ cos θs 4 2 which is the interfacial energy between the fluids and the solid substrate [39, 84]. In the fluid bulk phases, fw recovers σls and σgs , the liquid-solid and gas-solid interfacial tensions, which determine the static contact angle θ s via Young’s equation σ cos θs = σgs − σls . Equation (3.4) implies an equilibrium between the fluid components and the wall, such that the dynamic contact angle remains at the static value up to the leading order [80, 84]. The most important dimensionless control parameter is F = BD 2 /σ, where the characteristic length D is chosen to be the throat radius (R t ). This can be seen as an effective Bond number as B is an effective body force. Alternatively, if one thinks of the pressure gradient as producing a characteristic speed for the viscous flow in the channel, F may be viewed as a 18  3.3. Results capillary number. Other dimensionless parameters are the liquid saturation S, liquid-gas viscosity ratio M = µl /µg , and the length ratios α = Lc /Lt , β = Rc /Rt and γ = Lt /Rt . Additionally, two other parameters, the Cahn number Cn = /D and the parameter Λ = ld /D, come from the diffuse-interface model. In the simulations presented here, we have used Cn = 0.01 and Λ = 0.01 , as discussed in chapter 2.  3.3  Results  A large number of flow regimes have been predicted in our simulations, and hysteresis is a prominent feature in the transition between regimes. Related to hysteresis is the observation that under identical material, geometric and operating conditions, different regimes can be reached starting from different initial conditions. To present this rather complex picture in a systematic way, we will first describe all the predicted regimes reached from a certain initial condition at different operating parameters S and F . A flow regime map can thus be constructed in the (F, S) plane. Next, we will discuss hysteresis in the sense of dependence on the flow history, with F being decreased or increased gradually for a fixed S, and on the initial condition. This history dependence turns out to be closely related to our corrugated geometry, and we will explore the effect of replacing the sharp corners by rounded ones. Finally, we examine the role of solid wettability on the flow pattern.  3.3.1  Flow regimes  First we construct a map of flow regimes for a “baseline” setup of the problem. We fix the geometric parameters at α = 1, β = 2 and γ = 2, the wetting angle at θs = 135◦ , and the liquid-gas viscosity ratio at M = 18. The θs and M values correspond to typical operating conditions in the gasdiffuse medium of PEM fuel cells, with water and air at 80 ◦ C. Initially the liquid sits entirely in the pore chambers, as depicted in Fig. 3.1. The maximum realizable liquid saturation is then constrained by the depinning of a contact line at the sharp corner, or by the liquid meniscus reaching the axis of symmetry and thus forming a liquid bridge. As illustrated in Fig. 3.2, the Gibbs criterion for pinning at the 90 ◦ corner requires [e.g., 60]: θs − 90◦ < θ < θs .  (3.5)  19  3.3. Results  Figure 3.2: Schematic illustration of the condition for contact line pinning. The solid curve is the gas-liquid interface, and the dashed curves indicate its limiting positions before the contact line depins from the corner.  With the large contact angle θs = 135◦ in our case, the maximum saturation is S = 0.78 corresponding to the geometric condition of the liquid arc reaching the top of the computational domain. With less liquid, the interface lowers and may depin from the corner. So the minimum S is zero. At each prescribed S, the interface initially assumes a shape of a circular arc. Then a constant body force F is applied to initiate the flow, which evolves into a well-established pattern that is recorded as a regime. For the baseline setup of the problem, all the flow regimes predicted are depicted in Fig. 3.3 in the F –S parameter space. For the lowest values of the driving force F , one of two states prevails depending on the level of saturation. For low S, the liquid mostly stays within the pore chamber or groove of the channel. The gas flow in the middle is not strong enough to drive the liquid into the pore throat. This is the regime of gas flow, shown in Fig. 3.4. Pinning of the interface at the downstream corner of the pore chamber prevents the liquid from wetting the pore throat. For higher saturation, even the gentle driving force is able to cause the liquid meniscus to raise a crest, which merges in the center of the conduit to form a liquid bridge that completely blocks the gas flow. This blockage regime is depicted in Fig. 3.5. Though not obvious in the plot, the liquid bridge is slightly asymmetric fore and aft, with differing pinning angles at the corners. The net force due to interfacial tension balances the driving force F . 20  3.3. Results  Figure 3.3: Flow regime map for the baseline model. Note that the columns of symbols are spaced for clarity and not according to the exact S values. With increasing F , two routes of evolution are possible depending on S. For the lowest saturation (S = 0.2 in Fig. 3.3), the liquid remains in the groove even for the strongest driving force. The gas flow regime persists. For higher S, a series of flow regimes are predicted with increasing F . We will illustrate them by using S = 0.38. The first transition gives rise to the liquid flow regime depicted in Fig. 3.6. The liquid is drawn from the pore chamber into the throat and forms a liquid bridge spanning the cross section much like in the blockage regime. But now the greater F deforms it and moves it downstream until it merges with the liquid bridge of the next period. This produces a continuous liquid stream in the middle of the pore, with a gas pocket trapped in the chamber. At higher S, the blockage regime of Fig. 3.5 gives way to the liquid flow regime via a similar process. A still higher F produces the bubble-slug flow regime (Fig. 3.7). In the literature, bubbles are sometimes defined as gaseous blobs having an effective diameter below 75% of the tube diameter, and slugs as larger ones [17]. In 21  3.3. Results  Figure 3.4: The gas flow regime achieved for a low saturation (S = 0.38) under a gentle driving force F = 0.25. The meniscus is circular at the start of the simulation, with no flow (t = 0). It deforms moderately in the steady state (t = 43.6). The upper edge of the plots corresponds to the axis of the pore. Time is scaled by µl Rt /σ.  Figure 3.5: Development of complete blockage in the microchannel at F = 0.001 and S = 0.75.  our context there is no need to distinguish them. This regime resembles the liquid flow regime except for the discrete gas bubbles or slugs carried along the center of the conduit. Compared with Fig. 3.6, the stronger F elongates the liquid bridge axially so as to create a lamella oriented more or less along the flow direction (t = 12.7). When the undulation of the lamella causes a secondary coalescence in the center, a gas bubble or slug is entrapped (t = 17). The flow evolves into a time-periodic pattern with the passing of the bubble or slug through the pore. The gas pocket entrapped in the corner of the chamber maintains an essentially constant shape with recirculation inside. If F is increased further, the annular flow regime prevails for most of the saturation values tested (Fig. 3.8). Similar to the bubble-slug regime, the liquid moves out of the pore chamber into the pore throat. The fast gas flow spreads the liquid lamella downstream, keeping it close to the solid ridge and preventing the formation of a liquid bridge. When the lamella connects with the remaining liquid in the chamber (t = 10.9), a continuous liquid annulus is formed that encloses the gas flow in the core. Two variations of this are the droplet flow and annular-droplet flow regimes, which appear respectively for very low and very high saturations. At low S, there is not enough liquid to form a continuous annulus. Instead, some of the liquid is torn off by the gas stream to form droplets (Fig. 3.9a). At higher S, the 22  3.3. Results  Figure 3.6: Development of the liquid flow regime with F = 0.75 and S = 0.38. Because of the axisymmetry of the geometry, as the liquid moves from the chamber into the pore throat, its projected area on the meridian plane being plotted increases even though its volume is conserved.  Figure 3.7: Development of the bubble-slug regime at F = 1.75 and S = 0.38. abundance of liquid is such that both the annulus and droplets are created (Fig. 3.9b).  3.3.2  Hysteresis and dependence on initial configuration  In constructing the map of flow regimes in Fig. 3.3, we have always started from the initial condition depicted in Fig. 3.1, with a prescribed amount of liquid sitting in the pore chamber. With zero initial velocity, a prescribed body force is suddenly applied to both components, and the flow evolves till a robust pattern emerges, which may be steady or time-periodic. By hysteresis, we refer to how the flow regime realized depends on the history of varying the control parameters. If we start inside a certain regime, say annular flow, and decrease F in small steps, will we reach the bubbleslug flow and the other regimes successively that lie below in the map? What if we start from the gas flow regime at a small F , and then gradually increase F ? Will the boundary between regimes shift depending on whether a parameter is increased or decreased? If we start from different initial states, and suddenly impose the same control parameters, will we arrive at different final states? To probe these questions, we have tested two saturation values, S = 0.31 and S = 0.71. In each case, we ramp up the flow by increasing F gradually starting from an initial condition in the gas flow or blockage regimes. This is called the forward path. We then ramp down the flow, starting with the end state of the forward path, by decreasing F gradually. This will be the backward path. Figure 3.10 compares the forward and backward paths with 23  3.3. Results  Figure 3.8: Development of the annular flow regime at F = 3 and S = 0.38.  Figure 3.9: Variants of the annular flow regime: (a) droplet flow at F = 3 and S = 0.31; (b) annular-droplet flow at F = 3 and S = 0.75. Both are timeperiodic configurations, and the liquid film on the ridge in (a) has stopped moving downstream.  the regimes realized from the baseline model (i.e., starting from the standard initial condition of Fig. 3.1) taken from the regime map of Fig. 3.3. While the baseline model predicts multiple regimes, the forward path encounters only one transition: from gas flow (for S = 0.31) or blockage (for S = 0.71) to the liquid flow regime. Then liquid flow persists, not only till the highest F on the forward path, but also on the entire backward path, down to F as small as 10−3 . Therefore, the flow regimes show an exceedingly strong dependence on the history of changing F , so much so that the term “hysteresis” becomes somewhat inappropriate as it typically refers to the delay of a transition, not its utter absence. The dominance of the liquid flow regime can be rationalized from the interfacial morphology. Once the pore throat is completely filled with liquid and the gas is sealed within the pore chamber (last frame of Fig. 3.6), the throughput of the pore is entirely liquid. Varying the driving force F mostly changes the liquid flow rate, with a minor effect on the shape of the interface and the recirculation within the gas. In particular, the interface is pinned at the upstream corner, and there is no mechanism to liberate the gas from the confinement of the pore chamber. Therefore transition to other flow regimes is not possible. This contrasts with the setup of the baseline model, with gas initially in the central region of the pore. To further demonstrate the effect of the initial configuration on the final flow regimes, we start with a liquid bridge spanning the entire cross section of the pore chamber at different axial positions (first frames of Fig. 3.11a,b,c), 24  3.3. Results  (a)  (b)  Figure 3.10: Hysteresis or path dependence of the flow regimes at saturation S = 0.31 (a) and S = 0.71 (b). Other parameters are the same as have produced the flow regime map in Fig. 3.3.  all having the same saturation S = 0.38. When the same body force F = 3 is imposed, the flow evolves into three different regimes: droplet, bubbleslug and annular-droplet flows. None is identical to the annular flow regime obtained from the baseline setup (Fig. 3.8). In these particular cases, the final outcome depends on whether the initial liquid bridge adheres to parts of the chamber walls. In Fig. 3.11(a ), a liquid sheet remains attached to the upstream corner of the pore chamber and coats the pore throat as well. Figure 3.11(c) is similar but has the liquid sheet attached to the downstream side of the chamber. In Fig. 3.11(b), the liquid bridge never wets these parts of the chamber; it is readily stretched into the throat and eventually detaches entirely from the base of the chamber. Although these details are specific to the initial morphologies tested, their implication is general: initial morphology can affect the flow regimes  25  3.3. Results  (a)  (b)  (c) Figure 3.11: Dependence of flow regimes on the initial configuration. (a) Droplet flow develops from a liquid bridge initially wetting the upstream wall of the pore chamber. (b) Bubble-slug flow develops from a liquid bridge initially in the middle of the pore chamber. (c) Annular-droplet flow develops from a liquid bridge initially attached to the downstream wall of the pore chamber. F = 3, S = 0.38, and the other conditions are the same as in the baseline model. achieved under identical operating conditions. We have tested other initial conditions at different operating conditions with the same upshot. For brevity we will not show the detailed results. An overarching observation is that the hysteresis and missing flow regimes owe much to the corrugated geometry in our microchannel. The pinning of the interface at the corners of the ridge plays an important role. This explains why prior studies in smooth tubes and pipes of various sizes have not reported hysteresis [71].  3.3.3  Geometric effect: rounded corners  Given the prominent role of the corrugated geometry in determining the flow regimes and transitions, a natural question is what if we remove the sharp corners and make the constriction gradual and smooth. This has been tested by replacing the solid ridge in the baseline model, with α = 1, β = 2 and γ = 2, by a rounded (semi-circle) corners. All other conditions are the same as the baseline model with sharp corners. Starting from the same initial morphology with liquid resting in the pore chamber, imposing different F values leads to different regimes in the rounded-corner geometry. These regimes are compared in Fig. 3.12 with those in the sharp-corner geometry for a fixed S = 0.38. Whereas gas flow gives rise to liquid flow when F exceeds 0.5 with sharp corners, the rounded geometry produces a core-annular morphology with a liquid core enclosed in a gas shell, called the shell flow regime (Fig. 3.13). Absent the 26  3.3. Results sharp corner to pin the interface and anchor the liquid onto the ridge, as shown previously in Fig. 3.6, the liquid eventually detaches from the solid protrusion, which is hydrophobic with θ s = 135◦ . Thus a gas shell is formed that entirely insulates the liquid from the solid walls. This is reminiscent of the lubrication scenario in oil-water core-annular flows [41] and apparent slip on textured substrates [28]. As F increases to 1, the shell flow gives way to a droplet flow regime (Fig. 3.14). In the pore with sharp corners, a similar droplet regime (Fig. 3.9a) occurs at the highest F for relatively low S, above the bubble-slug regime (see Fig. 3.3). With the rounded corners and absent corner-pinning, on the other hand, the liquid detaches from the hydrophobic solid readily at relatively low flow rates. The detachment is complete in the sense that no residue film is left on the solid. As F increases further, the bubble-slug and annular flow regimes appear in succession. These are similar to those predicted in the sharp-cornered pores. To probe hysteresis in the geometry with rounded corners, we again follow a forward path by increasing F gradually, waiting for the flow regime to develop fully at each increment, and then a backward path by decreasing F gradually. Figure 3.15 compares the flow regimes encountered thus with those predicted by the baseline model with the fixed initial configuration. Starting from the gas flow regime (Fig. 3.15a), we have a transition to the shell flow regime at F = 0.5, as in the baseline model. However, no other transitions exist on the forward path. On the backward path, the shell flow regime persists until F = 0.1, below which it gives way to liquid flow. This picture resembles Fig. 3.10(a ), with two differences. First, the liquid flow regime is replaced by the shell flow regime. Second, there is an additional transition on the backward path, to the liquid flow regime. When starting from the droplet regime at F = 1.25 (Fig. 3.15b), no transition takes place on the forward path. On the backward path, the droplet regime persists down to F = 0.75 before turning into shell flow, which in turn transforms to liquid flow below F = 0.1. This last transition is the same as in Fig. 3.15(a ). To summarize the hysteresis in geometry with rounded corners, we see more transitions on the downward path than with sharp corners. In this sense, the dependence on history is lessened in comparison with the sharpcorner geometry. Without interface pinning on sharp corners, liquid flow is no longer the dominant regime. Instead, shell flow or droplet flow becomes more prominent. In all cases, the dominance of a flow regime is due to the fact that once it is formed, further variation of F tends only to vary the liquid and gas flow rate, with little effect on the interfacial morphology. 27  3.3. Results  Figure 3.12: Flow regimes in the baseline model with rounded corners, compared with those taken from the flow map of Fig. 3.3 with sharp corners. S = 0.38 and all other conditions are the same as in Fig. 3.3.  3.3.4  Wettability of the solid  The tendency for the liquid or gas to adhere to the solid wall is fundamentally important to the flow situations being studied. The wetting angle reflects the energetic cost of detaching each fluid component from the solid substrate, and directly affects the shape of the interface and more subtly the motion of the contact lines [80, 84]. The corrugated pore geometry amplifies these effects, e.g., in the pinning of the interface at corners, the detaching of the liquid from the pore chamber and its reattaching to the pore throat. The effect of wettability is highlighted in Fig. 3.16, where three wetting angles are compared in a pore with sharp corners. In the hydrophilic pore 28  3.4. Comparison with previous experiments  Figure 3.13: Development of the shell flow regime with rounded corners at F = 0.5 and S = 0.38.  Figure 3.14: Development of the droplet flow regime with rounded corners at F = 1 and S = 0.38. (a), the liquid spreads along the ridge and eventually forms a continuous film covering the solid, yielding the annular flow regime. In the neutrally wetting case (b), the spreading of the liquid film is slower. More importantly, once the liquid reaches the downstream corner of the ridge, the interface is unable to depin from the corner. Thus, a liquid finger is produced at the corner that cannot merge with the liquid in the pore chamber. A gas flow regime results. Finally, in the hydrophobic pore (c), detachment of a liquid drop occurs shortly after the start of the flow, leading to the droplet flow regime. The liquid attached to the wall assumes a morphology similar to that of (b).  3.4  Comparison with previous experiments  There exists a wealth of experimental observations on gas-liquid flows in smooth tubes and pipes, mostly of larger sizes, but more recently also of diameters below 1 mm. Before comparing our results with those, we note four important differences in the geometric and physical setup of our study. (a) Previous experiments control the gas and liquid flow rates and ascribe flow regimes to these rates. Our setup is such that a prescribed pressure gradient (or body force) is applied to drive both components. Therefore the gas and liquid flow rates are outcomes of the simulation rather than control parameters. 29  3.4. Comparison with previous experiments  (a)  (b)  Figure 3.15: Hysteresis in micropores with rounded corners. At S = 0.38, forward and backward paths are traced by increasing F gradually and then decreasing it, starting from (a) the gas flow regime and (b) the droplet flow regime. The regimes achieved from the baseline model with rounded corners (Fig. 3.12) are shown for comparison. Other parameters are the same as in the baseline model.  (b) Our geometry features a pore chamber and throat, in contrast to the smooth tubes and pipes studied before. This turns out to be an important factor in determining the interfacial behavior and the flow regimes. (c) Previous experiments often employ high flow rates for both components, sometimes with turbulent flows. Therefore, inertia is an important factor. Our computations concern much smaller length scales and inertia is negligible. Thus the flow regimes are all for creeping flows. (d ) Gravity is important in larger tubes, and vertical, horizontal and inclined pipes have to be studied separately. Here gravity has negligible 30  3.4. Comparison with previous experiments  (a)  (b)  (c) Figure 3.16: Effect of wall wettability on the interfacial evolution. (a) Annular flow in hydrophilic pore, θs = 70◦ ; (b) gas flow in neutrally wetting pore, θs = 90◦ ; (c) droplet flow in hydrophobic pore, θs = 135◦ . All other parameters are the same: F = 2, S = 0.6, with geometrical parameters α = 0.5, β = 2, γ = 2. effects, again thanks to the small length scales. A survey of the literature shows that bubbly, intermittent, churn and annular flows are the four basic regimes in smooth straight pipes [35, 70, 71]. For macroscopic horizontal pipes, a fifth regime, stratified flow, also appears. Of these, churn flows and stratified flows are not seen in our simulations. Even in microscopic smooth pipes, the lack of inertial and gravitational effects is known to suppress these regimes [27, 52]. Therefore, their absence in our simulations is not a surprise. All the other regimes appear in our microchannel flow, albeit in somewhat modified forms. Intermittent flows consist in gas slugs or plugs moving in the middle of the pipe, being carried by the continuous liquid phase. In our simulations, we have lumped these configurations and the bubbly flows into our bubble-slug regime. The distinction between “bubbles”, “slugs” and “plugs” is merely the size of the gaseous inclusion. Annular flows are also predicted in our simulations, along with its variant forms of droplet and annular droplet flows. Gas flow, liquid flow, and blockage flow types as well as shell flow regime are predicted here but not in previous studies. To a degree these novel flows have to do with the geometric features of our microchannel. For example, if the minor phase is entirely trapped in the pore chamber, it does not contribute to the total flow rate. Thus, only one component flows, as in 31  3.4. Comparison with previous experiments the gas flow or liquid flow regimes. If the driving force is weak, a liquid bridge covering the entire cross section can be anchored at the corners of the throat, thereby blocking flow of either component. These three flow types are facilitated by the pinning of the interface at sharp corners that preclude one or even both components from flowing. With rounded corners, the shell flow regime appears thanks to the lack of interface-pinning and the hydrophobicity of the solid. Although it is essentially an annular flow with a liquid core enclosed in a gas sheath, it has never been reported previously. In a straight tube, this configuration would be unsustainable because of the Rayleigh instability. As our control scheme prescribes F and S with the gas and liquid flow rates as outcome, it is interesting to compare quantitatively our map of flow regimes with those in the literature in terms of the superficial gas and liquid flow rates. We have attempted such a comparison for the sharp-cornered baseline model. Obviously this is feasible only for those regimes that appear in both settings, i.e., the bubble-slug and annular flow regimes. To obtain the gas and liquid flow rates QG and QL , we time-average one period of the unsteady bubble-slug regime. For steady annular flows, this is unnecessary. Defining an effective pore radius Re based on the volume of void in each period of the pore: Re2 =  αβ 2 + 1 Rc2 Lc + Rt2 Lt = Rt2 , Lc + L t α+1  we can compute the superficial gas and liquid velocities using the effective cross-section area Ae = πRe2 : VG = QG /Ae and VL = QL /Ae . Thus we may convert our dimensionless results into dimensional superficial velocities in the two flow regimes once the physical properties of the two fluids and the pore size are specified. For an experimental benchmark, we have chosen the study of [35]. These authors conducted air-water flow experiments in straight circular tubes of diameters on the order of 1 mm, and integrated their data with those in the literature to produce flow regime maps for horizontal and vertical microchannels of comparable sizes. For air and water flowing through a micropore with Re = 0.5 mm, our flow regimes are compared with the map of [35] for vertical tubes in Fig. 3.17. In our simulations so far, we have lumped bubbles and slugs into one regime. For the purpose of this comparison, we distinguish them by the criterion of [17]; slugs are thus gas blobs longer than 75% of the effective diameter of the channel. 32  3.4. Comparison with previous experiments  Figure 3.17: Flow regime map in terms of the superficial velocities: comparison of our numerical prediction with the experiment of [35]. Symbols are our results and solid lines indicate regime boundaries of [35] separating the regimes indicated at the bottom.  Thus defined, our bubble-slug boundary is somewhat to the right of the experimental boundary. With VG increasing further, the experimenters saw churn flow, which is absent in our inertialess computation. Instead we have annular flow, which occurs in the experiments further to the right. Owing to the many differences in the problem setup, the degree of similarity in Fig. 3.17 is intriguing. Note also that in our parameter space explored here (cf. Fig. 3.3), the gas velocity varies over a much wider range than the liquid velocity. As the same body force is applied to both phases, this discrepancy reflects the differing viscosities as well as the fact that in the two regimes plotted, the liquid stays in contact with the solid while the gas does not. Another important parameter in two-phase flows is the holdup ratio h, which is defined as the ratio of superficial velocities divided by the ratio of volume fractions. The gas holdup, for example, is hG =  VG /VL . (1 − S)/S  Obviously, when two phases are homogeneuously mixed and move with the same local velocity, the superficial velocities scale with the volume fractions and hG = hL = 1. Thus, the holdup ratio reflects the interfacial morphol33  3.5. Conclusion  Figure 3.18: Gas holdup ratio as a function of the homogeneous gas fraction: comparison between our numerical prediction and the experiment of [72]. The computations are for Re = 0.5 mm and the experimental data are for a circular tube of inner diameter 1.09 mm.  ogy and the relative motion of the two phases. Figure 3.18 compares our numerical predictions of the gas holdup with the experimental data of [72]. These authors ran air-water flows in straight circular and non-circular tubes of diameters around 1 mm, and documented the flow regimes, void fractions and frictional pressure drops. From the void fraction and superficial velocities, hG can be obtained as a function of the homogeneous gas fraction, HG = QG /(QG + QL ). Again there is close agreement between the simulation and measurement. In both cases, the gas holdup hG tends to increase with the homogeneous gas fraction HG , a trend well documented in the literature [12].  3.5  Conclusion  This computational study was motivated by the desire to understand the local hydrodynamics of gas-liquid flows through porous medium and microfluidic channels. To make the problem tractable, we simplified the geometry into axisymmetric tubes with periodic constrictions. Thus, we retained the features of pore chambers and pore throats, but disregarded connectivity and branching of pores. This simplification allowed us to carry out a detailed investigation of the interfacial evolution under creeping-flow condi34  3.5. Conclusion tions. Using our diffuse-interface model, we were able to compute interfacial deformation, breakup and coalescence from hydrodynamic principles, and delineate various flow regimes based on the control parameters. Within the parameter ranges explored in this work, the main results can be summarized as follows. (a) Starting from an initial condition with the liquid resting inside the pore chambers, the following flow regimes may appear depending on the pressure gradient, the liquid saturation and the pore geometry: gas flow, blockage, liquid flow, shell flow, bubble-slug flow and annular flow. Annular flow has two variants: droplet flow and droplet-annular flow. These may be represented by a flow regime map on a plane extended by the liquid saturation along one axis and the nondimensionalized pressure gradient along the other. (b) Bubble-slug flows and annular flows are similar to those previously observed in straight tubes and pipes of macroscopic and smaller sizes. But the gas flow, blockage, liquid flow and shell flow regimes are new. Their appearance reflects to a large extent the geometric features of the flow conduit. (c) For bubble-slug and annular flows, the current results are compared with previous experiments. Considering the differences in geometric setup, flow-control schemes and parameter ranges, the agreement is close in terms of the boundaries between flow regimes and the holdup ratio. (d ) The flow regimes are highly sensitive to initial conditions and flow history. Depending on whether the driving force is imposed abruptly or ramped up or down gradually, some of the flow regimes do not appear. This can be considered an extreme form of hysteresis. Different flow patterns can result from different initial conditions under identical control parameters. The dependence on flow history is largely due to interface pinning at corners, and the hysteresis effect is alleviated by rounding the sharp corners of the constriction. (e) Wettability of the solid surface plays an important role in shaping the final flow configuration. Generally, hydrophilicity toward the liquid encourages adherence of a liquid film along the walls, whereas hydrophobicity promotes detachment of the liquid and formation of drops in the axial region. 35  3.5. Conclusion The significance of these results lies in that they form the basis for understanding two-phase flows in complex geometries, and contribute to a more rational treatment of such flows in applications. More specifically, we envision improving the traditional Leverett approach to porous media flows, either by basing it on a firmer theoretical basis or by replacing certain ad hoc elements in it. Toward this aim, we did studied relative permeability of two phase flow in pore scale and probed how the concept of relative permeability can be interpreted and calculated from hydrodynamic principles, and how the capillary pressure, based locally on pore size and curvature of the meniscus, acts to drive liquid flow or resist it until an external pressure effects a capillary breakthrough.  36  Chapter 4  Relative permeability of two phase flow in model porous media 4.1  Introduction  As explained in Chapter 1 conventionally, two-phase flow in porous medium is modeled by generalizing Darcy’s law into a linear relationship between the velocity and pressure gradient for each phase [55], where fluid-fluid interaction across the interface is ignored, and relative permeability of each fluid Kri is taken to be only a function of the average volume fraction in the pores, known as the saturation S. However, as we mentioned earlier, because of the possibility of existence of non-continuse pathways for fluids [1, 6], viscous coupling [42], and transient flows the extended Darcy’s law has found low success. Alternatively, someone may view the extended Darcy’s law as the definition of the relative permeabilities. Then it comes as no surprise that measurement of K ri (S) has produced not only large quantitative scatter but sometimes qualitative contradictions in the literature [2, 4, 6, 15, 18, 79]. An obvious interpretation of this discrepancy is that two-phase transport in porous media depends strongly on the interfacial morphology and fluid dynamics near the interface. In turn these depend on an array of parameters, including the volume fraction of the fluids, pore geometry, capillary number, wetting angle, viscosity ratio, and flow history (e.g., imbibition vs. drainage). It is difficult to capture the interfacial morphology and flow behavior in real porous media. The pores typically have an irregular and complex geometry, and the solid is opaque. Direct flow visualization is therefore impossible. So far, experimental evidence on the flow field has come mostly from model porous media with regularized pore geometry. Recently, more realistic pore-scale computation by solving the NavierStokes equations has become possible. Accurate pore-scale knowledge of 37  4.2. Problem setup and averaging scheme the flow field and interfacial morphology would be the ultimate solution for two-phase transport in porous media. But there are formidable numerical challenges in such computations. First, the interface constitutes an inner boundary whose position has to be solved together with the Navier-Stokes equations for both phases. Second, the interfacial curvature and interfacial force have to be accurately computed using high resolution. Otherwise the interfacial deformation is subject to large errors [81] and even numerical instabilities [78]. So far, the lattice Boltzmann method has met with some success [e.g. 33, 34, 48, 79]. However, resolution of the interface remains a numerical bottleneck because the regular lattice precludes adaptive refinement at the interface. In the present study, we propose an alternative method diffuse-interface for computing the pore-scale flow and the relative permeability of a model porous media. By computing two-phase flows in a collection of such microchannels, the relative permeabilities can be evaluated by averaging over a prescribed pore size distribution. The strength of our method, compared with the lattice Boltzmann method, is the much higher resolution of the interfacial motion. The weakness is the relatively simple representation of the pore geometry. By systematically exploring the effects of saturation, capillary number, wetting angle, viscosity ratio and initial configuration, we have taken a step toward constructing a solid understanding of the interfacial hydrodynamics that underlies the relative permeability.  4.2  Problem setup and averaging scheme  The initial condition typically has the liquid and gas at rest in the domain of Fig. 4.1 with the desired saturation S and a certain initial configuration for the interface. Under the constant driving force B, both components start to flow and eventually a steady, time-periodic or quasi-periodic flow pattern establishes itself. The dimensionless parameters of the problem include an effective Bond number F = BD 2 /σ, which indicates the ratio between the driving force and the surface tension, the saturation S, non-wetting to wetting phase viscosity ratio M = µn /µw , contact angle θ and the length ratios α, β and γ. The length scale D will be specified below. A capillary number Ca can be defined for each phase using its average velocity (Eq. 4.3): Ca i = µi ui /σ. But this will be an outcome of the simulation, not a control parameter. The results consist mainly of the relative permeabilities as functions of these parameters. In addition, two other parameters, the Cahn number Cn = /D 38  4.2. Problem setup and averaging scheme  Figure 4.1: Schematics of the corrugated axisymmetric tube as a model for a pore. The computational domain is half of the meridian plane. The top of the domain is the axis of symmetry, and periodic boundary conditions are imposed between the left and right sides. The geometry is specified by three length ratios: α = Lc /Lt , β = Rc /Rt and γ = Lt /Rt . Initially the liquid forms a disc that blocks the cross-section of the groove representing the pore chamber.  and the parameter Λ = ld /D, come from the diffuse-interface model. Based on discussion in chapter 2 and using pore throat radius R t as characteristic length, we chose Cn = 0.01 and Λ = 0.01 for our computations in this section. We imagine the porous medium as a collection of parallel tubes, each having periodic contraction and expansion as shown in Fig. 4.1. Compared with the real geometry, this model retains the change in cross-sectional area along a pore through chambers and throats, but neglects branching and connectivity between the pores. Based on the volume of the voids, an effective radius can be defined for each corrugated tube: Re =  Rc2 Lc + Rt2 Lt Lc + L t  1/2  = Rt  αβ 2 + 1 α+1  1/2  ,  (4.1)  where α = Lc /Lt and β = Rc /Rt are length ratios. We assume that the collection of tubes is such that Re obeys a prescribed pore size distribution (PSD). To impose a certain saturation S to our model porous medium, we fill each pore with both species of fluid to the prescribed S. This is a simplification since one may imagine larger and smaller pores carrying different volume fractions in reality. Now a pressure drop is applied to the porous medium, and a two-phase flow develops in each pore subject to the common pressure gradient. Because of the periodicity along the flow direction, the saturation in each pore stays fixed at S. All boundary conditions are the same as with flow regime calculations. An averaging scheme is needed to compute macroscopic properties such 39  4.2. Problem setup and averaging scheme  Figure 4.2: Pore size distribution (PSD) of a model porous medium. The distribution function ψ(R) is converted from the data of [44] on pore volume distribution.  as permeability from pore-scale flow quantities. In each tube, we compute the flow until a steady or periodic flow pattern is achieved; in the latter case, the time-averaged flow rate is obtained. To illustrate the averaging among the tubes, let us consider a single-phase flow first. The flow rate through a tube of radius R is: Qi = πR2 ui , ui being the average velocity through this tube. The total flow rate through the collection of tubes is, therefore: Qv =  ψ(R)Qi dR,  (4.2)  where ψ is the PSD normalized such that ψ(R) dR = 1. As the total area A of the porous medium is related to the total pore area by the porosity ϕ: ϕA = ψ(R)πR2 dR, the average velocity over the entire porous medium can be computed: ϕ ψ(R)R2 ui dR Qv u= = . (4.3) A ψ(R)R2 dR Now the intrinsic permeability of the model porous medium, K, can be computed from u and the imposed pressure gradient: K = µu/|∇p|, µ being the viscosity. The same averaging scheme can be repeated for a two-phase flow simulation to yield the effective permeabilities for each phase. The ratio between the effective permeability and the intrinsic permeability gives the relative permeability for each phase. The porosity ϕ will cancel out and will not 40  4.3. Results affect the relative permeability. But it does enter the average velocity for each phase and the capillary numbers. In all results to be presented, we have used the PSD of Fig. 4.2 for a gas diffusion layer (Toray TGP-H-060) of hydrogen fuel cells [44]. The porosity for this porous medium, ϕ = 0.8, is used throughout the paper. The integrals in Eq. (4.3) are computed by summing over increments of the pore size; 25 pore sizes are computed in the range of 5.35–25 µm, corresponding to the measured data points in Fig. 4.2. We have chosen a problem setup with a set of baseline parameters: α = 1, β = 2 and γ = 2, wetting angle θ = 135◦ , and viscosity ratio M = 18. The θ and M values correspond to typical operating conditions in the gas diffusion layer (GDL) of PEM fuel cells, with water and air at 80 ◦ C. The solid wall is hydrophobic; water is the non-wetting phase while air is the wetting phase. For convenience, we will refer to the non-wetting phase as water or liquid, and to the wetting phase as air or gas. The conclusions drawn, of course, will be general and not restricted to this specific pair of fluids. In the following, S or “saturation” refers to that of the liquid, i.e. non-wetting, phase. Initially a liquid bridge of thickness Lc /2 occludes the pore chambers entirely (Fig. 4.1). This initial setup gives a liquid saturation S = 40%, which will remain fixed throughout each computation. In the subsections to follow, we will vary one of the parameters to examine its effects while keeping the other parameters at their baseline values.  4.3  Results  Since our porous medium consists of tubes of different sizes, and the same effective body force B is applied to all of them, the flow regime varies among the tubes. This is conveniently indicated by a local Bond number FR = BRe2 /σ defined for each tube. A number of flow regimes have been predicted with increasing FR : blockage, gas flow, liquid flow, bubble-slug, shell, annular and drop flow (Fig. 4.3). For more details see section 3.3.1. Once we have detailed flow data, we compute the relative permeability Krn and Krw from Eq. (1.1) by averaging over the PSD of a model GDL (Fig. 4.2). In this context, it is necessary to have an overall Bond number F = BD 2 /σ, in which the characteristic length is chosen to be D = 18.25 µm, the effective pore radius corresponding to the maximum of the PSD of Fig. 4.2, that is, the effective radius of the most common pore size. Now the relative permeabilities should be complex functions of the flow parameters [4]: Kri = Kri (S, Can , Caw , M, θ). (4.4) 41  4.3. Results  (a)  (b)  (c) Figure 4.3: Temporal development of three flow regimes for the baseline parameter values: length ratios α = 1, β = 2, γ = 2; liquid-gas viscosity ratio M = 18, liquid saturation S = 40% and wetting angle θ = 135◦ . Time is scaled by µn Re /σ. (a) Blockage at FR = 0.38, when the imposed pressure fails to dislodge the liquid bridge pinned at the corner of the pore chamber; (b) liquid flow regime at FR = 0.56, with the gas entirely trapped in the pore chamber; (c) drop flow regime at FR = 1.4, with large liquid drops being surrounded and conveyed by a continuous gas phase.  In addition, Kri may depend on the geometrical parameters, initial configuration and flow history. In our setup, the key control parameter is actually the driving force F ; the capillary numbers for the flow Ca i come out as a result. We list Cai instead of F in keeping with the conventional protocol in the literature. The dependence on each parameter will be examined below by systematically varying it relative to the baseline setup.  4.3.1  Capillary number effect  The premise for Darcy’s law is the linearity between the flow rate and the imposed pressure drop. In single-phase inertialess flow this holds. A fluid interface that deforms according to viscous and capillary forces introduces a geometric nonlinearity into the problem, and it is important to see if this compromises the linearity underlying the extended Darcy’s law. We vary the imposed body force B and compute the average velocity of each phase and a capillary number: Cai = µi ui /σ. Figure 4.4 plots the liquid and gas capillary numbers as functions of F . For both phases, the variation of Ca with F is decidedly nonlinear. For F < 0.29, the blockage regime prevails even in the largest pores, with neither phase moving (cf. Fig. 4.3a). With increasing F , capillary breakthrough takes place first in the largest pores [21]. This causes a transition to the liquid 42  4.3. Results  Figure 4.4: Dependence of the liquid and gas capillary numbers on the driving force F . All other parameters are at the baseline values given in Fig. 4.3.  flow regime (Fig. 4.3b), and produces a positive Ca n for the non-wetting liquid phase. The gas is still trapped in the pore chambers by the liquid at this stage. Gas flow starts around F = 0.84, when the liquid core breaks up into drops, liberating the gas pocket in the largest pores and causing a transition to the drop flow regime (Fig. 4.3c). Increasing F further activates gas flow in smaller pores, and causes Ca w to rise sharply. For even larger F , the drop flow regime prevails in most pores; further increase in F reduces the size of the liquid drops but produces no new flow regimes. Since the gas is much less viscous than the liquid, Ca w increases more steeply with F than Can . Based on the above, it is obvious that the evolution of interfacial morphology makes the overall flow nonlinear in the porous medium, and the initial rationale for extending Darcy’s law to two-phase flows was illconceived. It is interesting to compare our result with the lattice-Boltzmann computations of [33]. First, they have also observed a blockage regime at weak forcing and a subsequent breakthrough. However, they adopted a smoothed geometry with no interfacial pinning, and started from a core-annular initial morphology with the non-wetting fluid enveloped by the wetting one (see their Fig. 10). In their setup, therefore, it is the non-wetting species that is blocked at weak forcing. This comparison helps to illustrate the important role of the pore geometry and the initial configuration, which will be 43  4.3. Results examined at greater length in Sec. 4.3.5. Furthermore, they have recognized that the evolving interfacial morphology produces a nonlinear Ca ∼ F relationship, in qualitative agreement with our simulation. Finally, in the limit of strong forcing, they predicted drop flow for S = 10% and annular flow for S = 30%. These regimes have been seen in our simulations as well [see phase diagram in [1]. Interestingly, Gustensen and Rothman [33] observed that the Ca ∼ F relationship becomes more or less linear in this limit, since the two phases are flowing roughly in stratified layers with less coupling. This can also be said of our drop flow regime of Fig. 4.3(c). Therefore, a fast linear regime may exist in our Fig. 4.4 as well, though our data do not range to sufficiently large F values to provide clear evidence. The same nonlinearity can be illustrated by plotting the relative permeability of each phase as a function of their respective capillary number (Fig. 4.5). As expected, Kri is not a constant; it varies appreciably with the flow. This is particularly true for the wetting phase. The nearly linear rise of Krw with Caw mirrors the sharp upturn in Fig. 4.4; in the drop flow regime the gas transport increases more rapidly with increasing driving force than the liquid transport. For the non-wetting liquid phase, the initial rise of K rn with Can is due to the progressive formation of a liquid core lubricated by a gas pocket in the pore chamber, similar to Fig. 4.3(c). For higher Ca n , Krn saturates because even in the smaller pores, the liquid becomes fully insulated from the solid wall by a gas cushion, and lubrication cannot be further enhanced. That Krn is nearly flat for Can > 0.04 is consistent with the near linearity of Fig. 4.4(a) for larger Ca n , and can be likened to the fast linear flow regime of [33] mentioned above.  4.3.2  Saturation effect  To study the effect of the saturation level on the relative permeability, we use the same baseline geometry as illustrated in Fig. 4.1. By varying the initial thickness of the liquid bridge, we have computed liquid saturation S ranging from 10% to 80%. Figure 4.6 plots K ri as functions of S. Note that in our setup, the liquid is the non-wetting phase. In straight pores, it is intuitive to think that increasing the liquid saturation will increase the liquid-phase permeability and decrease the gas-phase permeability. In our geometry with pore throats and chambers, this intuition is mostly valid for the gas phase but not for the liquid. This can be understood from the flow regimes appearing in pores of different sizes at different saturations. We refer the reader to our earlier paper [1] for images of the various flow morphologies cited below. At low saturation levels S < 20%, 44  4.3. Results  (a)  (b) Figure 4.5: Variation of the relative permeability of (a) the wetting (gas) phase and (b) the non-wetting (liquid) phase with the respective capillary number. All other parameters are at the baseline values given for Fig. 4.3.  the drop flow regime prevails in almost all pores. Thus, increasing S leads to larger drops that flow in the central part of the pore, being lubricated by a gas layer that separates the drops from the solid walls. As a result, the liquid permeability Krn increases and the gas permeability K rw decreases with S. In this range, the behavior is consistent with convention in the literature [18]. Beyond S = 20%, the liquid drops become larger and start to come  45  4.3. Results  Figure 4.6: Relative permeabilities as functions of the liquid (non-wetting) saturation S. F = 1.4 and all other parameters are at their baseline values given for Fig. 4.3. into contact with the solid walls, especially in the smaller pores. This significantly suppresses the lubrication effect. Some of the liquid may even been trapped in the pore chamber and excluded from the liquid throughput altogether. The dominant flow regime becomes bubble-slug flow. Consequently, Krn declines with S. Meanwhile, in the smallest pores the gas starts to be entrapped in the pore chambers by liquid. Thus K rw continues its decline. For S > 40%, gas flux is maintained only in the largest pores, in the form of gas bubbles carried by a continuous liquid stream. In this stage, the liquid permeability levels off and becomes insensitive to S. Finally, as S → 1, the gas bubbles shrink and eventually disappear; K rn climbs back to unity. The distinctive features of Kri seen in our simulations are closely related to the geometry of our model porous medium. The interfacial morphology and flow regimes have much to do with interface pinning at sharp corners and with initial layout of the interface. Such detailed characteristics have not been recorded in prior experiments and computations. It is therefore particularly interesting to compare our results with prior studies. There is a wealth of experimental and numerical data in the literature on K ri as functions of S, as the saturation is the parameter most thoroughly studied. Unfortunately, there are great variations in the values of other parameters,  46  4.3. Results  (a)  (b) Figure 4.7: Comparison of the relative permeabilities Kri (S) with previous studies. (a) Krw of the wetting phase; (b) Krn of the non-wetting phase. Open symbols denote computational results while closed ones experimental results. For our result, the body force F = 1.4 and the other parameters are the same as in Fig. 4.3.  which make quantitative comparison very difficult. We have collected in Fig. 4.7 what appear to be reliable and representative data sets for a comparison with our numerical results. All studies have produced a Krw in close agreement with one another. In general, the gas permeability Krw decreases with S, which in our convention 47  4.3. Results denotes the saturation of the liquid, non-wetting phase. This robust trend can be attributed to the gradual narrowing of flow areas available to the gas with increasing S. Our Krw being relatively low may be because in our geometry, the gas tends to be trapped in the pore chamber in certain flow regimes. For Krn , there are much greater variations in the data. Previous results seem to segregate into two groups: the tube-bundle [23] and pore-network models [79] predict a higher Krn that exceeds unity, while lattice-Boltzmann computations in random porous media [44, 48] and experiments [2, 4] give a much lower Krn that increases monotonically with S. Our prediction falls into the first group for smaller S, and then crosses over to the second for larger S. In the first group, the larger K rn is clearly due to lubrication. In relatively simple geometries, stratification of the two phases occurs more consistently than in tortuous and random flow conduits. Besides, the lubrication effect is accentuated by a high non-wetting-to-wetting viscosity ratio M , which happens to be the case in the first group of studies (M ∼ 10 compared with M ∼ 1 for the second group). These conspire to produce the higher Krn . With increasing S, the water tends to touch the solid wall at the pore throat as mentioned above. This hampers lubrication and causes our Krn to decrease and join the second group of curves. In general, one can summarize the Krn data as follows. The non-wetting phase typically stays away from the solid walls, and its interfacial morphology is more sensitive to material and flow parameters than for the wetting phase. Hence, its relative permeability exhibits greater variation among different studies.  4.3.3  Wettability effect  Wettability of the pores is an important determinant of the relative permeabilities. This property affects the interfacial morphology and, if a threephase contact line appears, how fast the contact line moves on the walls. To examine this effect in our pore-scale calculations, we have used the same baseline setup (Fig. 4.1) and varied the contact angle θ from 60 ◦ to 165◦ . Figure 4.8 depicts the variation of K ri with θ. The main finding is that by making the pores more hydrophobic, with increasing θ, both relative permeabilities tend to decrease up to θ = 135 ◦ . With further increase in θ, the liquid relative permeability K rn increases while that for gas Krw stays near zero. Note that for θ ≤ 90◦ , the liquid becomes the wetting phase. For simplicity, however, we have not switched the subscript in Fig. 4.8; K rn refers to the liquid and Krw to the gas regardless of θ. If the solid is hydrophilic (θ < 90◦ ), the smaller pores are blocked by a 48  4.3. Results  Figure 4.8: Relative permeabilities as functions of the liquid wetting angle θ on the solid surface. F = 1.4, and the other parameters are at the baseline values of Fig. 4.3. The insets show dominant flow patterns among the pores responsible for the overall Kri .  liquid meniscus, but the dominant flow regime in larger pores is the liquid drop regime depicted by insets in the figure. Some liquid is retained in the pore chamber, while the rest is carried by the gas as large drops in the center of the pore. The liquid phase enjoys a large K rn because of this lubrication effect. Its decline from θ = 60 ◦ to 135◦ has to do with the narrow pore throat in our geometry. With increasing θ, more liquid is initially driven out of the pore chamber. It then makes contact with the wall in the throat and thus suppresses lubrication by the underlying gas. As a result, Krn decreases. While hydrophobicity rejects the liquid from the near-wall regions, the gas is increasingly pushed into the pore chambers. As a result, the gas permeability Krw continues to decline. At θ = 135◦ , the liquid tends to seal the gas completely in the chamber for smaller pores (lower inset) or to form liquid slugs separated by thin gas films in the larger pores (upper inset). The gas permeability approaches zero. At even higher θ, the liquid-solid contact becomes so energetically prohibitive that the liquid starts to detach from the wall, allowing a gas film to  49  4.3. Results  (a)  (b) Figure 4.9: Effect of the contact angle θ on the relative permeabilities K rw (a) and Krn (b): comparison with lattice-Boltzmann computations in the literature. wedge in between. In the end, the shell flow regime prevails in most of the pores, as illustrated by the inset for θ = 165 ◦ , with a liquid core enveloped by a thin gas sheath. This leads to the observed recovery of K rn , again thanks to the lubrication effect. Figure 4.9 compares our prediction of the wettability effect with several 50  4.3. Results lattice-Boltzmann computations. We have found no clear-cut experimental results; those that changed θ by using different materials typically had other parameters changed as well. Among the lattice-Boltzmann computations, the geometry ranges from relatively regular (pore-network) to highly random. The physical parameters differ widely as well. For K rw , the qualitative trend is the same for all studies: Krw decreases with θ. Note, however, that the prior computations are limited to relatively large θ, where the decrease of Krw with θ is mild. Our computation includes smaller θ where the effect is stronger. The difference in the numerical values among different studies has probably stemmed from the different physical and geometric parameters. For Krn , the picture is somewhat murkier. Li et al. [48] showed K rn to increase with θ, Huang and Lu [37] and Hao and Cheng [34] predicted a mild decrease, while Nguyen et al. [58] showed a non-monotonic trend. Our Krn first decreases and then increases with θ, and the magnitude is also much higher than the others. This can be attributed to our geometric and physical parameters amplifying the lubrication effect. For example, among the studies that have given the viscosity ratio M —M = 1 in [48], 12 in [37] and 18 in our study—Krn increases with M . This is due to the smaller friction on the non-wetting phase when the wetting phase viscosity is reduced (more on this in the next subsection). Besides, our geometry is regular and lacks tortuosity and cross-linkage among pores. It is well-known that lubrication effect is dampened in more complex and disordered pore geometry [79].  4.3.4  Viscosity ratio effect  The viscosity ratio affects the relative permeabilities mainly through lubrication effect. This has been noted in experiments [4] and computations in relatively simple geometries [23, 79]. We have used our baseline setup of Fig. 4.1 and computed a range of the viscosity ratio, from M = 0.2 to 18. Recall that M is the viscosity ratio between the non-wetting (liquid) and the wetting (gas) phase. Figure 4.10 shows that K rn increases monotonically with M while Krw decreases, both leveling off for large M . As explained below, this result confirms the key role played by lubrication effect. Since the relative permeability for each phase is computed using its own viscosity (Eq. 1.1), it is convenient to think of the M effect as due to changing the viscosity of the other phase. For instance, the increase of the liquid permeability Krn with M can be viewed as due to the gas viscosity getting lower. As Kri is computed by averaging over potentially different flow regimes in pores of different sizes, it is not easy to rationalize the variations 51  4.3. Results  Figure 4.10: Relative permeabilities as functions of the liquid-gas viscosity ratio. F = 1.4 and all other parameters are at the baseline values as given for Fig. 4.3. of Kri in terms of the flow patterns. However, a qualitative argument can be made by looking at the flow regime in the dominant pore sizes in the PSD, roughly between 15 and 20 µm. Figure 4.11 depicts the flow pattern in a pore of radius Rt = 18.25 µm for increasing M values. For smaller M values, shell flow prevails. Reducing the wetting-phase viscosity enhances the lubrication of the non-wetting fluid in the center. Hence K rn increases with M . The same effect is at work in the annular-droplet regime for larger M , with the liquid above the gas pocket being lubricated, although quantitatively it is not as pronounced. Similarly, in analyzing the gas-phase relative permeability K rw , we imagine that M is being increased by raising the liquid viscosity while the gas viscosity remains constant. As M increases from 0.2 to 1, Fig. 4.11 indicates a transition from the shell-flow to the bubble-slug regime. This traps gas in the pore chamber and reduces Krw sharply; all gas transport now occurs in the form of bubbles. With further increase in M , not only does the gas bubbles become smaller, but its velocity also decreases with increasing liquid viscosity. This explains the decline of K rw with M . The effect of the viscosity ratio on the relative permeability is an issue of some controversy in the literature [18]. Prior experiments and computations  52  4.3. Results  Figure 4.11: Flow patterns in a pore of radius Rt = 18.25 µm for increasing M values. Shell flow prevails for M = 0.2 while bubble-slug flow obtains for larger M . All other parameters are at the baseline values as given for Fig. 4.3. have been done at widely different values for the contact angle, saturation and capillary number. Nevertheless, we have compiled some data in Fig. 4.12 for a qualitative comparison. All computations, including ours, show a common trend that Krw decreases with M while Krn increases. This is mostly due to lubrication, although the effect should be tempered in real porous medium by the disordered nature of pore geometry [79]. The two experiments do not agree with the computational trend. [2] showed both K ri to decrease with M while [4] showed both to increase with M . We do not have an explanation for this discrepancy. Clearly more computational and experimental data are needed for a coherent picture to emerge.  4.3.5  Effects of pore geometry and initial configuration  Aside from the factors considered so far, [1] have demonstrated that the interface morphology and flow regimes are also sensitive to the pore geometry and the initial configuration. For one, sharp corners tend to pin interfaces. This gives rise to very strong hysteresis in flow regimes, i.e. dependence on initial condition and flow history. Consequently, the relative permeability is expected to depend on these as well [4]. Ahmadlouydarab et al. [1] have examined these effects in considerable detail by varying the initial configurations and approaching a given flow condition by ramping up or down the flow rate. In the following we will briefly explore their ramifications for the relative permeabilities. Figure 4.13 compares flow in micro-pores with sharp and rounded corners (a and b), and in the latter case, two different initial positions for the liquid  53  4.3. Results  (a)  (b) Figure 4.12: Comparison of (a) Krw and (b) Krn as functions of the viscosity ratio M with previous studies. Open symbols denote computations and filled ones experiments.  bridge (b and c). In Fig. 4.13(a), interface pinning at the two sharp corners slows down the movement of the liquid core and allows time for all the liquid to be drawn out of the pore chamber. In the end the liquid flow regime emerges, with the gas completed trapped in the pore chamber. In comparison, the flow develops much more rapidly with the rounded corners 54  4.3. Results  (a)  (b)  (c) Figure 4.13: Effect of pore geometry and initial interfacial configuration. (a) Pore with sharp corners and the liquid bridge initially at the upstream end of the pore chamber. (b) Pore with rounded corners and the liquid bridge initially at the upstream end of the pore chamber. (c) Pore with rounded corners and the liquid bridge initially at the downstream end of the pore chamber. The parameters are the same for all three cases: effective pore radius Re = 20 µm, α = 1, β = 3, γ = 2, S = 40.5%, FR = 1.4, θ = 135◦ and M = 18. (Fig. 4.13b), leading to the drop flow regime. To evaluate Kri for pores with sharp and rounded corners, we again average the flow of both phases over the PSD of Fig. 4.2. Obviously other flow regimes will appear in pores of larger and smaller sizes than that of Fig. 4.13, but the general trend seems to hold: sharp corners pin interfaces and generally tend to reduce both Kri relative to the rounded geometry. For the geometric and flow parameters tested, the sharp-cornered pores give Krw = 0.023 and Krn = 0.42, while the round-cornered ones have Krw = 0.069 and Krn = 1.34. In the latter, the large relative permeability for the liquid phase is thanks to the lubrication of the liquid drops by the gas in the larger pores (Fig. 4.13b). If the liquid bridge is initially next to the downstream end of the pore chamber (Fig. 4.13c), it wraps around the solid ridge and forms a continuous liquid film sealing some gas in the chamber. Meanwhile gas bubbles are transported in the middle in this bubble-slug regime. The difference between Fig. 4.13(b) and (c) is a special case of the hysteretic effect, where flow history affects the interfacial pattern that emerges eventually [1]. After 55  4.4. Conclusion averaging over the PSD, we obtain Krw = 0.013 and Krn = 0.53, comparable to the sharp-cornered geometry and much below those of Fig. 4.13(b). As in Fig. 4.13(a ), the liquid flow is hindered by wall friction and much of the gas is trapped. The lack of lubrication effect is the key factor in the lower Kri .  4.4  Conclusion  In this work we used pore-scale flow simulations to compute the relative permeability of a model porous medium and test the validity of the socalled extended Darcy’s law for gas-liquid two-phase flow. We used similar baseline model as our first study. The main conclusions can be summarized as follows. (a) The averaged gas and liquid flow rates both depend nonlinearly on the imposed pressure gradient. Thus the extended Darcy’s law does not hold in the usual sense. Instead, it can be viewed as the definition of the relative permeabilities, which are complex functions of the geometric and flow parameters of the system. (b) The relative permeability of the wetting phase K rw behaves more or less as expected from previous studies. It increases with the wettingphase capillary number and the wetting-phase saturation, but decreases with the wettability of the pore surface and with the nonwetting-to-wetting viscosity ratio. (c) The relative permeability of the non-wetting phase K rn has a much more complex behavior. It increases monotonically with the viscosity ratio, but exhibits non-monotonic dependence on the saturation and wetting angle. This is because the non-wetting phase tends to occupy the central part of the pore, and Krn is much more sensitive to the inteface morphology than Krw . In particular, lubrication of the nonwetting phase by a cushion of the wetting phase is responsible for elevated Krw , and loss of lubrication, e.g. when the non-wetting phase makes contact with the wall, leads to lower K rw . The strength of our work lies in the pore-scale resolution of the interfacial shape and its temporal evolution, as well as contact-line motion on the solid walls of the pore. Thus, the relative permeabilities can be directly related to the flow regimes and interfacial evolution in the two-phase flow. The 56  4.4. Conclusion role of interface pinning at sharp corners is highlighted, which leads to a dependence on flow history. The shortcoming of this work is the simplicity of its pore geometry. Periodic contraction (pore throat) and expansion (pore chamber) are accounted for, but connectivity among pores is neglected, as is the disordered nature of pore geometry in real porous media. This may have magnified the lubrication effect. Comparisons with experimental and computational results in the literature show qualitative agreement in cases where a consistent trend exists. In other cases, previous papers have presented conflicting data. Our study adds more data to this situation, but cannot resolve the existing discrepancies. More studies under carefully controlled conditions are needed in the future.  57  Chapter 5  Motion and coalescence of sessile drops 5.1  Introduction  This part of my research has been motivated by efforts to optimize the airwater transport through the gas diffusion layer (GDL) in proton exchange membrane (PEM) fuel cells. This is a porous medium through which oxygen flows to the reaction site while the water, a product of the reaction, is discharged [57, 75]. Water vapor condenses in the hydrophobic pores, and small droplets are believed to coalesce and then be driven out partly by wetting gradients. Meanwhile, the incoming air creates a counterflow, which tends to drive the water droplets in the opposite direction [32]. How do sessile drops move and coalesce on a substrate bearing a wetting gradient while being subject to an external flow? This is the question that we set out to investigate. Considerable work has been done on drop motion driven by wetting gradients. Early theoretical models such as lubrication approximation, wedge approximation theory have been tried to demonstrate the relevant physics by balancing the driving force —the net capillary force due to the asymmetry in contact angles—against the viscous friction to predict the velocity of drop motion [10, 30]. These models typically neglect drop deformation and assume the shape of a spherical cap for the drop, its footprint being a circle of radius R. The models differ mostly in the treatment of the contact-line singularity and determination of the viscous friction. One can estimate the friction locally from Tanner’s law [10], or calculate it using Cox’s wedge solution or a lubrication solution [68], and in the end, the difference amounts to a different prefactor. Recently, Xu and Qian[77] have reported detailed numerical solutions of the problem using a phase-field model. Although this is much more sophisticated than prior theoretical models, with the drop deformation properly calculated and the contact-line singularity regularized by Cahn-Hilliard diffusion, their results are qualitatively the same.  58  5.2. Problem setup and methodology Experimentally, the idea of moving droplets by wetting gradients was demonstrated by Chaudhury and Whiteside [11]. Based on this brief summary of theoretical, numerical and experimental work, the motion of a single droplet on a substrate having a wetting gradient is well understood. In contrast, little can be found in the literature on drop-drop interaction on wetting gradients. The most relevant studies dealt with coalescence between a drop moving on a wetting gradient and a stationary one resting on an adjoining substrate of uniform wettability [46, 76]. Drop motion driven by an external flow has been studied for some time [19, 20, 22, 43, 53, 65, 88], but only on a substrate with uniform wettability. Besides, to our knowledge, there have been no studies of drop coalescence driven by external flow. Here, this study tries to fill these gaps in knowledge. We carry out systematic 2D simulations of drop motion and coalescence on a substrate subject to wetting gradients and an external flow. For a single drop and a pair of drops, we consider three scenarios in which the drop motion is driven (a) by a wetting gradient, (b) by an external flow of the surrounding fluid parallel to the substrate, and (c) by a combination of the two. In particular, we investigate the competition between the two driving forces, and how they determine the coalescence between the two drops. For the simple case of a single drop on a wetting gradient, we also report a 3D simulation that helps to benchmark the numerical parameters used in 2D simulations.  5.2  Problem setup and methodology  We examine two related phenomena, the motion of a single sessile drop and the collision and coalescence of a pair of sessile drops. Most of the simulations to be reported are in a 2D planar geometry, with a rectangular computational domain of length L and height H (Fig. 5.1). From the area of the drop Ad , we define its effective diameter De = 8Ad /π, which will be used as the characteristic length throughout the paper. If the contact angle varies along the substrate, no initial drop shape will be at equilibrium. As soon as the drop makes contact with the substrate, it adjusts to the local contact angle, deforms and starts to move. For consistency, we have used an initial drop shape of a semicircle in most of the simulations. The only exception is Fig. 5.8, which examines the effect of the initial drop shape. Three domain sizes have been used. For single drop simulations, α = L/De = 6.67, β = H/De = 0.834. The two-drop simulations have a somewhat larger domain to accommodate the merged drop: α = 7.88, β = 0.985. The initial drop separation γ = D c /De will be varied for individual 59  5.2. Problem setup and methodology  Figure 5.1: Schematic of the initial configuration for 2D planar computations. The drops have an effective diameter De , and the geometry is specified by the length ratios α = L/De , β = H/De , and for two drops, γ = Dc /De .  runs. Numerical experiments show that the domain is long enough so α has negligible effect on the drop motion. The domain height β is chosen to be order one so as to model a small pore; its value affects drop motion driven by an external flow, but not that driven by wetting gradients on the substrate. Finally, to compare with experiments in §3.1, we have adopted the experimental geometry with α = 8.11 and β = 0.845. The problems to be computed involve the deformation and movement of the interface, a three-phase contact line, and non-uniform wettability. For droplets of sub-millimeter diameters, the Bond number is much below unity. Thus we have neglected gravity. Similarly, for the operating conditions that motivated this work, the Reynolds number is typically Re = O(10 −2 ) [44], and we have neglected inertia for most of the simulations. We impose stress-free boundary conditions on the top, the left and the right sides of the domain. On the solid substrate, no-slip conditions are used for velocity, and contact line motion is realized via Cahn-Hilliard diffusion [84]. As a consequence, the “contact-line velocity” is extracted not from the velocity field U , but from the movement of the point of contact, i.e. the point where the interface φ = 0 intersects the substrate. We have no-penetration condition for φ, and prescribe the contact angle via surface energies between each fluid and the solid [80, 84]. To simulate an external flow driven by a horizontal pressure gradient, we impose the body force B on both fluids. Thus a Poiseuille velocity profile develops at the entry and exit of the domain as they are sufficiently far from the drops. In presenting results, we use the effective drop diameter D e as the characteristic length, and the capillary time tc = µd De /σ as the characteristic time, µd being the drop viscosity. The dimensionless parameters of the problem include the length ratios indicated in Fig. 5.1, the drop-to-medium viscosity ratio  60  5.3. Results: the motion of a single sessile drop M = µd /µm and the density ratio ρd /ρm . As mentioned above, most of the simulations are inertialess; the density ratio is relevant only to one set of inertial results in §4.2. Besides, the wetting gradient W e = d(cos θ)/dx is made dimensionless as G = De W e, and the maximum contact angle at the hydrophobic end of the substrate is specified as θ m . When external flow is present, we have an effective Bond number F = BD e2 /σ that indicates the strength of the external flow relative to the surface tension. Beside these, the Cahn-Hilliard model introduces two mesoscopic dimensionless parameters: the Cahn number Cn = /D e and the diffusion parameter Λ = ld /De . As explained in Chapter 2, these parameters must be chosen judiciously. Cn should be small enough so the sharp-interface limit is approached [85, 90]. For 2D drop motion on a substrate, we have found that reducing Cn from 2 × 10−2 to 6.7 × 10−3 causes at most a 0.5% difference in the drop position throughout the simulation. Hence all subsequent 2D results are computed using Cn = 6.7 × 10 −3 . On the other hand as mentioned in previous chapters, l d is the counterpart of the slip length l s commonly used in sharp-interface models [80, 84], Λ represents the strength of Cahn-Hilliard diffusion in moving the contact line, and is closely related to the contact line speed. Thus, it should in principle be determined by fitting an experimental datum for the specific fluids and substrate material [80]. We discuss the choice of this parameter next in §3.1.  5.3 5.3.1  Results: the motion of a single sessile drop Single drop driven by wettability gradient  This subsection has two objectives. The first is to use the prior theoretical and experimental results [54, 68] to determine the parameter Λ for our simulations, and the second is to examine the motion of a single sessile drop driven by a uniform wetting gradient. Moumen et al. [54] conducted an experiment on the motion of sessile drops on a substrate bearing a prescribed profile of the wetting angle, which is reproduced in Fig. 5.2 along with the dimensionless wetting gradient G. They compared their experimental results with the prediction of slip-based theoretical models of [68] applied to the wetting gradient of Fig. 5.2. We will simulate the drop motion under the same wetting gradient. As we discussed in Chapter 2, our diffusion length l d is related to the slip length ls by ls ≈ 2.5ld . This is the basis for evaluating our l d and Λ. Figure 5.3 compares our 2D simulations using two values of Λ with the model predictions of [68] based on Cox’s wedge-flow approximation. Inertia 61  5.3. Results: the motion of a single sessile drop  Figure 5.2: Experimental conditions of Moumen et al. [54] study used to benchmark our simulation. (a) Contact angle profile along the solid substrate. (b) Profile of the dimensionless wettability gradient G.  Figure 5.3: Comparison of the droplet center velocity, non-dimensionalized into a capillary number, between our 2D computations and the wedge-flow model predictions of [68]. For each Λ value, a slip length ls = 2.5ld = 2.5ΛDe is used in the theoretical model. is neglected in all four calculations. The drop center x c is the midpoint between the two contact lines, scaled by the effective drop diameter D e . The instantaneous drop velocity u, taken to be the average between the velocity of the two contact lines, is non-dimensionalized into a capillary number Ca = µd u/σ. Our diffusion lengths are matched to their slip lengths by l s = 2.5 ld . Thus our computations at Λ = 2.14 × 10 −3 and 8.33 × 10−3 correspond, respectively, to slip length of ls = 5.35 × 10−3 De and 2.83 × 10−2 De . The drop velocity roughly mirrors the profile of the wetting gradient G 62  5.3. Results: the motion of a single sessile drop  Figure 5.4: Comparison of the drop center velocity among our 3D computation, experimental measurement and theoretical predictions of the wedge-flow model [54]. Our computation is based on Cn = 2.03 × 10−2 and Λ = 6.76 × 10−3 , and one of the wedge-flow predictions is at the corresponding slip length. A much smaller ls /De = O(10−7 ) is needed to predict the experimental data.  (Fig. 5.2b); the drop accelerates first, and then decelerates as it spreads out on the increasingly hydrophilic substrate. As expected, the drop moves faster at the larger Λ value or slip length. The slip-model predictions are somewhat lower than our numerical results. Given the various simplifications and approximations in the model, e.g., the assumption of a spherical cap and the wedge-flow approximation, this quantitative discrepancy is not surprising. This comparison confirms the correspondence between l d and ls and offers a basis on which to choose our numerical diffusion length. It remains to be determined whether the Λ values tested above are realistic. As noted before [25, 80], the slip length for real materials are usually in the nanometer range, much below the typical resolution of continuum computations. Thus, previous computations using sharp-interface and diffuse-interface models alike have adopted computable slip lengths ls /De = O(10−2 ) [22, 65, 77], comparable to those in Fig. 5.3. To illustrate this limitation, we have done a 3D simulation and compared it with the experimental result of [54] (Fig. 5.4). The 3D computational domain has a width of W = 2De , and the drop is located in the center plane. The front half of the box is computed with symmetry conditions imposed on the center plane that cuts the drop in half. At Λ = 6.76 × 10 −3 , close to the minimum 63  5.3. Results: the motion of a single sessile drop  Figure 5.5: Instantaneous velocity of the sessile drop on a substrate with a constant wettability gradient for four values of the viscosity ratio M = 1, 10, 50 and 100. G = 0.144, θm = 92.3◦ . diffusion length that we can readily compute in 3D, our phase-field simulation overpredicts the drop velocity roughly by a factor of 5. The wedge-flow model overpredicts by a smaller amount at the corresponding slip length. [54] also tested the model at a much smaller slip length, l s = 0.5 nm, corresponding to Λ = 1.35 × 10−7 . Now the prediction is within a factor of 2 of the measurement. Thus, it is clear that we will have to settle for a larger, computable Λ value as others have done. To achieve quantitative agreement with experiments, Yue and Feng [80] introduced wall-energy relaxation in the Cahn-Hilliard model to compensate for using an unrealistically large diffusion length, and recommended choosing the relaxation parameter by fitting an experimental datum. In this study we will not strive to match a specific experiment, and thus have used Λ = 8.33 × 10 −3 for all the 2D simulations presented hereafter. This needs to be kept in mind when interpreting the numerical results quantitatively. We now turn to the more general question of the behavior of a sessile drop on a constant wetting gradient G. As the drop experiences viscous friction on the substrate as well as on its upper surface, its motion turns out to be sensitive to the viscosity of both the drop and the surrounding fluid. Figure 5.5 shows the evolution of the instantaneous drop velocity for four different M values as the drop traverses the substrate with a fixed 64  5.3. Results: the motion of a single sessile drop G. Ca exhibits a gentle acceleration for M = 1 but a deceleration for large M . As inertia has been neglected, the drop velocity is determined by the balance between two forces: the driving force due to the differing contact angles at the leading and trailing edges, and the hydrodynamic drag. The latter consists of an internal friction between the drop and the substrate, and an external drag due to the surrounding medium. As the drop moves downstream into more hydrophilic areas, the driving force increases with the lengthening footprint of the drop. The change in the drag, on the other hand, depends on M . For M = 1, the medium is as viscous as the drop and the external drag is significant. As the drop spreads, the internal friction increases not only because the footprint of the drop lengthens, but also because the decreasing height induces greater shear rates inside. In the meantime, the external drag declines with the height of the spreading drop. The balance between the driving force and the total drag is such that the drop sees a slight acceleration in this case. For larger M , achieved by reducing the outer viscosity, the external drag becomes insignificant and the drop moves faster. Now the increasing internal friction overwhelms the increasing driving force to produce the gentle deceleration. Note also that the effect of M diminishes for larger M ; the medium viscosity becomes negligible. For all subsequent results, we have used a fixed M = 50 to approximate the viscosity ratio between water and air at room temperature. For such an M , the external viscous friction can be neglected in analyzing the results. To explore the effect of the wetting gradient G, we track the average velocity of the drop from the inception of its motion at x c = 2.08 till xc = 4.17, and plot the capillary number based on this average velocity as a function of the wetting gradient G (Fig. 5.6). As expected, the drop velocity increases with G. The driving force, due to the different contact angles at the two contact lines, is proportional to G. As the drop moves downstream, however, its changing shape affects the viscous friction in a nonlinear fashion. Thus, the Ca ∼ G relationship is not strictly linear. This is supported by the wedge-flow approximation, whose prediction is also plotted for comparison.  5.3.2  Single drop driven by external flow  The motion of a sessile drop driven by an external flow has been studied by several groups before, for example [65] and [20], and the physics is more or less well understood. The objective of this brief subsection is to reprise the basic features of the motion as they are relevant to analyzing more complex scenarios later in the paper. 65  5.3. Results: the motion of a single sessile drop  Figure 5.6: The average velocity of the drop measured between x c = 2.08 to xc = 4.17 as a function of the dimensionless wettability gradient G. Our computation is compared with the wedge-flow model of [54].  Consider a drop deposited on an ideal substrate with a uniform contact angle θ. In the absence of contact angle hysteresis, such a drop moves with the slightest external driving force. We impose a constant body force on both the drop and the surrounding fluid; this is equivalent to applying a fixed pressure drop across the length of our computational domain. A shear flow develops and sweeps over the drop, which deforms and starts to move. In roughly 10tc both the external flow and the drop motion approach a steady state. Our geometric setup of Fig. 5.1 is such that a parabolic velocity profile develops at the entry and the exit. Because of the additional drag incurred by the sessile drop, the steady-state centerline velocity is lower than that expected in a Poiseuille flow, which is u m = BH 2 /(2µm ) in dimensional form, and the deficit depends on the drop size and the contact angle. The steady-state velocity of the drop, in terms of the capillary number Ca = µd u/σ, is plotted in Fig. 5.7 against the dimensionless driving force F for two contact angles, θ = 45◦ and 135◦ . Note that an “external capillary number” can be defined using the characteristic velocity u m of the external flow, the medium viscosity µm and the effective drop diameter De : Cam = µm um /σ = 12 F β 2 , which ranges from 0 to 5.22 × 10−2 for the F values in Fig. 5.7. In our inertialess flow, the drop velocity is determined by a balance between the driving force on the drop surface due to the external flow and  66  5.3. Results: the motion of a single sessile drop  Figure 5.7: Dependence of the steady-state drop velocity, in terms of the capillary number Ca, on the driving force F for two contact angles. The insets indicate the degree of drop deformation at the largest F value in each case.  the viscous friction on the substrate. The deformation of the drop, which is appreciable in the parameter range of Fig. 5.7 (see inset), should have introduced nonlinearity into the problem. Yet surprisingly, Ca scales more or less linearly with the driving force F . Similar linear dependence has been reported in the literature [20, 65]. Furthermore, the drop moves faster on a more hydrophobic substrate; it has a smaller footprint and is also taller, and thus experiences a larger driving force and a smaller friction. At higher flow rates and with inertia, the drop may become entrained by the flow and pinch off from the wall [43, 65], but such regimes are not directly relevant to the current study.  5.3.3  Competition between external flow and wetting gradient  If the external flow and the wetting gradient push the drop in the same direction, the two effects reinforce each other and the outcome is faster drop motion. If the two oppose each other, their competition can be quite subtle; the drop dynamics exhibits a sensitivity to initial conditions and forcing history, and is therefore hysteretic. This can perhaps be anticipated from the fact that both the capillary force and the external force depend intimately on the shape of the drop. For example, the drop’s footprint along 67  5.3. Results: the motion of a single sessile drop the flow direction determines the magnitude of the capillary force, and its hump determines how much hydrodynamic drag it receives from the external flow. In addition, the viscous friction depends on the shape of the drop as well. Thus, depending on the initial drop configuration and the flow history, the outcome of the competition may vary. We have probed this hysteresis by numerical experiments. Figure 5.8 compares the trajectory of three drops starting from the center of the domain (initial drop center position at x0 = 3.33) with different initial shapes, the instantaneous contact angle being 45 ◦ , 90◦ and 135◦ . The initially tallest drop experiences the greatest drag from the external flow, and starts to move to the left. As the contact lines adjust to the local wetting angles, however, the drop spreads out. Consequently the hydrodynamic drag on the drop declines while the wetting force increases. The two eventually reaches a balance when the drop comes to rest at x r = 3.25, having moved about 8.4% of its effective diameter D e . Similarly, the initially most spreadout drop favors the wetting gradient and moves to the right (final position xr = 3.39), and the intermediate drop winds up at an intermediate position (xr = 3.31). If F and G are grossly mismatched, the drop may find no stationary state at all. This scenario will be examined shortly. We have also tested dependence on the history of the control parameters, and confirmed that depending on how F varies in time, the drop can be driven to different locations and shapes even though the final F is the same. This hysteretic behavior, i.e. dependence on initial configuration and history of the control parameters, is rooted in the fact that it takes a finite time—on the order of the capillary time t c — for the drop to adjust its shape. It is a fundamental feature of the F ∼ G competition in the current flow setup. Perhaps unfortunately, this implies that the quantitative results of the simulations, if not the qualitative trends, are specific to the initial conditions used. In the following, we will use a consistent initial condition to investigate the F ∼ G competition, with the drop being a semicircle centered at x0 = 3.33, the substrate bearing a constant G toward the right, and the external flow at a fixed F being turned on at t = 0 in the opposite direction. When comparing such numerical results with experiments, one needs to beware of the initial conditions in the latter, and to keep the hysteretic effects in mind. Starting from this initial condition, we simulate the motion of the drop at different combinations of (F, G) values. The final resting position of the drop center, xr , is plotted in Fig. 5.9 as functions of F for a series of G values. This figure paints a rather delicate picture about the F ∼ G competition. For each G value, there is a narrow range of F within which a 68  5.3. Results: the motion of a single sessile drop  Figure 5.8: Trajectory of drops starting from three initial shapes with instantaneous contact angles of 45◦ , 90◦ and 135◦. The drop is initially centered at x0 = 3.33. The constant wetting gradient, with G = 0.245 and θm = 130◦, drives the drop toward the right, while the external flow at F = 6.12 × 10−2 goes to the left. The insets illustrate the initial and final shapes of the drop. stationary state can be achieved. Within this range, the drop stops further to the hydrophilic area (larger xr ) for larger G and smaller F , and to the hydrophobic area (smaller xr ) for smaller G and larger F . The stationary drop assumes such a shape that the hydrodynamic and capillary forces balance each other. It is a fragile balance. For an F value outside the narrow range, the hydrodynamic force will never be balanced by the capillary force. The drop moves continually to the left or right, and no state of rest is achieved. It is somewhat surprising that the drop attains a stationary state at all. One may imagine shifting the drop slightly to the left by external perturbation. Then the substrate becomes more hydrophobic and the drop is supposed to bulge upward. This would reduce the capillary force and increase the hydrodynamic drag, and the drop would then continue to move to the left. This argument seems to suggest a linear instability of any resting position. What makes the rest state possible is again the hysteretic nature of the system. After the imagined slight shift to the left, the drop shape does not immediately adapt to the local contact angle; that process takes a finite time on the order of the capillary time t c . The ensuing change in drop shape 69  5.4. Results: drop-drop coalescence  Figure 5.9: Final position for a single drop that has come to rest under the opposing actions of F and G. The y-axis shows the displacement of the drop center scaled by the effective drop diameter. Initially the drop is centered at x0 = 3.33.  and hence the hydrodynamic drag cannot be anticipated from the simple “equilibrium” argument posed above. In fact, there is no equilibrium shape to relax to on a wetting gradient. We have carried out further numerical experimentation on the delicate balance between F and G. Taking one of the curves in Fig. 5.9 corresponding to a fixed G, we vary F in small increments, each time after a new stationary position has been established. This way, we have confirmed that the stationary state can be extended to a much wider F range than indicated in Fig. 5.9, where all stationary states have been reached from the same initial position under a fixed F .  5.4 5.4.1  Results: drop-drop coalescence Coalescence driven by wettability gradient  Merging of two sessile drops driven by a constant wetting gradient on the substrate is a relatively simple affair. Figure 5.5 has shown that in a lowviscosity medium (M ≥ 10), a single drop will decelerate as it moves into increasingly hydrophilic regions, thanks to the increasing viscous friction on the substrate. When two identical drops are placed on such a substrate 70  5.4. Results: drop-drop coalescence  Figure 5.10: Coalescence of two identical drops on a substrate with a constant wetting gradient G = 0.244 and maximum contact angle θm = 164◦ . The drops are initially centered at x = 3.30 and 4.44 (γ = 1.14), and the two solid curves indicate the trajectories of the front contact line of the trailing drop and the rear contact line of the leading drop. After merging at t = 30, the center of the new drop moves according to the dotted line. with an initial separation, the one initially on the more hydrophobic area will catch up with the other, and coalesce with it. Figure 5.10 demonstrates such a process. In fact, this merging can hardly be called “interaction” since it arises mostly from the behavior of single drops at different positions on the substrate. Since the surrounding fluid has low viscosity, the two drops move more or less independently of each other. As one may expect, the time for catch-up increases with initial separation and decreases with the wetting gradient G.  5.4.2  Coalescence driven by external flow  Now we place two drops on a substrate of uniform wetting angle θ = 105 ◦ , at an initial separation of γ = 1.02. After turning on an external flow driven by F = 0.309, the two drops both start to move downstream (Fig. 5.11). Based on the characteristic velocity of the external flow u m , the medium viscosity µm and the effective drop diameter De , the external capillary number for each drop is Cam = 0.150, large enough to produce considerable drop 71  5.4. Results: drop-drop coalescence  Figure 5.11: Coalescence of two identical drops on a substrate with uniform contact angle θ = 105◦ driven by an external flow at F = 0.309 toward the right. The drops are initially centered at x = 1.46 and 2.48 (γ = 1.02), and the curves indicate drop positions in the same way as in Fig. 5.10. Coalescence starts at t = 164. To magnify the view of the curves near the point of coalescence, we have omitted the initial period until t = 150. deformation. In time, the trailing drop catches up with the leader, and the two merge into one. The merged drop moves with a steady-state velocity u that corresponds to a capillary number Ca = µ d u/σ = 0.147. Note that the catch-up takes rather a long time, and the drops have moved through the length of the domain a few times. To accommodate the long simulation, we have adopted periodic boundary conditions at the two ends of the domain. The coalescence implies that, despite the absence of inertia, the trailing drop experiences a greater hydrodynamic drag than the leading one. To investigate this, we have analyzed the velocity, pressure and stress distributions around the two drops. As the drops start to move, both deform with the top tilting downstream. This creates a relatively narrow gap between the drops, where the velocity of the surrounding fluid is suppressed by viscosity. This is demonstrated in Fig. 5.12(a ) by comparing profiles of ux , the x component of the velocity, around the two drops. At the contact lines, ux exhibits very large gradients within an effective “slip layer” produced by Cahn-Hilliard diffusion [84]. But the feature of interest here is  72  5.4. Results: drop-drop coalescence  (a)  (b) Figure 5.12: Comparison of (a) the x component of the velocity and (b) the pressure distribution on two drops at a dimensionless time t = 123. The velocity and pressure values are taken along the contour of φ = 0.9, just outside the interface, and are plotted against the polar angle θ measured from the midpoint between the two contact lines. Velocity and pressure are made dimensionless by σ/µd and σ/De , respectively. that the leading drop experiences a lower velocity on its upstream side than the trailing drop. This can be thought of as a viscous shielding effect due to the asymmetric shape of the trailing drop. In Stokes flow around an object with fore-aft symmetry, no such shielding is expected. A direct consequence of the shielding is a lower pressure on the upstream side of the leading drop 73  5.4. Results: drop-drop coalescence (Fig. 5.12b). Note that the pressure also varies sharply within the slip layer. Since this layer contributes little to the overall hydrodynamic forcing of the drops, we have matched the two pressure profiles at a small distance of 4.16 off the substrate at the upstream contact line. This distance is roughly the interfacial thickness [81], and corresponds to a polar angle of 3.4◦ off the substrate. The pressure being set to zero at this point, the rest of the profile shows clearly that the shielded drop experiences much reduced dynamic pressure on its upstream side. This explains the greater driving force on the trailing drop, and the eventual catching-up and merging of the two. We have also analyzed the viscous shear and normal stresses around the drops. The pressure contributes roughly 85% of the total drag, and is therefore the dominant factor. At the instant analyzed in Fig. 5.12, the pressure force on the trailing drop is 35% greater than that on the leading one, and the total drag is 32% greater. Numerical experiments show that if the two drops are initially farther apart, the shielding effect will be weaker, and the flow-driven coalescence will take longer. In principle, two drops will always interact hydrodynamically because of the elliptic nature of the Stokes flow. Recall that our study was motivated by drop motion in hydrogen fuel cells, and the typical operating conditions correspond to negligible inertia. But since inertia is known to affect hydrodynamic interaction among particles and drops [24, 59], we carried out a set of simulations to explore the effect of finite Reynolds numbers. At Re = ρ m um D/µm = 13, with a dropto-medium density ratio of 100 and a viscosity ratio of 50, the merging of the two drops happens much more quickly. Under the same conditions of Fig. 5.11, the two drops first touch at a dimensionless time of 17.0, down from 164 without inertia. Evidently, the shielding effect is much enhanced by inertia. For the rest of the results to be presented, inertia will again be omitted.  5.4.3  Competition between external flow and wetting gradient  Let us consider the interaction of two sessile drops when a wetting gradient G and an external flow F are both at work. As in §3.3, we expect the drop dynamics to exhibit hysteresis, i.e. sensitivity to initial conditions and forcing history. For simplicity and brevity, we will ignore hysteresis in this subsection, and use the same initial drop shape (semi-circle) and constant F and G forcing in the following. The situation is simpler if F and G are co-current, i.e. oriented in the 74  5.4. Results: drop-drop coalescence  Figure 5.13: Drop coalescence under co-current F and G: the catch-up time t and catch-up distance d as functions of F for a fixed G = 0.183 and θm = 136◦ . Initially the drops are centered at x0 = 3.30 and 4.44, with initial separation γ = 1.14. The catch-up distance d is defined as the distance travelled by the center of the trailing drop up to the time of coalescence.  same direction. Both F and G promote catch-up in this case, and it is no surprise that the catch-up time decreases monotonically when either driving force is increased. Figure 5.13 demonstrates this effect for an increasing F . Note that F exerts a larger force on the trailing drop not only because of viscous shielding as discussed in Fig. 5.12, but also because the trailing drop sits in a more hydrophobic area, and thus presents a taller hump to the external flow. On the other hand, the catch-up distance d, defined as the distance travelled by the center of the trailing drop up to the instant of coalescence, shows a non-monotonic dependence on F (Fig. 5.13). A larger F produces faster motion on both drops. For the G value used, increasing F beyond roughly 0.05 will cause the catch-up distance to increase as the gain in drop velocity more than compensates for the loss in catch-up time. This non-monotonicity may have been present for the F -only scenario of Fig. 5.11 already. Since coalescence happens slowly when driven solely by F , we have not carried out the numerous long runs required to confirm this. We have also tested the effect of varying G at a fixed F . For increasing G, both the catch-time and the catch-up distance decrease monotonically. This is because the leading drop spreads out more rapidly downstream, and thus decelerates more quickly. 75  5.4. Results: drop-drop coalescence  Figure 5.14: Drop coalescence under opposing effects of the external flow (toward the left) and the wetting gradient (toward the right). G = 0.244, θm = 164◦ , and F = 0.01. The initial configuration is identical to that of Fig. 5.10, and the curves have the same meanings. Coalescence occurs at t = 46.  Now we consider the counter-current situation with the external flow toward the left and the wetting gradient driving the drop toward the right. In this case, the presence of the wetting gradient fundamentally modifies the shielding effect, and the F ∼ G competition cannot be anticipated from the effect of each alone. Keeping G fixed and increasing F gradually, we have predicted two distinct scenarios, with coalescence or separation of the drops. Coalescence occurs for sufficiently small F such that F weakens but does not overwhelm the G-driven catch-up and coalescence discussed in §4.1. An example is shown in Fig. 5.14 for G = 0.244, θ m = 164◦ and F = 0.01. Both drops are driven by G to move toward the right. Since the leading drop sits on a more hydrophilic surface, it spreads out more and consequently experiences less hydrodynamic drag from the external flow. Thus, the shielding effect is reversed; the trailing drop (on the left) receives a greater drag due to the external flow, which hinders its motion and delays the catch-up and coalescence. The catch-up time and distance both increase monotonically with F in this scenario (Fig. 5.15). Since F exerts a larger drag on the left drop, a sufficiently large F will prevent coalescence altogether. This is the second scenario, drop separation. 76  5.5. Conclusion  Figure 5.15: The catch-up time t and catch-up distance d both increase with F , other conditions being the same as in Fig. 5.14.  Depending on the relative strength of F and G, each of the two drops may continue to move to the right, come to rest, or reverse course and move to the left, according to the single-drop dynamics of Fig. 5.9. In principle, therefore, one may envision several sub-scenarios of non-coalescence. Figure 5.16 illustrates two of those. At the moderately large F = 0.051 (Fig. 5.16a), the drag on the left drop overpowers the G-induced driving force. It stops, and then reverses direction and moves slowly to the left. In the meantime, the right drop continues to move to the right, and the two separate. For the even larger F = 0.154 (Fig. 5.16b), the hydrodynamic drag on both drops is able to overcome the wetting force. Thus, both drop are driven toward the left. Since the right drop spreads out more, it receives a smaller driving force from F , a greater resistance from G, as well as a larger viscous friction on the substrate. Thus it lags behind the left drop, and the two separate in time.  5.5  Conclusion  Using a diffuse-interface formalism, we have carried out a numerical investigation of the motion and coalescence of sessile drops driven by a wettability gradient, an external flow along the substrate, or both. Even with one of the two driving agents at work, the drops exhibit interesting and subtle dynamics. When the two are set against each other, their competition can 77  5.5. Conclusion  Figure 5.16: Drop separation under opposing effects of the external flow (toward the left) and the wetting gradient (toward the right). G = 0.244, θm = 164◦ . (a) At F = 0.051, the right drop continues to the right but the left one moves to the left. (b) At F = 0.154, both drops move to the left while separating from each other. produce subtle and sometimes unexpected effects on the dynamics of sessile drops, including hysteresis and either coalescence or separation of the drops. The new findings of the study, within the parameter ranges investigated, are summarized as follows: (a) A single sessile drop on a constant wetting gradient decelerates in its motion provided that the viscosity of the surrounding medium is much lower than that of the drop. This is owing to the increasing viscous friction on the drop as it moves into more hydrophilic areas. (b) When the wetting gradient is opposed by an external flow, the motion of a single drop exhibits a strong dependence on the initial condition and the history of forcing. The hysteresis arises because both the wetting force and the hydrodynamic drag depend strongly on the drop shape, and the drop shape adapts to contact angles and external flows not instantaneously but on a finite time scale. (c) Two identical drops driven by a constant wetting gradient will coalesce into one. This is a direct result of the deceleration of each drop, with the leading one suffering a greater reduction in speed. (d ) Two identical drops driven by an external flow will coalesce, thanks to a viscous shielding effect associated with the asymmetric shape of the upstream drop. This occurs even at zero Reynolds number, but is greatly enhanced by inertia. 78  5.5. Conclusion (e) When the wetting gradient is opposed by an external flow, two drops may coalesce or separate, depending on the relative strength of the two forcing agents. Coalescence is favored by a stronger wettability gradient, and separation by a stronger external flow. These results may have practical implications for microfluidic devices and micro-fabrication processes that rely on the manipulation of droplets. For example, the sensitivity to initial conditions and forcing history may lead to new strategies of controlling drop motion, coalescence and deposition. By balancing the effects of external flow and substrate wettability, one may design gas diffusion layers for optimal air-water transport in hydrogen fuel cells, the engineering application that has served as our initial motivation. In the meantime, we point out several limitations of the present work that should be kept in mind when interpreting the numerical results quantitatively. First, most of the simulations are in two dimensions (2D). Comparisons with our own limited 3D simulations and with 3D experiments show qualitative agreement but quantitative differences. One may well expect, for instance, the viscous shielding effect to be attenuated in 3D relative to its 2D counterpart. Second, we do not yet have a satisfactory physical model for the moving contact line, and its numerical simulation invariably requires additional ad hoc input. In our case, it is the Cahn-Hilliard diffusivity, which produces a diffusion length that is the counterpart of the slip length in sharp-interface models. Finally, we have not strived for a comprehensive description of the hysteretic effect. After testing a few scenarios to demonstrate the sensitivity of drop motion to initial conditions and forcing history, we have limited subsequent simulations to a uniform initial condition for ease of analysis. Potentially there may exist other, more exotic, regimes of drop dynamics that can be realized through different initial conditions and forcing history. Such scenarios remain to be explored in future work.  79  Chapter 6  Conclusion and recommendations The overarching theme of the research conducted in this thesis is to use numerical simulations to explore pore-scale hydrodynamics. This is a necessary and significant step in studying two-phase flows in porous media. Up to now, most of our knowledge has been obtained from observations and measurements of overall quantities on the macroscopic scale. While valuable as guidelines for engineering designs and developments, such knowledge has lacked a firm scientific foundation, and thus has been susceptible to misinterpretations and incorrect generalization. For example, the so-called extended Darcy’s law has been widely used, despite contravening evidence. Various microscopic flow regimes and scenarios have been hypothesized. But high resolution imaging has been difficult in small dimensions. Therefore, numerical simulations fill a much needed role in this context. Several theoretical development and computational innovations have made the current research possible. The use of a diffuse-interface theory not only allows straightforward interface-capturing, but also provides a physically reasonable phenomenology for the interfacial tension. The latter idea, combined with the Cahn-Hilliard diffusion and transport equation, further serves to regularize the contact line singularity. And the contact line turns out to play important roles in defining the pore-scale hydrodynamics. Numerically, the computation has been facilitated by the use of a flexible adaptive mesh generator, without which high-resolution computation would be prohibitively expensive. In view of the outcome of the three projects, one can claim that the pore-scale simulations have succeeded in identifying and clarifying the local flow mechanisms underlying macroscopic observations. In addition, we have uncovered novel physics in the transition between flow regimes and hysteretic behavior in sessile drop dynamics. In the following, I will first summarize the key findings, and then reflect on their significance and limitations.  80  6.1. Summary of key findings  6.1 6.1.1  Summary of key findings Interfacial flow regimes in corrugated microchannels  To thoroughly understand the local hydrodynamics of two-phase flow in micropores, we made the problem tractable by simplifying the geometry into axisymmetric tubes with periodic constrictions. We computed interfacial deformation, breakup and coalescence from hydrodynamic principles, and delineate various flow regimes based on the control parameters. Within the parameter ranges explored in this work, depending on the body force, capillary force and the level of liquid saturation, a number of flow regimes were seen in the corrugated microchannel: gas flow, blockage, liquid flow, bubble-slug flow, droplet flow, annular flow and annular-droplet flow. We constructed a 2D map of flow regimes for a set of geometric and flow parameters starting from a prescribed initial configuration. Some of the regimes such as gas flow, blockage, liquid flow and shell flow regimes are new, and there is significant hysteresis when transitions happen on the flow map. Different flow patterns can result from different initial conditions under identical control parameters. Considering the differences in geometric setup, flow-control schemes and parameter ranges, the numerical simulations are in remarkably close agreement with experiments in terms of the boundaries between flow regimes and the holdup ratio. Perhaps the most important findings of this study are the new flow regimes and the hysteresis, which highlight the importance of geometric features of the flow conduit, e.g. corners in micropore where the interface tends to pin.  6.1.2  Relative permeability of interfacial flows in a model porous media  Complex flow pathways, viscous coupling and transient flows are factors that challenge the validity of the extended Darcy’s law, which assumes a linear relationship between the velocity and pressure gradient for each phases. Accordingly, the relative permeabilities K ri should be functions of several parameters in addition to the saturation level. Because of this multi-variate dependence, experimental data show large quantitative scatters as well as qualitative contradictions in the literature. An obvious interpretation of this discrepancy is that two-phase transport in porous media depends strongly on the interfacial morphology and fluid dynamics near the interface. Using an accurate numerical simulation to demonstrate the connection between the  81  6.2. Significance and limitations relative permeability and the interfacial dynamics is the main contribution of this project. In this study we used pore-scale flow simulations to compute the relative permeability of a corrugated tube bundle model as a porous medium. Results have shown a generally nonlinear relationship between the pressure gradient and the velocity of each fluid phase. Lubrication of the non-wetting phase by the wetting phase is the main reason for the experimentally observed dependence of the relative permeability on the viscosity ratio. Besides, the simultaneous appearance of different flow regimes in pores of various sizes causes non-monotonic dependence of the relative permeability of the non-wetting phase on the saturation level. Variation in the hydrophobicity of the pore walls is another factor. Additionally, different initial configuration and/or pore geometry may lead to different relative permeabilities for each fluid, which demonstrates a dependence of the relative permeability on the flow history. All these factors contribute to viscous coupling, and invalidate the extended Darcy’s law.  6.1.3  Sessile drops motion and coalescence  Our motivation for studying sessile drop behavior under the simultaneous action of a wettability gradient and an external flow comes from air-water transport in hydrogen fuel cells, although similar dynamics takes place in many other natural and engineering processes, including water-oil flow in enhanced oil recovery. A salient feature of the drop dynamics is the hysteresis in the final resting position of drops when they move on solid substrates bearing wetting gradient and exposed to an external flow. Beside the competition among the three forces: external flow force, wetting gradient force, and viscose friction, the initial drop configuration is another important factor in determining the drop shape and motion. To have a successful coalescence for a pair of drops, the external flow should be co-current with the drop motion induced by the wetting gradient. In the the case of the two forces being counter-current to each other, the external flow should be weak compared with the wetting gradient force. If the external flow is strong enough, no coalescence will be seen; the drops will separate from each other, with one or both moving against the wetting gradient.  6.2  Significance and limitations  The insights gained from this research are potentially useful in two general ways. First, the novel phenomena that we have simulated may inspire 82  6.2. Significance and limitations further, more in-depth research. Although we have uncovered certain interfacial dynamics, we have not at all achieved a complete understanding of the physics involved. This is due in part to the limitations of our methodology (more on this below), and also to the finite scope of a Ph.D. thesis. Nevertheless, the findings point to new directions of research which deserve to be explored in depth. Second, a thorough understanding of interfacial flow characteristics in micropores is a fundamental step towards the rapid development of microfluidic and lab-on-chip devices. For example, optimization of the watermanagement strategy in micropores as well as in gas channel of PEM fuel cells promises to bring forth significant improvement in current hydrogen fuel cell technology. Based on our results, a gas diffusion medium with an engineered pattern of wetting gradient from the catalyst layer towards the gas channels would help to remove liquid water more efficiently. Current GDL technology uses PTFE coating to make the pores more hydrophobic. To build a wetting gradient into the GDL requires a more refined coating procedure or a gradient in pore size. In the meantime, to coat the pores with a large amount of PTFE may have adverse effects on water removal because of pore size reduction. Another potential application of our research is in guiding the miniaturization of microelectronic components that incurs high-flux heat generation from small areas. For such applications, the development of microscale heat exchangers using two-phase flows is a promising technology. The limitations of this research can be summarized as follows: 1. Two dimensionality. A real porous media is a complex 3D structure, but software and hardware capability have limited us to mostly 2D computations (planar or axisymmetric). For one, in simulating sessile drop dynamics, almost all of our computations were in planar twodimensional domains. In this way we failed to include effects of the surrounding fluid that flows around the sides of the drops. 2. Simplicity of the pore geometry. In calculating the interfacial flow regimes and relative permeabilities, we used a corrugated microchannel as a representative geometry. This does not include some potentially important geometric features of real pores, including branching and connection of micropores. This would call for a truly three-dimensional computation. 3. Limitations of the theoretical model. Currently, there are no satisfactory, self-contained contact line model, and our Cahn-Hilliard model 83  6.2. Significance and limitations has to be supplemented by ad hoc input. Besides, we cannot compute realistic, nanoscale, diffusion (or slip) lengths but have to settle for much larger lengths. A side effect of Cahn-Hilliard diffusion is drop shrinkage, specially when a small drop is surrounded by a large domain of the other fluid. This requires careful choice of the model parameters to avoid undesirable diffusion and shrinkage effects while capturing the key physics of the problem. 4. Numerical method. We have used a fully implicit numerical scheme for stability, but the inversion of large matrices turns out to be very time consuming, and cannot be easily parallelized. Thus generalization to 3D computation has been difficult and limited in scope. There have been promising developments in algorithms that may mitigate this limitation in the future. 5. Limited range of initial setup. Although hysteresis has been shown to be an important effect, we generally avoided an exhaustive test of the dependence of the dynamics on initial conditions, mostly owing to the finite scope of the thesis. Besides, some simplifications may have introduced unrealistic conditions into the physics of the problem. For example, in the relative permeability study (Chapter 4), we assumed the same local saturation for all different pore sizes of the model porous media, but in reality capillary forces would produce different level of saturations in different pores. This summary makes it clear that there are myriad new problems that remain to be investigated in two-phase flows in micropores. It is our hope that the findings of this thesis will inspire other researchers to contribute to this area of research, perhaps by bringing new tools and strategies to the as yet unresolved questions.  84  Bibliography [1] M. Ahmadlouydarab, Z. S Liu, and J. J. Feng. Interfacial flows in corrugated microchannels: Flow regims, transitions and hysteresis. Int. J. Multiphase Flow, 37:1266–1276, 2011. [2] J. O. Amaefule and L. L. Handy. The effect of interfacial tensions on relative oil/water permeabilities of consolidated porous media. Society of Petroleum Eng. J., 22:371–381, 1982. [3] Javier Atencia and David J. Beebe. Controlled microfluidic interfaces. Nature, 437:648–655, 2005. [4] D. G. Avraam and A. C. Payatakes. Flow regimes and relative permeabilities during steady-state two-phase flow in porous media. J. of Fluid Mech., 293:207–236, 1995. [5] D. G. Avraam and A. C. Payatakes. Generalized relative permeability coefficients during steady-state two-phase flow in porous media, and correlation with the flow mechanisms. Transport in Porous Media, 20:135–168, 1995. [6] D. G. Avraam and A. C. Payatakes. Flow mechanics, relative permeabilities, and coupling effects in steady-state two-phase flow through porous media, the case of strong wettability. Ind. Eng. Chem. Res., 38:778–786, 1999. [7] A. M. Barajas and R.L. Panton. The effects of the contact angle on twophase flow in small diameter pipes. Int. J. Multiphase Flow, 19:337–346, 1993. [8] D. Barnea, O. Shoham, Y. Taitel, and A. E. Dukler. Flow pattern transition for gas-liquid flow in horizontal and inclined pipes: comparison of experimental data with theory. Int. J. Multiphase Flow, 6:217–225, 1980.  85  Bibliography [9] M. J. Blunt, M. D. Jackson, M. Piri, and P. H. Valvatne. Detailed physics, predictive capabilities and macroscopic consequences for porenetwork models of multiphase flow. Adv. Water Resour., 25:1069–1089, 2002. [10] F. Brochard. Motions of droplets on the solid surfaces induced by chemical or thermal gradients. Langmuir, 5:432–438, 1989. [11] M. K. Chaudhury and G. M. Whitesides. How to make water run uphill. Science, 256:1539–1541, 1992. [12] Thomas Cubaud and Chih-Ming Ho. Transport of bubbles in square microchannels. Phys. Fluids, 16:4575–4585, 2004. [13] Thomas Cubaud, Umberto Ulmanella, and Chih-Ming Ho. Two-phase flow in microchannels with surface modifications. Fluid Dyn. Res., 38:772–786, 2006. [14] C. A. Damianides and J. W. Westwater. Two-phase patterns in a compact heat exchanger and in small tubes. In Proc. 2nd UK National Conference on Heat Transfer, pages 1257–1268, London, 1988. Institute of Mechanical Engineers. [15] E. Dana and F. Skoczylas. Experimental study of two-phase flow in three sandstones. I measuring relative permeabilities during two-phase steady-state experiments. Int. J. Multiphase Flow, 293:1719–1736, 2002. [16] S. Daniel, S. Sircar, J. Gliem, and M. K. Chaudhury. Ratcheting motion of liquid drops on gradient surfaces. Langmuir, 20:4085–4092, 2004. [17] S. R. A de Loos, R. M. Tiggerlaar J. van der Schaaf, T. A Nijhuis, M. H. J. M. de Croon, and J. C. Schouten. Gas-liquid dynamics at low Reynolds numbers in pillared rectangular micro channels. Microfluid. Nanofluid., 9:131–144, 2010. [18] A. H. Demond and P. V. Roberts. An examination of relative permeability relations for two-phase flow in porous media. Water Resour. Bull., 23:617–628, 1987. [19] P. Dimitrakopoulos and J. J. L. Higdon. On the displacement of threedimensional fluid droplets adhering to a plane wall in viscous pressuredriven flows. J. Fluid Mech., 435:327–350, 2001. 86  Bibliography [20] Hang Ding, Mohammad N. H. Gilani, and Peter D. M. Spelt. Sliding, pinch-off and detachment of a droplet on a wall in shear flow. J. Fluid Mech., 644:217–244, 2010. [21] N. Djilali. Computational modelling of polymer electrolyte membrane (PEM) fuel cells: Challenges and opportunities. Energy, 32:269–280, 2007. [22] Jean-Baptiste Dupont and Dominique Legendre. Numerical simulation of static and sliding drop with contact angle hysteresis. J. Comput. Phys., 229:2453 – 2478, 2010. [23] R. Ehrlich. Viscous coupling in two-phase flow in porous media and its effect on relative permeabilities. Transport in Porous Media, 11:201– 218, 1993. [24] J. Feng, H. H. Hu, and D. D. Joseph. Direct simulation of initial value problems for the motion of solid bodies in a newtonian fluid. Part 1. Sedimentation. J. Fluid Mech., 261:95–134, 1994. [25] M. Fermigier and P. Jenffer. An experimental investigation of the dynamic contact angle in liquid-liquid systems. J. Colloid Interface Sci., 146:226–241, 1991. [26] L. A. Freitag and C. F. Ollivier-Gooch. Tetrahedral mesh improvement using swapping and smoothing. Int. J. Numer. Methods Eng., 40:3979– 4002, 1997. [27] T. Fukano and A. Kariyasaki. Characteristics of gas-liquid two-phase flow in a capillary tube. Nucl. Eng. Des., 141:59–68, 1993. [28] P. Gao and J. J. Feng. Enhanced slip on a patterned substrate due to depinning of contact line. Phys. Fluids, 21:102102, 2009. [29] G. R. Gerauld and S. J. Salter. The effects of pore-structure on hysteresis in relative permeability and capillary pressure. Transport in Porous Media, 5:103–151, 1990. [30] H. P. Greenspan. On the motion of a small viscous droplet that wets a surface. J. Fluid Mech., 84:125–143, 1978. [31] A. K. Gunstensen and D. H. Rothman. Lattice-boltzmann studies of immiscible two-phase flow through porous media. J. of Geophysical Research, 98:6431–6441, 1993. 87  Bibliography [32] V. Gurau and J. A. Mann, Jr. A critical overview of computational fluid dynamics multiphase models for proton exchange membrane fuel cells. SIAM J. Appl. Math., 70:410–454, 2009. [33] A. K. Gustensen and D. H. Rothman. Lattice-boltzmann studies of immiscible two-phase flow through porous media. J. Geophys. Res., 98:6431–6441, 1993. [34] L. Hao and P. Cheng. Pore-scale simulations on relative permeabilities of porous media by lattice boltzmann method. Int. J. Heat and Mass Transfer, 53:1908–1913, 2010. [35] I. Hassan, M. Vaillancourt, and K. Pehlivan. Two-phase flow regime transitions in microchannels: A comprehensive experimental study. Nanoscale Microscale Thermophys. Eng., 9:165–182, 2005. [36] H. H. Hu, N. A. Patankar, and M. Y. Zhu. Direct numerical simulations of fluid-solid systems using the arbitrary Lagrangian-Eulerian technique. J. Comput. Phys., 169:427–462, 2001. [37] H. Huang and X. Lu. Relative permeabilities and coupling effects in steady-state gas-liquid flow in porous media: A lattice Boltzmann study. Phys. Fluids, 21:092104, 2009. [38] Y. Ito, M. Heydari, A. hashimoto, T. Konno, A. Hirasawa, S. Hori, K. Kurita, and A. Nakajima. The movement of a water droplet on a gradient surface prepared by photodegradiation. Langmuir, 23:1845– 1850, 2007. [39] D. Jacqmin. Contact-line dynamics of a diffuse fluid interface. J. of Fluid Mech., 402:57–88, 2000. [40] M. Joanicot and A. Ajdari. Droplet control for microfluidics. Science, 309:887–888, 2005. [41] D. D. Joseph and Y. Y. Renardy. Fundamentals of Two-Fluid Dynamics. Part II: Lubricated Transport, Drops and Miscible Liquids. Springer-Verlag, New York, 1993. [42] F. Kalaydjian. Origin and quantification of coupling between relative permeabilities for two-phase flows in porous media. Transport in Porous Media, 5:215–229, 1990.  88  Bibliography [43] QJ Kang, DX Zhang, and SY Chen. Displacement of a threedimensional immiscible droplet in a duct. J. Fluid Mech., 545:41–66, 2005. [44] T. Koido, T. Furusawa, and K. Moriyama. An approach to modeling two-phase transport in the gas diffusion layer of a proton exchange membrane fuel cell. J. Power Sources, 175:127–136, 2008. [45] E. C. Kumbur, K. V. Sharp, and M. M. Mench. On the effectiveness of Leverett approach for describing the water transport in fuel cell diffusion media. J. of Power Sources, 168:356–368, 2007. [46] Yu-Hsuan Lai, Miao-Hsing Hsu, and Jing-Tang Yang. Enhanced mixing of droplets during coalescence on a surface with a wettability gradient. Lab Chip, 10:3149–3156, 2010. [47] R. Lenormand, E. Touboul, and C. Zarcone. Numerical models and experiments on immiscible displacements in porous media. J. of Fluid Mech., 189:165–187, 1988. [48] H. Li, C. Pan, and C. T. Miller. Pore-scale investigation of viscous coupling effects for two-phase flow in porous media. Phys. Rev. E, 72:026705, 2005. [49] P. L. Maffettone and M. Minale. Equation of change for ellipsoidal drops in viscous flow. J. Non-Newtonian Fluid Mech., 78:227–241, 1998. [50] J. M. Mandhane, G. A. Gregory, and K. Aziz. A flow pattern map for gas liquid flow in horizontal pipes. Int. J. Multiphase Flow, 1:537–553, 1974. [51] C. Marle. Multiphase flow in porous media. Editions Technip, Paris, 1981. [52] K. Mishima, T. Hibiki, and H. Nishihara. Some characteristics of gas-liquid flow in narrow rectangular ducts. Int. J. Multiphase Flow, 19:115–124, 1993. [53] B. M. Mognetti, H. Kusumaatmaja, and J. M. Yeomans. Drop dynamics on hydrophobic and superhydrophobic surfaces. FARADAY DISCUSSIONS, 146:153–165, 2010. [54] N. Moumen, R. S. Subramanian, and J. B. McLaughlin. Experiments on the motion of drops on a horizontal solid surface due to a wettability gradient. Langmuir, 22:2682–2690, 2006. 89  Bibliography [55] M. Muskat and M. W. Meres. The flow of heterogeneous fluids through porous media. Physics, 7:346–363, 1936. [56] J. H. Nam and M. Kaviany. Effective diffusivity and water-saturation distribution in single- and two-layer PEMFC diffusion medium. Int. J. Heat and Mass Transfer, 46:4595–4611, 2003. [57] J.-H. Nam, K.-J. Lee, G.-S. Hwang, C.-J. Kim, and M. Kaviany. Microporous layer for water morphology control in PEMFC. Int. J. Heat and Mass Transfer, 52:2779–2791, 2009. [58] V. H. Nguyen, A. P. Sheppard, M. A. Knackstedt, and W. V. Pinczewski. The effect of displacement rate on imbibition relative permeability and residual saturation. J. Pet. Sci. Eng., 52:54–70, 2006. [59] P. O. Olapade, R. K. Singh, and K. Sarkar. Pairwise interactions between deformable drops in free shear at finite inertia. Phys. Fluids, 21:063302, 2009. [60] J. F. Oliver, C. Huh, and S. G. Mason. Resistance to spreading of liquids by sharp edges. J. Colloid Interface Sci., 59:568–581, 1977. [61] C. F. Ollivier-Gooch and C. Boivin. Guaranteed-quality simplical mesh generation with cell size and grading control. Engineering with Computers, 17:269–286, 2001. [62] W. Rose. Petroleum reservoir engineering at the crossroads (ways of thinking and doing). The Iran Pet. Inst. Bull., 46:23–27, 1972. [63] W. Rose. Second thoughts on darcy’s law. The Iran Pet. Inst. Bull., 48:25–30, 1974. [64] J. Rupert. A Delaunay refinement algorithm for quality 2-dimensional mesh generation. J. Algorithm, 18:548–585, 1995. [65] AD Schleizer and RT Bonnecaze. Displacement of a two-dimensional immiscible droplet adhering to a wall in shear and pressure-driven flows. J. Fluid Mech., 383:29–54, 1999. [66] V. P. Schulz, J. Becker, A. Wiegmann, P. P. Mukherjee, and C. Y. Wang. Modeling of two-phase behavior in the gas diffusion medium of PEFCs via full morphology approach. J. Electrochem. Soc., 154:B419– B426, 2007. 90  Bibliography [67] E. J. Soares, M. S. Carvalho, and P. R. S. Mendes. Liquid water removal from a polymer electrolyte fuel cell. J. Fluids Eng., 127:24–31, 2005. [68] R. S. Subramanian, N. Moumen, and J. B. McLaughlin. Motion of a drop on a solid surface due to a wettability gradient. Langmuir, 21:11844–11849, 2005. [69] Qian T., Wang X. P., and Sheng P. Molecular hydrodynamics of the moving contact line in two-phase immiscible flows. Comm. Comput. Phys., 1:1–52, 2006. [70] Y. Taitel, D. Barnea, and A.E. Dukler. Modeling flow pattern transitions for steady upward gas-liquid flow in vertical tubes. AIChE J., 26:345–354, 1980. [71] Y. Taitel and A. E. Dukler. A model for predicting flow regime transitions in horizontal and near horizontal gas-liquid flow. AIChE J., 22:47–55, 1976. [72] K. A. Triplett, S.M. Ghiaasiaan, S.I. Abdel-Khalik, A. LeMouel, and B. N. McCord. Gas-liquid two-phase flow regimes in microchannels. Part II: void fraction and pressure drop. Int. J. Multiphase Flow, 25:395–410, 1999. [73] Cahn J. W. Critical-point wetting. J. Comput. Phys., 66:3667, 1977. [74] Cahn J. W. and Hilliard J.E. Free energy of a nonuniform system.I. interfacial free energy. J. Comput. Phys., 28:258, 1958. [75] C.-Y. Wang. Fundamental models for fuel cell engineering. Chem. Rev., 104:4727–4766, 2004. [76] H. Wang, Q. Liao, X. Zhu, J. Li, and X. Tian. Experimental studies of liquid droplet coalescence on the gradient surface. J. Supercond. Nov. Magn., 23:1165–1168, 2010. [77] X. Xu and T. Qian. Droplet motion in one-component fluids on solid substrates with wettability gradients. Phys. Rev. E, 85:051601, 2012. [78] T. Ye, W. Shyy, and J. N. Chung. A fixed-grid, sharp-interface method for bubble dynamics and phase change. J. Comput. Phys., 174:781–815, 2001.  91  Bibliography [79] A. G. Yiotis, J. Psihogios, M. E. Kainourgiakis, A. Papaioannou, and A. K. Stubos. A lattice Boltzmann study of viscous coupling effects in immiscible two-phase flow in porous media. Colloids and Surfaces A: Physicochem. Eng. Aspects, 300:35–49, 2007. [80] P. Yue and J. J. Feng. Wall energy relaxation in the Cahn-Hilliard model for moving contact lines. Phys. Fluids, 23:012106, 2011. [81] P. Yue, J. J. Feng, C. Liu, and J. Shen. A diffuse-interface method for simulating two-phase flows of complex fluids. J. Fluid Mech., 515:293– 317, 2004. [82] P. Yue, C. Zhou, and J. J. Feng. A computational study of the coalescence between a drop and an interface in Newtonian and viscoelastic fluids. Phys. Fluids, 18:102102, 2006. [83] P. Yue, C. Zhou, and J. J. Feng. Spontaneous shrinkage of drops and mass conservation in phase-field simulations. J. Comput. Phys., 223:1– 9, 2007. [84] P. Yue, C. Zhou, and J. J. Feng. Sharp-interface limit of the CahnHilliard model for moving contact lines. J. Fluid Mech., 645:279–294, 2010. [85] P. Yue, C. Zhou, J. J. Feng, C. F. Ollivier-Gooch, and H. H. Hu. Phasefield simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing. J. Comput. Phys., 219:47–67, 2006. [86] P. Yue, C. Zhou, and J.J. Feng. Phase-field simulation of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing. J. Comput. Phys., 219:47–67, 2006. [87] F. Y. Zhang, X. G. Yang, and C. Y. Wang. Liquid water removal from a polymer electrolyte fuel cell. J. Electrochem. Soc., 153:A225–A232, 2006. [88] J. Zhang, M. J. Miksis, and S. G. Bankoff. Nonlinear dynamics of a two-dimensional viscous drop under shear flow. Phys. Fluids, 18:072106, 2006. [89] C. Zhou, P. Yue, and J. J. Feng. The rise of Newtonian drops in a nematic liquid crystal. J. Fluid Mech., 593:385–404, 2007.  92  Bibliography [90] C. Zhou, P. Yue, J. J. Feng, C. F. Ollivier-Gooch, and H. H. Hu. 3D phase-field simulations of interfacial dynamics in Newtonian and viscoelastic fluids. J. Comput. Phys., 229:498–511, 2010.  93  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items