RIVULET FLOW A N D STABILITY by MARGARET HONGXIA LIANG B A S c . (Engineering Physics) U n i v e r s i t y of B r i t i s h C o l u m b i a , 2000 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of M a t h e m a t i c s Institute of A p p l i e d M a t h e m a t i c s We accept this thesis as. conforming to the required stanclar/i THE UNIVERSITY OF BRITISH COLUMBIA December 2002 (c) Margaret H o n g x i a L i a n g , 2002 In presenting this degree at the thesis in partial fulfilment of University of British Columbia, I agree freely available for reference copying of department this or publication of and study. his or her representatives. may be It is requirements that the for an Department The University of British Columbia Vancouver, Canada advanced Library shall make it that permission for extensive granted by the understood head that this thesis for financial gain shall not be allowed without permission. DE-6 (2/88) I further agree thesis for scholarly purposes by the of my copying or my written Abstract T w o dimensional, steady-state rivulet flow d o w n an i n c l i n e d plane, a n d combined r i v u l e t / a i r flow i n a circular pipe are studied. B o t h asymptotic analysis a n d n u m e r i c a l methods are used to solve the Navier-Stoke equations for rivulet flow d o w n an i n c l i n e d plane. A s y m p t o t i c approxi m a t i o n , w h i c h is v a l i d only for rivulet of large cross-sectional aspect ratios, is substantiated by exact n u m e r i c a l solutions up to specified error tolerance. S t a n d a r d shooting techniques are used to solve O D E for rivulet shape and finite element methods are used to approximate the flow velocity. T h e stability of the rivulet flow is also studied based on energy criterion. W e find that pure gravity driven rivulets are subject to break up a n d pure shear driven rivulets are stable i n that sense. For the combined a i r / r i v u l e t flow, the relationship between pressure gradients a n d fluxes is investigated using the annulus rivulet a n d circular arc r i v u l e t models. T h e annulus rivulet m o d e l is solved analytically, a n d the circular arc r i v u l e t is c o m p u t e d using F E M L A B . It is found that pressure gradient needed to drive the r i v u l e t flow is m u c h more sensitive to change i n channel size t h a n to change i n contact angle. F i n a l l y , linear stability analysis of the combined a i r / r i v u l e t flow i n a rectangular d o m a i n is formulated. n Table of Contents Abstract ii Table of Contents iii List of Figures iv List of Tables vi Acknowledgement Chapter 1. 1.1 1.2 1.3 vii Introduction 1 Background A p p l i c a t i o n to fuel cells M a t h e m a t i c a l formulation - Chapter 2. Rivulet Flow Down an Inclined Plane - Asymptotic Analysis Chapter 3. Rivulet Flow On A n Inclined Plane - Numerical M e t h o d 3.1 3.2 Numerical method Stability Analysis Based O n Energy Criterion Chapter 4. 1 2 4 7 15 15 24 Combined A i r and Rivulet Flow 28 4.1 A i r Only Channel Flow 4.2 Annulus Rivulet 30 4.3 Circular A r c Rivulet . . : 32 4.4 4.3.1 T h e direct p r o b l e m 4.3.2 T h e inverse p r o b l e m Pressure gradient vs. flux of water 33 34 38 Chapter 5. : Stability Analysis of Combined A i r and Rivulet Flow 30 40 5.1 S t a b i l i t y of the annulus m o d e l 40 5.2 L i n e a r stability of a rectangular geometry 40 Chapter 6. Conclusions and Future Work Appendix .1 48 50 Verification of femlab c o m p u t a t i o n using a rectangular d o m a i n Bibliography 50 51 iii List of Tables 3 . 1 N u m e r i c a l estimates of Q c and L c for gravity d r i v e n rivulet 4.1Pressure gradient A ( P a / m ) a n d size of rivulet for different channel sizes iv 27 36 List of Figures 1.1 Cross section of a P E M F C 1.2 Cross-section of a rivulet on a flat surface. T h e flow velocity w(x,y) is perpendicular the page 4 2.1 Cross-section of a rivulet on a flat surface 3.1 3.2 3.3 3.4 3.5 3.6 3 to 7 Cross-section of a rivulet on a flat surface R i v u l e t cross sections for L = 5mm, 5cm, 9 = 15° R i v u l e t cross sections for L = 1cm, 9 = 15°, 40° A s y m p t o t i c velocity contours i n m / s for L = 5mm, T = 0,9 = 15° N u m e r i c a l velocity contours i n m / s for L = 5mm, T = 0,9 = 15° A s y m p t o t i c velocity contours i n m / s for L = 5 m m , T = 0.5,9 = 15° 16 17 18 20 20 20 3.7 N u m e r i c a l velocity contours i n m / s for L = 5mm, T = 0.5,9 = 15° 3.8 A s y m p t o t i c velocity contours i n m / s for L = 5mm, T = 1,9 = 15° 3.9 N u m e r i c a l velocity contours i n m / s for L = 5mm,T = 1,9 = 15° 3 . 1 0 A s y m p t o t i c velocity contours i n m / s for L = 5mm, T = —1,9 = 40° 3.11Numerical velocity contours i n m / s for L = 5 m m , T = —1, 9 = 40° 3 . 1 2 F l u x of rivulet for L = lem, T = 0,a = 15°, 9 = 15° 3 . 1 3 F l u x of rivulet for L = lem, T = 0, a = 15°, 9 = 40° 3 . 1 4 F l u x o f r i v u l e t f o r L = l c m , T = l , a = 0,(9=15 3 . 1 5 F l u x of rivulet for L = 1cm, T = 1, a = 0,9 = 40° 3.16Energy of r i v u l e t for L = lem, T = 0, a = 15,9 = 15° 3.17Energy of r i v u l e t for L = lem, T = 0, a = 15,9 = 40° 3.18Energy of r i v u l e t for L = lem, T = 0, a = 0, 9 = 15° 3.19Energy of r i v u l e t for L = 1cm, T = 1, a = 0, 9 = 40° 20 21 21 21 21 22 22 '....22 22 23 23 23 23 o 3.20Energy vs. flux of rivulet for L = lem, T = 0, a = 15,9 = 15° 3.21Energy vs. flux of rivulet for L = lem, T = 0, a = 15, 9 = 40° 3.22Energy vs. flux of rivulet for L = 1cm, T = 1, a = 0,9 = 15° 25 25 25 3.23Energy vs. flux of rivulet for L = lem, T = 1, a = 0,9 = 40° 25 4.1 Geometry of a fuel cell channel 29 4.2 Cross section geometry of a ring rivulet channel • 4.3 Cross section of a circular arc rivulet channel 4.4 Triangular mesh of the cross section domain 30 ' 4.5 Contour plot of fluxes 4.6 shape of rivulet for r = l m m , 6 = 100° 4.7 shape of rivulet for r=.25mm, 6 = 60° 33 34 34 - 37 37 4.8 Rivulet size ip vs. Qratio for R = l m m , 6 = 84° 4.9 A i r flux vs. rp for R = l m m , 6 = 84°, and A = 1 4.10Pressure gradient vs. Qw for annulus and circular arc rivulet 38 38 38 5.1 rectangular geometry 41 1 50 Verification geometry v Analytical and numerical solutions Acknowledgement I a m grateful for the supervision a n d support of D r . B r i a n W e t t o n . I w o u l d also like to thank T i m M y e r s , who develeped the asymptotoics of s m a l l aspect ratio for rivulets described i n the b e g i n n i n g of chapter two. T h i s work is p a r t l y funded by M I T A C S . I also like to t h a n k a l l the friends, colleagues a n d faculty members who have helped me so m u c h along the way. vii Chapter 1 Introduction 1.1 Background W h e n l i q u i d flows d o w n an i n c l i n e d flat surface, several different flow regimes can be observed. T h e flow regimes depend m o s t l y on the flow rate, the p h y s i c a l properties of the liquid-solid surface system a n d the i n c l i n a t i o n of the surface. T h e flow p a t t e r n c o u l d be droplets, linear rivulets, meanders, p e n d u l u m rivulets, or even i n sheets [3]. W e w i l l investigate water flow i n the linear r i v u l e t region i n this thesis. A r i v u l e t is a stream of l i q u i d flowing on a solid surface under a gas w i t h (basically) u n i d i rectional flow p e r p e n d i c u l a r to the r i v u l e t cross section. T h e m a t h e m a t i c a l interest begins w i t h the study of the m o t i o n of a t h i n continuous film d r i v e n by b o d y or surface forces over a d r y substrate. A t some stage, the continuous film w i l l break up at the leading edge into a series of distinct rivulets or droplets. M o d e l i n g this behavior requires the e x a m i n a t i o n of distinct regions, namely, the t h i n continuous film, the m o v i n g contact line a n d the rivulet[2]. T h i n film flow, w i t h one free surface, is already relatively well u n d e r s t o o d by a p p l y i n g l u b r i cation theory[1], However, as the w i d t h of a film decreases into the r i v u l e t region, t h i n film a p p r o x i m a t i o n becomes i n v a l i d . I n this thesis, we w i l l study the steady flow of a single rivulet, driven by inter-facial shear a n d gravity. T h e p h y s i c a l quantities that are u s u a l l y considered are: g r a v i t a t i o n a l forces, incompressble viscous fluid mechanics (often i n a t h i n film a p p r o x i m a t i o n ) , surface tension of the l i q u i d gas interface, and the contact angle at the solid-liquid-gas interface. T h e r e is no "front" to the r i v u l e t a n d therefore no need to address the issue of m o v i n g contact line. 1 R i v u l e t flow plays a n i m p o r t a n t role i n a number of i n d u s t r i a l processes, i n c l u d i n g gas-liquid contacting equipment i n d i s t i l l a t i o n and absorption, d r y p a t c h f o r m a t i o n o n heated surfaces, l i q u i d film drainage from steam t u r b i n e stator blades, etc. Specifically, i n heat exchanges, interfacial shear, caused by the over l y i n g gas flow, plays a significant role. Here, the break-up of a film into rivulets m a y have disastrous consequences [2]. T h e aerospace i n d u s t r y also has a n interest' i n r i v u l e t flow d r i v e n by air shear, w h e n considering the de-icing of air-plane wings. The p r o b l e m occurs because ice formation or snow a c c u m u l a t i o n on airplane wings poses a threat to safety of flight. D u r i n g a p e r i o d after take-off a n d before the l i q u i d is finally b l o w n off, it forms rivulets and flows as "run-back" water w h i c h affect the a e r o d y n a m i c behavior of the wings. K n o w l e d g e of this flow a n d consequently how far back the aerofoil requires heating is essential i n developing a n a n t i - i c i n g strategy. The most recent i n d u s t r i a l interest to the understanding of a r i v u l e t flow arises from m o d e l l i n g the p r o t o n exchange membrane fuel cell ( P E M F C ) , w h i c h is also the i n i t i a l m o t i v a t i o n of this thesis. T h e connections between fuel cells a n d r i v u l e t flow are described i n details below. I n all of the above areas, a knowledge of the fluid m o t i o n a n d s t a b i l i t y characteristics of rivulet flow is required to assess the efficiency of the process. 1.2 Application to fuel cells A fuel cell basically converts electro-chemical into electrical energy w i t h h i g h efficiency, and i n some cases, w i t h zero p o l l u t i o n . T h e r e are m a n y different types of fuel cells, such as the p r o t o n exchange membrane hydrogen fuel cell, alkaline fuel cell, m o l t e n carbonate fuel cell, and direct m e t h a n o l fuel cell etc. A l t h o u g h different types of fuel cells a l l have their unique characteristics, they are a l l environment friendly power generators. F u e l cell technology existed since 1839, but only recently have fuel cells gained p o p u l a r recognition a n d come under serious consideration as a n economically a n d technically viable power source. D u e to increased concerns of p o l l u t i o n caused by automobiles and power plants, several companies have p u t a lot of effort into the development of fuel cells. For example, some car manufacturers are h o p i n g to c o m m e r c i a l l y produce fuel cell powered cars i n 2004. 2 P E M F C s require hydrogen a n d oxygen as inputs, oxygen is often s u p p l i e d by ambient air. T h e reaction is c a t a l y z e d by a t h i n layer of P l a t i n u m , to w h i c h the reactant gases are delivered v i a a series of pressurized flow fields on either side of N a t i o n membrane. T h e membrane is permeable to only s m a l l positive ions, like the H + ions, a n d is s u p p o r t e d by a gas diffusion e l e c t r o d e ( G D E ) w h i c h is currently c o m m o n l y made from a porous teflonated c a r b o n fibre paper. T h e G D E allows reactant gases to reach the active catalyst sites on the fuel cell membrane a n d carries current away from the sites. F i g u r e (1.1) below illustrates the structure of a P E M F C o n the cathode side. O n top of the P l a t i n u m coating, there is the anode G D E a n d then the Pkimiiuun Nation coating Oxygen flow filed plate F i g u r e 1.1: Cross section of a P E M F C . hydrogen channel. D u r i n g the reaction process, water is formed on the cathode side a n d accumulates i n the oxygen channel. T h o u g h water is i m p o r t a n t to m a i n t a i n the h u m i d i t y level of the channel so the membrane r e m a i n well hydrated, excess amount of l i q u i d water i n the gas flow channels can reduce the efficiency of the fuel cell significantly by b l o c k i n g the channel, a n d so i n h i b i t i n g the diffusion of gas t h r o u g h the m e d i u m , thereby prevents reaction from o c c u r r i n g . T h u s "water management" inside the fuel cell is very i m p o r t a n t to achieve h i g h efficiency. Since the excess amount of water i n channels is carried out using the gas shear flow, it places an e x t r a load o n the compressor used to p u m p reactant gasses throught the fuel cell. T h i s is an i m p o r t a n t design consideration. 3 1.3 Mathematical formulation Consider a N e w t o n i a n rivulet flow d o w n a channel w i t h cross section b o u n d a r y T. F i u g r e 1.2. T h e surface of the rivulet h(x) shown i n divides the d o m a i n into region 1 a n d region 2. r F i g u r e 1.2: Cross-section of a rivulet on a flat surface. T h e flow velocity w(x, y) is perpendicular to the page. T h e d y n a m i c s of a viscous incompressible fluid are governed by the Navier-Stokes equation [6]. P^r + (u • V ) u = / - Vp + fiV u 2 together w i t h the continuity equation V • u = 0. T h e terms are: • p: density of fluid • u: velocity of fluid • p,: viscosity of fluid • / : external force ( i n c l u d i n g gravity) • p: pressure gradient T h e rivulet flow is modeled under the following assumptions: 1. the fluid is incompressible, 4 (1.1) 2. the flow has reached a steady-state, 3. the flow is fully developed, that is, it is u n i d i r e c t i o n a l a n d the derivatives i n the direction of flow are negligible. A s s u m p t i o n s 2 a n d 3 indicate the velocity field is specified by u = (0,0, w(x,y)), a n d we also assume that gravity is the o n l y external force, so the Navier-Stokes equation reduces to - V p + MV u + pG = 0 (1.2) 2 where G = (G\, G , G3) = pg(0,g -y,g • z) after decomposing the g r a v i t a t i o n a l force into y and 2 z directions (z denotes the unit vector i n z direction). C o m p o n e n t wise, we have: dx - ^ +G dy 4- A- 11 +liAw dz dp 2 +G 3 = 0 (1.3) = 0 (1.4) = 0 (1.5) Note that the continuity equation is a u t o m a t i c a l l y satisfied. E q u a t i o n s (1.3) a n d (1.4) can be solved to obtain: , p = G y + C (z) + d. 2 0 (1.6) T h i s can be substituted into equation (1.5), resulting i n : -C' (z) + pAw + G = 0. 0 3 (1.7) Since w is independent of z, and Co is a function of z only, b o t h terms have to equal to a constant i n order for the equation to h o l d . L e t —C' (z) = A, then 0 \ + pAw + G = 0 (1.8) p = -Xz + G y + C (1.9) 3 Also, 2 A p p r o p r i a t e b o u n d a r y and interface conditions for the system are: 5 (a) N o - s l i p at the fluid-solid b o u n d a r y T, (1.10) w{x,y) = 0, (b) C o n t i n u i t y of velocity at free surface y = h(x), (1.11) [w] — W\ — W2 = 0 (c) C o n t i n u i t y of shear stress at free surface, H[w ) = /J,(w n nl Wn2) = 0 (1.12) (d) B a l a n c e of n o r m a l stress at interface, \p]=P2-Pi~ a K (1.13) where K is the curvature of the surface. If it contacts b o u n d a r y w i t h contact angle 0\ = #2, then K = (1 + ^)3/2 We w i l l consider special cases of rivulet flow i n the later chapters. I n chapters 2 a n d 3, we w i l l study r i v u l e t flow d o w n a n i n c l i n e d plane, using b o t h n u m e r i c a l a n d a s y m p t o t i c techniques. T h e a n a l y t i c a l s o l u t i o n requires the r i v u l e t aspect ratio to be s m a l l . T h e n u m e r i c a l investigation is carried out w i t h o u t these restrictions. T h i s allows us to verify the a n a l y t i c a l solutions a n d also ascertain their region of validity. I n chapter 4, we assume r i v u l e t flow d o w n a flat circular channel, a n d investigate the relationship between fluxes a n d pressure gradient needed to drive the flow, w h i c h is of p a r t i c u l a r interest to fuel cell design. M A T L A B a n d F E M L A B [11] are used to solve the problems numerically. F i n a l l y , the i m p o r t a n t question of whether a rivulet w i l l break-up into two or more sub-rivulets is addressed. A l i m i t e d s t a b i l i t y study based o n energy c r i t e r i o n is done i n chapter 3, a n d a more rigorous linear s t a b i l i t y analysis is shown i n chapter 6. 6 Chapter 2 Rivulet Flow Down an Inclined Plane Asymptotic Analysis Consider a r i v u l e t of w i d t h L and height y = h(x) flowing d o w n an i n c l i n e d surface i n c l i n e d at angle a to the horizontal as shown i n figure (2.1) below. L F i g u r e 2.1: Cross-section of a rivulet o n a flat surface. I n this case, the p r o b l e m is simplified i n the following sense. 1. It is assumed that air provides the tangential stress T only, so b o u n d a r y c o n d i t i o n (1.11) is not a p p l i e d . C o n d i t i o n (1.12) becomes \xw n 2. A i r pressure is constant, so p = p a p = Gy 2 = T. + OK for a l l z. T h u s equation (1.9) becomes + C = G h(x) + C = p 2 a + <JK (2.1) Differentiate the equation above, to o b t a i n the L a p l a c e - Y o u n g equation, r — - - —( 8x~ dx\(l 2 a \ + hl)Wj hxx (2.2) w i t h G2 = pgcos(a). Besides the previous assumptions, we also assume the aspect ratio e = H/L height scale) is small, so that terms of 0(e ) 2 (where H is the are negligible, a n d do a n a s y m p t o t i c analysis of the p r o b l e m . F i r s t , we non-dimensionalize the components of the Navier-Stokes equations, a n d they are: Px w - B + ew 2 yy xx p y W h e r e B — pgH sma/(p,W) 2 = 0, (2.3) = 0 , (2.4) = - 1 . (2.5) is the B o n d number, represents the ratio of g r a v i t a t i o n a l to viscous forces. A l l variables are related to the dimensional quantities (denoted by hats) by p —> pgH cos a p y -» Hy x-^Lx w —>• Ww T = ^ - T T h e choice of velocity scale, W, depends on whether gravity or surface shear is the p r i m a r y d r i v i n g force. If gravity dominates the flow then W = pgH s'm a/and the B o n d number 2 5 = 1. If air shear drives the flow then W = (HT/p,). T o keep the most general m o d e l , the velocity scale w i l l be left unspecified i n the following analysis. A p p r o p r i a t e b o u n d a r y conditions become: / ~\ ~ 2 1 dw (dh ,1 , 2 dy \ dx I dhdw " \dxdx dw dy 0 p = -C- h [l-e h 1 2 xx 2 x ^ m y=h + ---] dx . =0 (2.6) x=0 (2.7) These represent no-slip at the solid surface, continuity of shear stress at the free surface, symmetry about the center-line a n d continuity of n o r m a l stress at the free surface respectively. T h e non-dimensional number C~ l = a/(pgL 2 cos a) is the inverse c a p i l l a r y number, where C resp- resents the ratio of viscous to g r a v i t a t i o n a l forces. Note, we assumed eh x expansion w h i c h requires a relatively s m a l l l i q u i d / s o l i d contact angle. <C 1 i n the pressure T h e form of the system (2.3)-(2.7) indicates there exists a series s o l u t i o n i n powers of e for 2 h,w a n d p. Integration of (2.5), subject to (2.7), immediately, gives the series for pressure i n terms of the u n k n o w n functions hf. p=(h -y)- C - % 0 x x +e (2.8) + C- {h hl - h )] +••• . 2 1 0xx lxx T h e terms i n the film height series are determined by s u b s t i t u t i n g for pressure i n (2.3), ho ~Ch = 0, hi -Ch = ^[h h ] xxx 0x xxx lx (2.9) 0x (2.10) • 2 0xx T h e appropriate b o u n d a r y conditions are, Mo) = o, M i ) = o (2.ii) and either jf = 2^ff , J\idx = 0 (2.12) or hoM) = - ^ - , h (l)=0. (2.13) lx C o n d i t i o n s (2.11) state there is s y m m e t r y about the m i d - p o i n t , at the edge the height is zero. C o n d i t i o n (2.12) relies on the cross-sectional area, A, b e i n g a measurable quantity, while (2.13) uses the contact angle. Solve for hi a n d get: ho = /ii = K(COS\I(VCX) - c o s h ( v C ) ) , (2.14) / (cosh(VCx) - c o s h ( v C ) ) + f{x) - / ( l ) , (2.15) / ai where K C wh(VCx){V^x 3 S /(*) /g(l) 1 ~ _ /c sinh{y x) ^c cosh{ x)) ) ( 2.i ) 6 i^. i j\ VCsmh(VC) ' (2.18) and K = —Ltan9/HVCs'mh(VC). T h e height-scale H is chosen such that /IQ(0) = 1, i n w h i c h case, Ltan6>(cosh(y C) - 1) / 7= Tl — K — ' '" ;= , fCsmh(VC) (2.19) cosh( C) - 1 ' / v T h e fluid velocity is found by integrating (2.4), subject to (2.6) a n d a p p l y o n y = ho + e h\-\ , 2 w(x,y) = b y + B,y— + e hy + Bh 0 — 0xx + (2.20) where b = T-Bho h = 0 - B h , 1 (2.21) - ? ^ - B h o h l 2 T -•-•-u, ^ h + (2.22) 2 T o highlight the effect of the d r i v i n g forces on the fluid flux a n d p e r m i t a relatively simple a n a l y t i c a l solution, o n l y the leading order terms w i l l be considered from now on. T h e fluid flux d o w n the slope, Qo = Qo/{LHW), is defined as -1 rh (x) 0 w (x,y)dAo = - 2 / / Jo Jo w (x,y)dydx 0 Jo 0 T%x B- (2.23) (2.24) T h i s integrates to Use 2s SVC + ' QVC 3c 2 2 Qo = T} BK " J - TK — cr 2 3cs 1 2V3 2VC ' 22J ' + ( 2 - 2 5 ) where s = s i n h s/C a n d c = cosh \[C. L i m i t i n g values for s m a l l C (narrow rivulets, low density or high surface tension) a n d large C can now be calculated. For s m a l l C , (2.26) for large C Ke^ 2 Qo = 2 2nBe^ 24 + 3r (2.27) T h e role of inter facial shear on flux can be clearly seen from the preceding equations. For example, i n the s m a l l C l i m i t (and note that K < 0) a positive air shear, T > 0, acts against 10 gravity a n d so decreases the flux d o w n the slope. W h e n T = —2KBC/7 the a i r shear exactly balances gravity a n d the net flux is zero. A value of T below this w i l l result i n a fluid flux d o w n the slope. In d i m e n s i o n a l form, the l i m i t i n g expressions (2.26), (2.27) for flux become = -^^^LH^e-^LH^e Q \x 105 = 5 Q ^ (2.28) 15// ^ < - ^ ' ^ ' . <«»> respectively. It is interesting to note that the s m a l l C solution, equation (2.28), associated w i t h h i g h surface tension, i n fact does not involve surface tension. A n o t h e r point of interest is that the above results depend either u p o n t a n 9 or t a n 9, i n d i c a t i n g that accurate measurement of the contact 3 2 angle is crucial, p a r t i c u l a r l y for values of 0 m u c h less or m u c h greater t h a n 4 5 ° . As L — > oo, the height-scale H is identified as the a p p r o x i m a t e l y constant central height (the l i m i t i n g value of (2.19a)), g i v i n g H = tan 9 a / (pg cos a) a n d C = L t a n 9/H . 2 2 2 Substituting this into (2.29) leads to a n expression for the d i m e n s i o n a l f l u x / u n i t w i d t h : Qo pgH sma TH 3 2L 2 2u ' 3p [ ' 1 w h i c h is the result obtained from the two d i m e n s i o n a l analysis for a fluid sheet given i n M y e r s & T h o m p s o n [5]. The e n e r g y / u n i t length of a r i v u l e t involves kinetic, p o t e n t i a l a n d surface energies a n d also viscous d i s s i p a t i o n . T h e s u m of these must r e m a i n constant, a l t h o u g h the i n d i v i d u a l terms may alter. U n d e r the current m o d e l assumptions a n d p r o v i d e d the flow is isothermal, it is clear that as the r i v u l e t moves d o w n the slope its k i n e t i c a n d surface energies r e m a i n constant since w a n d h are independent of z. P o t e n t i a l energy w i l l be lost a n d this must be balanced b y viscous d i s s i p a t i o n . T h e s u m of kinetic a n d surface energies per u n i t length is, E = E pLHV 1 f 1 2 = - f h J q 2 j J , Pa + v dydx + 11 2(a -a ) sl sg , (2.31) where A is the cross-sectional area, a \ a n d a s are the s u r f a c e / l i q u i d a n d surface/gas interfacial sg tensions, P is the perimeter of the liquid-gas interface. T a k i n g the leading order velocity terms from equation (2.20) a n d setting h = ho the kinetic energy per u n i t length, K, is a s y m p t o t i c a l l y rho "1 / „ 2 \ 2 = Ji11° fa + j) K B = f jhl +^h* b 10 + v d (- ) dx 2 32 ^hldx. (2.33) , (2.34) S u b s t i t u t i n g for ho (x) a n d integrating leads to K=(a + T + aT) 2 0 a i 2 where - 6 4 s - 6 0 7 c s + 120c VC 2 BK 2 a 0 — 5 5 900 55«; 4 72c - 2 7 4 c s + 600c VC + VC 4 3 225cVC + 9 VC - 5 0 s c + 24c y/C - 55cs 2 3 A 288 /•c 6 c m L - 4 s + 9c\/C- llc s 18 y/C 3 a-2 3 2 T h e perimeter length, to leading order, is the same as the r i v u l e t w i d t h P = 2 [ y/1 + e h%dx = 2 + 0(e ) . Jo l 2 2 (2.35) T h e Y o u n g - D u p r e s equation[12]provides a relation between the surface tensions: 0~sl ~ sg = —cr cos (9 . (2.36) a A l t h o u g h , as is frequently pointed out, this implies a unique contact angle w h i c h m a y only be achieved o n extremely s m o o t h surfaces, hence it does not account for hysteresis. T h e leading order a n a l y t i c expression for the energy E is therefore E =K 2a + ^(l-cose) T pHV* where K is specified by equation (2.34). 12 , (2.37) It is at this stage that the advantage of the current m e t h o d over l u b r i c a t i o n theory becomes apparent; the k i n e t i c energy t e r m enters the energy balance at leading order. It w o u l d only do so i n the l u b r i c a t i o n a p p r o x i m a t i o n p r o v i d e d certain other restrictions are satisfied. T o a p p l y — 8 pVL/p l u b r i c a t i o n theory the reduced R e y n o l d s number S R 2 2 e ~ pHLV 2 •< p /(8 p) 2 3 (in order to satisfy 5 R 2 e <C 1. T h e k i n e t i c energy <C 1). T h e surface energy ~ aL. T a k i n g t y p i c a l values for a water r i v u l e t of dimensions L x H = l c m x 1 m m shows the k i n e t i c energy must then be significantly less t h a n 1 0 . T h e surface energy is 0 ( 1 C T ) . T h e k i n e t i c energy w o u l d - 6 3 therefore not enter the leading order energy balance i n this case. I n fact, for water flow the energy c a l c u l a t i o n w i l l only be v a l i d under l u b r i c a t i o n theory w h e n 8 H <C 1 0 2 _ 8 m. For b o t h shear d r i v e n ( a = 0) a n d gravity d r i v e n ( T = 0) large L flow, the a s y m p t o t i c fluxes can be derived from equations (2.28) a n d (2.29), a n d the a s y m p t o t i c energy behaviour c a n be derived from equation (2.37) a n d we find Q oc L , a n d E oc L , so we conclude: £ oc Q for b o t h shear a n d gravity driven, large L flow. T h e results also make sense physically. For L large, the r i v u l e t w i l l be like a flat sheet for most of the region (see F i g u r e 3.2 for r i v u l e t shape when L = 5 c m ) . So E a n d Q b o t h vary linearly w i t h length L. T h e leading order a s y m p t o t i c behaviour is different for s m a l l L. T h e equations derived above won't a p p l y i n this case since they are only v a l i d for e s m a l l w h i c h corresponds to large L. We c a n get the a s y m p t o t i c results for s m a l l L as follows. F o r L s m a l l , the cross section of a rivulet looks like a semi-circle. T h e dominate energy is the surface energy w h i c h p r o p o r t i o n a l to arclength w h i c h is TTL i n this case. So we have E oc L w h i c h is v a l i d for b o t h shear a n d gravity driven flow. For fluxes we need to consider the two cases separately. F i r s t nondimensionalize the general equation pAw = A. L e t x —> Lx, y —> Ly, where x a n d y are n o n d i m e n s i o n a l quantities. T h u s dw dx 2 2 dw dy 1 2 L 2 2 13 dw dx 2 2 dw dy 2 2 For gravity d r i v e n flow, we are solving Aw conditions, where G oc area = irL . 2 = G i n the d o m a i n w i t h Dirichelet boundary T o get flux, we need to integrate w over the d o m a i n , the area of the cross section. It results another factor of L . 2 Q oc L T h u s we get 4 for this case. For shear d r i v e n flow, we solve Aw = 0 i n the d o m a i n a n d dw/dn = TL o n the circular arc and Dirichelet c o n d i t i o n o n the b o t t o m . So we have Q oc L 3 I n summary, we should expect E oc Q 1 / 4 for gravity d r i v e n s m a l l L flow, a n d E oc < 2 1/3 for shear d r i v e n s m a l l L flow. These results are reflected i n the n u m e r i c a l solutions i n chapter 3, the asymptotics are also needed to fill i n the gaps i n the n u m e r i c a l results. 14 Chapter 3 Rivulet Flow On A n Inclined Plane Numerical Method 3.1 Numerical method Consider a rivulet flows d o w n an inclined plane shown i n F i g u r e 2.1. T o solve the p r o b l e m numerically, we compute the cross section of the rivulet first by a p p l y i n g the L a p l a c e - Y o u n g equation (2.2). pg cos(a)h = -a4z x ( (3.1) w h i c h governs the height h(x)of the rivulet. Integrate the equation w i t h respect to x, we obtain: h pgcos{a) xx (1 + hi) / a 3 2 h + C (3.2) where C is a constant that is to be determined along w i t h the profile h(x) for given symmetric half length L [h(±L) = 0] a n d surface contact angle 0 [tan6 = h (±L)]. x Note that the left h a n d side of (3.2) is the curvature of the rivulet. E q u a t i o n (3.2) is a second order differential equation. T o reduce it to first order, we parametrize the equation i n arclength s. L e t x = x(s) and y = h = y(s). \{x ,y )\ y = C Then, (3.3) a Also, |(x',y')| = l 15 (3.4) Differentiate (3.4) get: x x + yy (a; ) + (y ) 2 2 = 0 = (C + (3.5) . y) 2 (3.6) a Introduce a new variable <f>, the angle of the tangent line o n the surface of the r i v u l e t , shown i n F i g u r e 3.1. -L 0 F i g u r e 3.1: Cross-section of a rivulet on a flat surface. T h e n x = cos(c/>) a n d y — sin(c/>). T h e curvature K satisfies: K 2 = (x") + (y") = ( - s i n ( c 6 ) ^ ) + (cos(ci)cA') 2 (3.7) 2 2 2 (3.8) (3.9) T h u s equation (3.2) becomes: cos(c/>) y I = I 1 \<P J sin(c/.) f~i 1 (3.10) po cos a „ \ C. + vs-z—y J w i t h the i n i t i a l conditions(s=0): ' x(0) \ V y(0) #0) ( -L ^ = ) 0 \ 0 (3.11) ) T h i s arc-length p a r a m e t r i z e d formulation not only simplifies the e q u a t i o n to first order, but also enables us to compute r i v u l e t shapes w i t h contact angles greater t h a n 90° w i t h o u t a d d i t i o n a l 16 work, a l t h o u g h this c a p a b i l i t y is not used i n this s t u d y where we are a i m i n g at c o m p a r i n g to asymptotics suited for s m a l l angles. Cross section T h e O D E solver c a p a b i l i t y i n M A T L A B is used to compute the cross section geometry. We used a shooting m e t h o d to determine C such that <f>{s*) = 0 where s* satisfies x(s*) = 0, a s y m m e t r y c o n d i t i o n that w o u l d i m p l y the original c o n d i t i o n at x = L . F i g u r e 3.2 below shows the r i v u l e t cross-section obtained from b o t h the a n a l y t i c a l formula, equation(2.19) a n d the n u m e r i c a l m e t h o d described above, for water flowing d o w n a plane w i t h a = 6 = 15°. L takes the values of 5 m m and 5 c m . N o results are shown for L < 5 m m since these are v i r t u a l l y i d e n t i c a l to the 5 m m graph, s i m i l a r l y wide r i v u l e t s a l l resemble the 5cm graph. T h e figure shows that narrow rivulets, shaped by surface tension, assume a n a p p r o x i m a t e l y c i r c u l a r profile. A s the w i d t h increases a flat zone develops i n the centre, where gravity determines the shape a n d surface tension plays a significant role only near the edges. ' 1 1 1 1 ' 1 0,9 0.8 —. | J — 1 1 L=5cm ^~v>. numerical j asymptotic | 0.7 L=5mm y 0.6 I 05 C - \1 0.4 0.3 • 0.2 • 0 0.1 0.2 0.3 0.4 0.5 lenglh(m) 0.6 0.7 0.8 0.9 1 F i g u r e 3.2: R i v u l e t cross sections for L = 5 m m , 5 c m , 0 = 15°. F i g u r e (3.3) shows a comparison of the n u m e r i c a l a n d a s y m p t o t i c results of the cross-sections of rivulets for different contact angles w i t h a = 15° and L = 5 m m . T h e a s y m p t o t i c results follow the leading order height, equation (2.19). W i t h a contact angle 9 = 15° the two solutions are i n good agreement, w i t h a m a x i m u m error at L = 0 of a p p r o x i m a t e l y 3%. A s the contact 17 angle increases, the solutions diverge. W h e n 9 = 4 0 ° the m a x i m u m error increases to around 16%. T h i s provides a guideline for w h e n the a s y m p t o t i c s o l u t i o n is v a l i d . F i g u r e 3.3: R i v u l e t cross sections for L = lcm,6 = 1 5 ° , 4 0 ° Velocity contour N e x t , we used the P D E t o o l b o x i n M A T L A B to compute the velocity w given the P D E : Aw = pg sin a (3.12) with boundary condition w(x,0) = 0 and c o n t i n u i t y of tangential stress T dw dn at the free surface y = h(x). T h e P D E t o o l b o x a u t o m a t i c a l l y grids the specified d o m a i n a n d w i l l s u b d i v i d e the o r i g i n a l g r i d for more accurate c o m p u t a t i o n s as requested. Numerical results are c o m p u t e d to at least 1% accuracy based o n a refinement study. For example, i n the L = 1cm, 9 = 15° r i v u l e t c o m p u t a t i o n , a t o t a l of 128 t r i a n g u l a r grids are used i n the d o m a i n to give this accuracy. T h e t o o l b o x uses standard linear finite element m e t h o d . We consider two representative (T ^ 0,a cases i n the computations. O n e is the shear d r i v e n case = 0), a n d the other is the gravity d r i v e n case w i t h zero shear (T = 0, a ^ 18 0). T h e a s y m p t o t i c a n d n u m e r i c a l velocity contours are shown i n Figures 3.4-3.11 below. T h e a s y m p o t o t i c results follow equation (2.20). T h e axes of the a s y m p t o t i c velocity contours are non-dimensional, w h i l e they are dimensional for the n u m e r i c a l values. However, the velocity contour values are d i m e n s i o n a l for b o t h a s y m p t o t i c a n d n u m e r i c a l computations. I n a d d i t i o n , velocity contours of o n l y half of the r i v u l e t are shown for the a s y m p t o t i c computations, while the contours are shown for the whole rivulet for the n u m e r i c a l ones. O n e i m p o r t a n t t h i n g to notice is that the velocity values obtained from the a s y m p t o t i c formula a n d those obtained from the n u m e r i c a l calculations are very close to each other. However, this is not the case for fluxes and energies w h i c h w i l l be shown later. T h e parameter values used i n the computations are: • p: density of water (lOOOfcg/m ) 3 • /J,: viscosity of water (1.0 * I0~ kg/m/s) 3 • a: surface tension between water a n d air (72.4 * 10~ kg/s ) 3 2 • g: gravity acceleration ( 9 . 8 1 m / s ) 2 • a: angle of i n c l i n a t i o n (15°) For the 5 m m r i v u l e t , velocity contours for zero inter-facial shear are shown i n Figures 3.4 a n d 3.5. T h e velocities decrease from a m a x i m u m of —0.5 (i.e. i n the d i r e c t i o n of gravity) at the apex to 0 at the substrate. T h e velocity contour values are negative since the flow is d r i v e n b y gravity w h i c h is i n the negative direction. Note, that at the liquid-gas interface, the contour lines are almost vertical, i n keeping w i t h the c o n d i t i o n dw/dy = T + 0(e ) ss 0. 2 19 io" 4 Velocity Contour for T=0. alpha=15 Velocity Contour T=0, alpha= 15 55T 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 F i g u r e 3.4: A s y m p t o t i c velocity contours i n m / s F i g u r e 3.5: N u m e r i c a l velocity contours i n m / s for L = 5 m m , T = 0,9 = 15° for L = 5mm, T = 0,9 = 15° If the inter-facial shear is increased, then the surface components of the fluid w i l l move against gravity, a n d a sufficiently h i g h values of T pushes a l l the fluid against gravity. T h i s resembles the familar s i t u a t i o n of rivulets being b l o w n up the windscreen of a fast m o v i n g car. O n F i g u r e 3.6 a n d 3.7 a n example is shown w i t h relatively low shear ( T = 0.5), so the central part of the rivulet moves w i t h gravity. T h e zero velocity contour, w h i c h starts at the apex, can be seen to split the rivulet; at the edges the fluid flows i n the direction of the air, near the centre gravity dominates a n d the flow is d o w n the plane. Velocity Contour T=0.5, alpha=15 -, \ \ ' ' ' _ 0. O. 0_ V \ \ m .-0.0620.043- 0 0 1 02 03 04 05 0.6 07 08 09 F i g u r e 3.6: A s y m p t o t i c velocity contours i n m / s F i g u r e 3.7: N u m e r i c a l velocity contours i n m / s for L = 5 m m , T = 0.5,9 = 15° for L = 5mm, T = 0.5,9 = 15° Increasing T further the gravity d r i v e n region becomes smaller a n d finally disappears. 20 Subse- quently the velocity contours become m u c h flatter, approaching that of a t h i n sheet rather t h a n a rivulet, i.e. the velocity w ~ w(y). Increasing T further strengthens this effect, u n t i l the x dependence becomes negligible, except for i n a b o u n d a r y layer at the rivulet edge. A sample is shown i n Figures 3.8 a n d 3.9 for T = 1. Velocity Contour T=0, alpha=15 Velocity Contour T= 1, alpha= 15 -4 -5 -3 F i g u r e 3.8: A s y m p t o t i c velocity contours i n m / s F i g u r e 3.9: N u m e r i c a l velocity contours i n m / for L = 5 m m , T = 1,9 = 15° for L = 5mm,T =1,9 = 15° For T < 0, i.e. negative shear stress, w h i c h is i n the same d i r e c t i o n as gravity, we expect larger negative velocity values compare w i t h zero shear stress case. F o r example, w i t h T = —1, the m a x i m u m negative velocity is - 1 . 0 , as shown i n Figures 3.10 a n d 3.11. „ Velocity Contour T=-l, alpha=15 $ f 4 1 0 Velocity Contour T=-l,alpha=15 » / ? 0 88 .--an''' -0.66 - -0.56 -0 55 -0.44 A* -0.33 -0.33 ^ -0.22 0.11 a M 0 0.1 0.2 F i g u r e 3.10: 03 0.4 0.5 0.6 0.7 OB 0.9 A s y m p t o t i c velocity contours i n m / s for L = 5mm,T = -1,9 = 40° F i g u r e 3.11: N u m e r i c a l velocity contours i n m / for L = 5mm,T = -1,9 = 40° 21 F l u x and energy T h e asymptotic a n d n u m e r i c a l flux results are shown i n Figures 3.12-3.15 as a function of rivulet length L. T h e asymptotic flux follows equation o 0.001 0.002 (2.25). o. F i g u r e 3.12: F l u x of rivulet for L = lem, T F i g u r e 3.13: F l u x of r i v u l e t for L = lem, 0,a = 15°, 9 = 15° 1 | 0 F i g u r e 3.14: F l u x of r i v u l e t for L = lem, l , a = 0,9 = 15° T 0,a = 15°, 0 = 40° T 0.001 0.002 0,003 0,004 0,005 length 0.006 0.007 0.006 asymptotic 1 numerical f 0.009 0.01 F i g u r e 3.15: F l u x of r i v u l e t for L = lem, T l , a = 0,6 = 40° I n F i g u r e 3.12 the flux v a r i a t i o n w i t h w i d t h is shown for a r i v u l e t d r i v e n solely by gravity, w i t h a contact angle 8 = 15°. T h e two curves are very close but slowly diverging as L increases. I n this case the m a x i m u m h a l f - w i d t h is 1cm, at w h i c h stage the percentage difference between the curves is 9.98%. I n F i g u r e 3.13 the contact angle increases to 4 0 ° . T h e two curves diverges m u c h faster compare to F i g u r e 3.12. T h i s shows strong dependence of 9 value i n the 22 asymptotics. Figures 3.14 a n d 3.15 show the flux curves for an inter facial shear d r i v e n r i v u l e t ( T = 1) w i t h the contact angle t a k i n g values of 9 = 15°, 4 0 ° . In each case the curves diverges w i t h increasing length. W h e n 9 = 15° the m a x i m u m error is around 0.4%, a n d the m a x i m u m error increases to 11.3% for 6 = 4 0 ° . T h e energy results as a function of length L are shown below i n Figures 3.16-3.19. T h e asymptotic energy follows equation (2.37). F i g u r e 3.16: E n e r g y of rivulet for L = 1cm, T = 0,a = 15,61 = 15° F i g u r e 3.17: E n e r g y of r i v u l e t for L = l e m , T 0, a = 15,9 = 40° :| — :| * asymptotic 1 numerical | i i i i : ]| i,—i—. asymptotic numsrical | y'. \ y ' / y' / ; y 0 0.001 0.002 0.003 0.004 0.005 length 0.006 0.007 0.008 0.009 0.01 F i g u r e 3.18: E n e r g y of rivulet for L = l e m , T = 0, a = 0,9 = 15° F i g u r e 3.19: E n e r g y of r i v u l e t for L = lem, T 1, a = 0,9 = 40° T h e curves resembles the curves of flux vs. length for the gravity d r i v e n energy plot (Figures 23 3.16 a n d 3.17), i.e. the error increases w i t h increasing length a n d 9 value. However, the results are not quite as expected for the shear d r i v e n case (Figures 3.18 a n d 3.19). T h e error diverges as length increases up to a c r i t i c a l value a r o u n d L = 3mm, then there seems to be a constant offset between asymptotic energy a n d n u m e r i c a l energy as L increases further. T h i s is due to the fact that surface energy is dominant w h e n there is no gravity, w h i c h is p r o p o r t i o n a l to length L of the rivulet for L large. 3.2 Stability Analysis Based On Energy Criterion In 1990, P . S c h m u k i a n d M . L a s o [3] studied the s t a b i l i t y of rivulet flow using a simplified rivulet m o d e l , i.e. treating the shape of the rivulet as a circular arc. T h e m o d e l is based on a n energy criterion: the rivulet ceases to be stable when the t o t a l energy of the single rivulet can be lowered by its decaying into several smaller ones (subrivulets). T h e strategy is to formulate the t o t a l energy of the system a n d then to find the m i n i m u m value of this energy as a function of the number of sub-rivulets for a given t o t a l flow rate. We w i l l borrow the idea a n d study the s t a b i l i t y of a single rivulet vs. two sub-rivulets. For a rivulet of a specified flux Q, it w i l l be unstable to breakup into two smaller rivulets i f E(Q) > E(Q - q) + E(q) for some 0 < q < Q. W e denote by Q the m i n i m u m of a l l such fluxes Q a l l o w i n g energetically c favorable breakup a n d let L c length less t h a n L c be the corresponding rivulet length. C l e a r l y then, rivulets of are stable to breakup a n d those w i t h length greater t h a n L c are unstable. For gravity driven, zero shear flow, we compute the rivulet energy as a function of flux as shown i n Figures 3.20 a n d 3.21. Note that E oc Q for large L a n d the dependence of E a n d Q / 1 s m a l l L can be seen clearly from the plots. B u t we should expect E oc Q / 1 4 4 for as predicted by a s y m p t o t i c analysis i n chapter 2. For the shear d r i v e n a n d no gravity flow, the energy vs. flux results are shown i n Figures 3.22 and 3.23. T h e qualitative behaviour also fits our a s y m p t o t i c analysis i n chapter 2, i.e. E oc Q 24 1 / / 3 F i g u r e 3.20: E n e r g y vs. flux of rivulet for L = F i g u r e 3.21: E n e r g y vs. flux of rivulet for L l e m , T = 0, a = 15,0 = 15° 1cm, T = 0, a = 15,6 = 40° for L s m a l l , a n d E oc Q for large L. x 10" F i g u r e 3.22: E n e r g y vs. flux of rivulet for L = F i g u r e 3.23: E n e r g y vs. flux of rivulet for L lcm,T = l,a = 0,6 = 15° l e m , T = I,a = 0,0 = 40° E x a m i n e Figures 3.20-3.23, we notice that the concavity of the E(Q) g r a p h for the gravity d r i v e n case changes from concave up to concave d o w n near the origin, w h i l e the g r a p h is always concave d o w n for the shear d r i v e n case. W e can conclude that the gravity d r i v e n case has p o t e n t i a l for breaking up into smaller rivulets, but it is never energetically favourable to break up into smaller rivulets for the shear d r i v e n rivulets. T o be precise, we emphasize that only one k i n d of i n s t a b i l i t y is considered i n this analysis (break up into m u l t i p l e rivulets) a n d that other 25 mechanisms (meander, breakup into droplets) c o u l d lead to instabilities i n the shear d r i v e n case, a n d may lead to earlier instabilities i n the gravity d r i v e n case. Asymptotic inflection points Since change of concavity leads to instability, and concavity changes at inflection points, we search for the inflection point using the asymptotic energy formula. T h e n we compute energies for different r i v u l e t sizes to see i f it is energetically favorable for a r i v u l e t to break up. Since we only have formulas for E a n d Q as a function of L , a n d it is impossible to solve E as a function of Q e x p l i c i t l y due to the c o m p l e x i t y of the equations, W e a p p l y the chain rule to find the inflection points w h i c h corresponds to Q values such that E (Q) = 0. Since, Q'(L)E"(L) E i { Q ) = - E'(L)Q"(L) f ~ ( 3 - 1 3 ) so E"(Q) = 0 Q' (L)E" (L) - E' (L)Q" (L) = 0 (3.14) M a p l e is used to do the root searching, a n d find that the inflection point occurs at about L c = 1.4mm for 9 = 4 0 ° , a n d L c = 2 . 6 7 m m for 9 = 15° for the gravity d r i v e n case. This indicates that rivulet is more stable for s m a l l contact angles since the inflection point occurs at larger values of L. T h i s is what we expected. N e x t we go o n to find the c r i t i c a l point where the rivulet breaks. Numerical and asymptotic critical points W e d i d a spline fit using the d a t a E(a) and Q(a) (both n u m e r i c a l a n d a n a l y t i c a l data) to o b t a i n a n n u m e r i c a l a p p r o x i m a t i o n of the function E(Q). intervals, a n d compute E(Qi) for each Qi,i T h e n we d i v i d e Q into N equally spaced = 1 • • • N. T h e inflection point Qj at specific L was located. S t a r t i n g from Qi+i upwards we searched for the smallest c r i t i c a l flux Q such that c the following equation is satisfied: Q c = i n f { Q : E(Q) > E(q) + E{Q-q),0<q< 26 Q} (3.15) c 9 method Q 10° asymptotic 4.9513e-07 7.3mm numerical 6.6490e-07 8.5mm asymptotic 3.0942e-07 4.0mm numerical 5.0016e-07 4.8mm asymptotic 5.3161e-07 3.6mm numerical 7.2570e-07 4.1mm asymptotic 6.4747e-07 1.9mm numerical 9.1391e-07 1.9mm asymptotic 8.5269e-07 3.3mm numerical 1.1602e-06 4.5mm asymptotic 7.6853e-07 1.9mm numerical 1.1721e-06 2.5mm 15° 20° 30° 40° 45° L c Table 3.1: N u m e r i c a l estimates of Q c and L c c for g r a v i t y d r i v e n r i v u l e t . For the gravity d r i v e n case, the c r i t i c a l flux Q c a n d the corresponding c r i t i c a l length L c for 6 = 15° a n d 9 = 40° are listed i n table 3.1.The d a t a shows that b o t h Q c and L c get smaller as 9 increases. T h u s rivulets w i t h bigger contact angles are more likely to break up t h a n rivulets w i t h smaller contact angles, w h i c h makes sense. 27 Chapter 4 Combined Air and Rivulet Flow Two-phase flow c a n exist i n a variety of regimes, such as b u b b l y flow, slug flow, churn flow, a n d annular flow. T h i s is described i n detail by Fowler[13]. Fowler also i n t r o d u c e d a flow regime map, w h i c h classifies the flow patterns according to density, velocity, a n d gas volume fractions. However, this classification doesn't include the case of the two-phase flow i n a t y p i c a l s m a l l dimensional fuel cell. Those fuel cells are i n the order of 1 m m diameter, a n d the velocity of flow is t y p i c a l l y i n the order of l m / s for air a n d 1 0 m / s for water. So we propose a new m o d e l to _ 4 study the flow inside a fuel cell. We m o d e l the fuel cell channel as a circular pipe w i t h a given diameter a n d s t u d y the flow i n 2-D. W e also assume that the flow is i n rivulet form. F o r a s m a l l r i v u l e t w i t h a s m a l l c a p i l l a r y number (equals to 0.14 for L = 1 m m ) , pressure and surface tension are the dominant forces that affect the m o t i o n of a r i v u l e t . G r a v i t y can be ignored to simplify the c o m p u t a t i o n . F i r s t , we do the direct p r o b l e m , i.e. compute the fluxes of the coupled w a t e r / a i r flow for a given pressure gradient a n d configuration. T h e n we do the inverse p r o b l e m , where fluxes are given, a n d we compute the pressure gradient needed to drive the flow a n d configuration of the r i v u l e t . T h e n we a p p l y the p r o b l e m to a specific situation, m o d e l i n g water movement inside a fuel cell. M a n y different aspects of research have been done to understand the science b e h i n d a p r o t o n exchange m e m b r a n e ( P E M ) fuel cell. S u c h as gas flow i n the electrodes, gas flow w i t h condensation a n d other t h e r m a l effects, and l i q u i d (water) m o t i o n d r i v e n by air flow. O n e aspect that affects the performance of P E M fuel cells is the a c c u m u l a t i o n of l i q u i d water i n the oxygen channel. Presence of certain amount of water is necessary for ionic t r a n s p o r t a t i o n i n the 28 membrane. However, i f there is too m u c h water inside the electrodes, it could block the flow of oxygen, thus reducing the efficiency of the fuel cell. I n some fuel cell designs, l i q u i d water is pushed out into the channels w i t h the flow of reactant gases. T h u s , a n u n d e r s t a n d i n g of the two phase (water a n d air) flow is essential to improve the performance of fuel cell. T h e geometry of a n idealized flow channel is shown i n F i g u r e 4.1. F o l l o w i n g the same as- F i g u r e 4.1: Geometry of a fuel cell channel. sumptions listed i n chapter 1 a n d neglect gravity i n a d d i t i o n , the system is governed by u = (0,0,w(x,y)) w i t h p(z) = Xz . E q u a t i o n (1.1) becomes fj,Aw - A = 0 (4.1) Physically, the pressure drop along x - y d i r e c t i o n is governed b y e q u a t i o n (1.9). Since G 2 = Gz = 0, we get Pw-Pa = -crk (4.2) where a is the surface tension of water a n d h is the height of the r i v u l e t . P is the air pressure a and P w is water pressure. Since gravity is ignored, the pressure d r o p i n the x - y plane is constant. We get a constant curvature i n the plane. Constant curvature corresponds to the geometry of circular arcs. C i r c u l a r arcs exist i n two states. O n e is the full circle w h i c h corresponds to annulus flow, a n d the other geometry is a n arc of a circle w h i c h corresponds to circular arc flow. These two k i n d s of flows are investigated i n d e t a i l below. B u t first we s t u d y the circular channel flow w i t h only a i r i n i t . 29 4.1 A i r Only Channel Flow A s s u m e there is no water present i n the circular channel, we solve the poission equation p Aw a = A i n polar-coordinates w i t h b o u n d a r y c o n d i t i o n w = 0 . L e t the radius of the channel be R. T h e n the velocity of air is: W (r) = -^(-r a 2 + R) 2 (4.3) To get the flux of air, we integrate the veloxity over the surface area, i.e. R f Q = 2irp a a Jo W (r)-rdr a (4.4) w h i c h gives us the following linear relationship between the pressure gradient a n d flux of air. A= -fjr-Qa (4-5) irR^Pa W h i c h agrees w i t h the result of [6]. N e x t we allow b o t h water a n d air coexist i n the channel, a n d do a s i m i l a r analysis. 4.2 Annulus Rivulet F i r s t , assume the rivulet clings itself to the w a l l of the channel s y m m e t r i c a l l y , w i t h air d o w n i n the core as shown i n F i g u r e 4.2. 30 flowing B o t h the water a n d air domains are governed b y similar equations, a n d b y the same b o u n d a r y conditions as for the air only case. W e have p Aw = A i n the air region, a n d p Aw a w = A i n the water region. Together w i t h the b o u n d a r y a n d interface conditions below Boundary Condition: w = 0 Interface C o n d i t i o n : [w] = 0 a n d [fj,w ] = 0 n we compute the velocity profile inside the channel. T h e interface conditions correspond to continuity of velocity a n d stress at the water-air interface 7 . [•] denotes the difference i n the quantities across the interface 7 . A n d w is the n o r m a l velocity. n For a given channel size R, a n d radius of air flow r , we find the velocity of air a n d water at radius r are: W {r) a = A A - — ( f - R ) + — (r - f ) 4A*a W (r) = -^-(r -R ) w 2 2 2 2 2 (0 < r < r) (f<r<R) 2 Integrating velocity over the area of each corresponding d o m a i n , we get the fluxes of water a n d air. Q w = j W (r)27rrdr = TT/O^A w -R -r 4 Qa = f Jo =-- 7rp A 1 8 R (R -r ) 4 2 2 2 Mi W {r)2Trrdr a a ,4 8/J, j_ a 4 \ \H W pJ a We are interested i n the relationship between pressure gradient a n d water flux for fixed air flux. T o find i t , we w i l l compute the fluxes u p to leading order for s m a l l h = R — f a n d fixed gas 31 flux. f Q w = R / W {r)2nrdr w JR-h = -¥l± * (h ) + 3 h + 0 - rR-h Qa = W (r)2jrrdr a Jo 7TA r . R^ R^ R^ . _ ., 9. T , Let C\ = — T^- a n d Co = Jr- + £ 7 - , write A = A + A' where A is the base solution, a n d A' is 0 e a s m a l l p e r t u b a t i o n from the base solution. T h e n a linear analysis for s m a l l h gives: Qa = l(X + \')(C +C h = ^(X C 0 1 0 1 + ---) + \'C 2 1 + XCh 0 2 + ---) For fixed air flux, the last two terms must cancel. T h u s we have A' OC h. Since Q oc h , we get 2 w the pressure gradient increases as the square root of the increase i n water flux, i.e. A' 4.3 OC (4.6) Circular Arc Rivulet For the case that water accumulates at the b o t o m of the channel, we assume it forms a rivulet w i t h circular arc interface there, shown i n F i g u r e 4.3. 6 is the contact angle, a n d ip is the angle that determines the size of the circular arc. T h e governing equations are the same as the annulus rivulet case, except that the geometry is different. Since there is no a n a l y t i c equations to describe the flow field w as a function of 9 and ip, we are not able to solve the coupled p r o b l e m analytically. Instead, we used a software package called F E M L A B to calculate the velocity profile. F E M L A B allows the user to specify the geometry a n d the P D E s to be solved, together w i t h b o u n d a r y a n d interface conditions. T h e n it generates triangular meshes a n d uses finite element m e t h o d to solved the P D E . F i g u r e 4.4 below shows the meshes generated by F E M L A B after refinement. 32 F i g u r e 4.3: Cross section of a circular arc rivulet channel. 4.3.1 T h e direct problem F i r s t we solve the direct problem, i.e. w i t h ip a n d A given, find the size of the rivulet, compute the velocities of water a n d air, a n d then integrate to get the fluxes. T h e following d a t a are used i n the calculations. • fj, = 1.8 * I O - 5 a • H w = 1.0 * 1 0 - 3 kg/m/s kg/m/s • R a d i u s of channel: 1 m m • Pressure gradient: 100 P a / m • C o n t a c t angle (9): 84° (between graphite a n d water)[10] T h e velocity contours are shown i n figure 4.5. N o t i c e that the velocity contours have different slopes at the interface due to continuity of shear stress a n d the difference i n viscosity of air a n d water. T h e n u m e r i c a l solutions given by F E M L A B are verified using a simpler rectangular geometry where a n a l y t i c a l solutions can be obtained easily, (see A p p e n d i x for details) 33 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 F i g u r e 4.4: Triangular mesh of the cross section domain 1.031 0.9792 0.9276 0.8761 0.8246 0.773 0.7215 HO-67 0.6184 0.5669 0.5154 n B 0.4633 0.4123 0.3607 0.3092 0.2577 0.2061 0.1546 0.1031 0.0515 F i g u r e 4.5: Contour plot of fluxes 4.3.2 The inverse problem T h e inverse p r o b l e m is: given fluxes of a i r ( Q ) , water (Q ) a n d size of the channel a w determine the size of the rivulet (ip) a n d the pressure gradient (A) needed to drive the flow. Consider the general equation fiAw = A. L e t x —> L £ , y —> L y , where x a n d y are nondimensional quantities. Thus ^ dw dx 2 2 dw dy 1 2 L 2 2 dw dx dw dy 2 2 2 2 Since Q ~ f wdA, for fixed Q , we have A varies like 1 / r . 4 34 ( A factor of 1 / r 2 comes from w a n d another factor of 1 / r 2 comes from A , the area of the cross section). F o r c o m p u t a t i o n a l simplicity, let R = 1, and also A = 1. W e w i l l use A = 1 a n d R = 1 i n our proceeding calculations of the reference water a n d air fluxes. N o t e that the equation fj,Aw = A is linear i n A. T h i s property is e x p l o i t e d below to reach a scalar inverse p r o b l e m . T h e c o m p u t a t i o n is carried out i n the following steps: 1. B a s e d o n the r a d i u s ( R ) of the channel a n d contact angle(#) between water a n d the channel, calculate the size of the rivulet that gives Q tio ra = Q using A = 1 a n d R = 1 by m o d i f y i n g 2- the angle ip i n successive iterations, a n d terminate w h e n the specified error tolerance is reached. 2. Use either equation A = 3Q!°A=I) ' = 3QI°A=I) OR A t o c a l c u l a t e t h e pressure gradient a n d fluxes for the specified channel size and contact angle. N e x t we a p p l y the m e t h o d above to a specific case. Suppose we have a 1 m m diameter circular channel w i t h length l m . T h e n the effective area of the channel is diameter times length, w h i c h 1 0 c m . A l s o suppose the current density is lA/cm , 2 then the t o t a l current goes t h r o u g h the 2 channel is lOC/s. T h e electro-chemistry inside the fuel cell is governed by: H -> 2H + 2 0 2 + 2e~ + 4e~ + 45+ -» 2H 0 2 T h u s 1 mole of water corresponds to 2 moles of electron. T h e current I corresponds to moles of water, where F is the Faraday's constant is: Q w = 18 * 10/2Fkg/s, (9.632 * 10 C/mol). A a n d 18 is the molar mass of water. (32 + .78/.21 * 28) * 2 * 10/4Fkg/s, I/2F T h e flux of water T h e flux of air is: Q a = where 32 is the molar mass of oxygen a n d 28 is the molar mass of nitrogen, a n d O x y g e n occupies 21% of the volume i n air a n d nitrogen occupies 78%. O t h e r s m a l l percentage of gases that makes the air c o m p o s i t i o n are neglected. 2 is the usual stoichiometric flow factor w h i c h is the ratio of the amount of i n p u t air to the m i n i m u m amount of air needed for reaction. So the targeting fluxes of air a n d water are: Q = 9.3 * 10- kg/s A w 35 Table 4.1: Pressure gradient A ( P a / m ) a n d size of rivulet for different channel sizes e R 1mm 0.5mm 0.25mm 60° A = 2.83248e + 01 A = 4.529674e + 02 A = 7.2475e + 03 17.5° 84° A = 2.82130e + 01 A = 4.506544e + 02 A = 7.2105e + 03 12° 100° A = 2.82906e + 01 A = 4.526498e + 02 A = 7.2424e + 03 10° Q = 7.1 * lO'^kg/s a U s i n g these data, we find that the ratio of flux of air a n d water is: Q tio ra = = -6. 7 F E M L A B is used to carry out the steps 1-3 above, a n d compute the results i n T a b l e 4.1. T h e d a t a shows that A (pressure drop) is more sensitive to change i n size of channel(radius) t h a n to change i n contact angle. 36 Figures 4.6 a n d 4.7 show rivulet sizes for different specifications of contact angle and radius of channel. T h e r i v u l e t is concave i n F i g u r e 4.6, and convex i n F i g u r e 4.7. T h e qualitative shape of a rivulet depends o n a number of factors. F o r example, the contact angle between water a n d the channel, radius of the channel a n d the flux ratio. Notice that the velocity contours of water are not shown i n F i g u r e 4.7. T h i s is because the velocity of water is m u c h smaller t h a n the velocity of air. « 10" F i g u r e 4.6: shape of rivulet for r = l m m , 8 = 100°. F i g u r e 4.7: shape of rivulet for r=.25mm, 6 = 60°. In a n a c t u a l fuel cell, fluxes vary d o w n the channel. It w i l l be more useful from a p r a c t i c a l point of view that we have a way to find the pressure gradient for a n a r b i t r a r y Q tio- So we ra get the following universal curves. F i g u r e 4.8 gives rivulet size ip for any Qratio, a n d F i g u r e 4.9 gives the air flux (Q*) for A = 1. W i t h this information, we c a n estimate the pressure gradient for a n a r b i t r a r y given air flux Q a and radius r of the channel, i.e. Engineers c o u l d use these universal curves to recompute each specific configuration d o w n each unit section of the channel a n d integrate the local pressure gradient to o b t a i n a n a p p r o x i m a t i o n of the net pressure drop d o w n the channel. 37 0 F i g u r e 4.8: Rivulet size tp s- Qratio for R = l m m , 9 - 84°. y 4.4 10 20 30 40 50 70 80 90 F i g u r e 4.9: A i r flux vs. tp for R = l m m , 9 = 84°, and A = 1. Pressure gradient vs. flux of water We know that A a C\^Q W up to leading order for relatively s m a l l Q w with Q a fixed for annulus rivulet from asymptotic analysis, as indicated i n equation (4.6). However, it is impossible to get asymptotic results for the circular arc rivulet due to algebraic c o m p l e x i t y of the problem. So we d i d a comparision of pressure gradient vs. flux of water of the two cases numerically, and the results are shown i n F i g u r e 4.10. - O - arc rivulet • ©• annulus rivulet F i g u r e 4.10: Pressure gradient vs. Qw for annulus and circular arc rivulet We conclude that the pressure gradient of the annulus m o d e l is more sensitive to change i n 38 water flux t h a n that of the circular arc m o d e l , b u t the qualitative behaviours are similar. 39 Chapter 5 Stability Analysis of Combined A i r and Rivulet Flow 5.1 Stability of the annulus model T h e stability of a steady, a x i s y m m e t r i c , laminar, nondiffusive, p r i m a r y flow composed of two fluids flowing concentrically i n a straight circular tube has already been investigated by Hickox[8]. T h e s t a b i l i t y of the flow is determined by the m e t h o d of s m a l l p e r t u r b a t i o n s . B o t h a x i s y m metric a n d a s y m m e t r i c disturbances to the p r i m a r y flow are considered. O n l y infinitesimal disturbances to the p r i m a r y flow are considered; hence, the analysis is a linear stability analysis. H e found that the parallel flow is always unstable i n the sense that a s m a l l p e r t u r b a t i o n to the flow w i l l i n i t i a l l y grow w i t h time at an exponential rate, regardless of the size of the R e y n o l d s number. It was found that the most i m p o r t a n t single cause of i n s t a b i l i t y is the difference i n viscosity between the fluid regions. Density v a r i a t i o n a n d surface tension have less pronounced effects. However, only long wave numbers are considered i n his analysis. 5.2 Linear stability of a rectangular geometry There has been no study of the stability of our arc rivulets. A full scale analysis of the circular m o d e l is quite complicated and numerically involved. W e w i l l formulate a linear stability analysis of the free b o u n d a r y coupled flow i n a simple two d i m e n s i o n a l rectangular d o m a i n instead to get the flavour of the technique, w h i c h c o u l d be a p p l i e d to the c o m p u t a t i o n a l stability analysis of the arc rivulet. 40 Consider water flows down the unbounded channel in x-direction in region 1 with viscosity p , w and air flows in region 2 with viscosity p . a u=0 y=l : y=ho y=0 u=0 F i g u r e 5.1: rectangular geometry Steady state Given the pressure gradient A and the height of water h , we can calculate the velocity profile 0 in the two regions for a steady state flow. The flow is governed by steady state equation: fiV u 2 = V p = (-A,0), 0 o (5.1) where u = (u (y),0). 0 0 Integrate (5.1) with u = 0 on y = 0,1. (5.2) = (5.3) 0 then uoi Ay - —y , 2 P-w u 02 = B(y-l) + —(y-l) , 2 (5.4) Ma (5.5) po = —Xx. where A and B are determined by uoi(h ) 0 = u 2(h ), 0 — (/i ) == .(J-a-Q^Hw—z.— (ho) Ma~5 0 41 (5.6) 0 ("o)- (5.7) Solve the system a n d get A !±B, = (5.8) Pw B = C = — , a (5.9) A ^ Pa Pw ( 5 1 0 ) ^Pa Linearized perturbation Introduce a s m a l l p e r t u r b a t i o n to the steady state solution by w r i t i n g u = uo + u[, where = ( i, i), u v a n d p = p + Pi • T h e n substitute into the Navier-Stoke equation (1.1), ignoring 0 external forces (including gravity), i.e. p ^ + (u-V)u = pV u-Vp, (5.11) 2 (so where (u • V)u = (uu x + vu , uv y x + vv ). (5-12) y Substitute a n d remove higher order terms, get the linearized system for u\ u ui 0 x + VlU 0y uoviy \ ^ f V ui \ / \ VV \ pity 2 J J p ix (5.13) Dispersion relation W e are interested i n how the s m a l l disturbance v\ grows i n time. E x p e r i e n c e w i t h the m e t h o d of separation of variables a n d Laplace transforms suggests that the general solution can be expressed i n terms of: ii! = ui(v)e hi = hi{y)e (5.15) Pi = Pi(y)e (5.16) lkx+at ikx+at lkx+at (5.14) W e are interested i n k n o w i n g the sign of a i n solving these equation. If a > 0, the disturbance w i l l be amplified, growing exponentially w i t h t i m e u n t i l it is so large that nonlinearity becomes 42 significant. T h e n nonlinear s t a b i l i t y analysis w i l l be compulsory. If a = 0, it is said to be n e u t r a l l y stable. A n d i f a < 0, it is said to be stable or a s y m p t o t i c a l l y stable since the disturbance w i l l r e m a i n s m a l l for a l l time. Since a = a(k), where k is called the mode of the disturbance wave or the wave number, a n d it is generally true that a s m a l l disturbance of the basic flow w i l l excite a l l modes, i f a < 0 for at least one mode t h e n the flow is unstable. T h i s analysis follows the m e t h o d of n o r m a l modes described i n D r a z i n [7]: Analysis Substitute 5.14-5.16 into (5.13), get ikui +v[ •= 0 (5.17) opu\ + ikuoUi + V\UQ = p,(—k u\ apvi + uov[ = n(-k vi — ikpi 2 (5.18) + v'i) — p\. 2 (5.19) w i t h u\ = v\ = 0 on y = 0 , 1 . Specifically, un,u\2,«n, v\2 s h o u l d satisfy the following linearized b o u n d a r y a n d interface con- ditions. Uij denotes velocity component u i n region j(j = 1,2) w i t h i = 0,1 i n d i c a t i n g the base solution a n d the p e r t u r b e d s o l u t i o n . 1. N o slip c o n d i t i o n at the b o u n d a r y : •3(0) = 0=>un(0) =«n(0) =0 (5.20) u(l) = 0 => « i ( l ) = v {l) = 0 (5.21) 2 12 2. C o n t i n u i t y of velocity at interface: u[(h) = v {ho) = v {h ) u U2(h) = 12 0 (5.22) { u (ho) + h-iu'^iho) = ui {h ) + n 2 hiu' (h ) 0 Q2 0 3. C o n t i n u i t y of shear stress: ^ Hi{u' {h ) n 0 ( h ) = + hiu' \(h )) 0 0 43 = m ^ ( h /j, {u (h ) 2 12 0 ) ( 5 + hiu {h )) 02 0 . 2 3 ) (5.24) Derivation: n= « i + (/i') + 2 So the linear t e r m for f a n d h is f«(l,/i') and n « ( - / i ' , l ) (5.25) W e only keep the linearized terms i n the following computations. ua,T T-un = d U t h T dn = (l,h') • (un,vn) (5.26) = (l,/i')-(« (5.27) = U01+U11 + - - - = n-V« 0 1 +«n,0-l-t;ii) (5.28) (5.29) i l l T I u 1 \ /dun dun - ^ x - ^ 1 = { h A T ) T { ] ( 5 - 3 0 ) = (-h\l)-(0,u' +u ) (5.31) = u (5.32) 0l 01 u +u n Thus, (ft) = u' (h ) + it'ii(/^o) + ^ i o i H (5-33) u dn 01 0 Similarly, (ft) = u' (h ) + u' {h ) + hiu dn 02 0 12 0 Q2 +• • ' (- ) 5 34 T h u s we get equation (5.24) since Mi%f- = M2%^f at / i = fto4. Balance of n o r m a l stress: CTK ^ g ^ ' " (ft) ~ Mi d g ^ ' " (ft) = (p ~ P i i ) + n v' {h ) - Miu'u(fto) = ( P 1 2 - P n ) - crfc fti(ft ) 2 l2 0 44 i 2 2 0 (5.35) (5.36) Derivation: K = (1 + ^ 2 ) 3 / 2 = h (l-\ = un,n = n-Uii -k hije° 2 = (-h = (-h',1) = du il,n dn ) xx , 1) • (uii,vn) • (UQI + u i i , 0 + u i i ) -h'u i+vn 0 = n • Vuji = (-h',1) • (0,-h"u i = -h'luoi -y'iu' i n - h'u' 0 Q1 +v' 0l n S i m i l a r l y , we get dUi2 n II I I = -h u -yiu t dn 1 0 2 0 2 I + V l2 P l u g (5.46) a n d (5.47) into (5.35) get (5.36). 5. K i n e m a t i c free surface condition: h = vi(h) - ah\ = v{h§) — ikh\Uo(ho) t h u(h) x 6. C o n t i n u i t y equation i n d o m a i n : ikun +vn ikui + Vi 2 2 45 = 0 = 0 +v' ) n Stream function T o simplify the system, we introduce a stream function ip(x, y, t) such that the two components of the disturbance velocity are given by u\ = dtp/dy and v\ = —dip/dx If we next let iP(x,y,t) = f>(yy ° kx+ t ( then ui = (p and vi = —ikcp (5.52) disregarding the exponential part. T h u s the continuity equation (5.17) is a u t o m a t i c a l l y satisfied, a n d equation (5.18) a n d (5.19) become: —iku (p + (op + iku 0 — k p)(p — pep = 2 0 — (ikop + ik p)<p — ikuo<p +ihp<p = 3 —ikpi (5.53) —p (5.54) 1 Solve pi from (5.53), then substitute into (5.54), we o b t a i n a fourth-order Orr-Sommerfeld equation: (u' ' + ikop + ik p)cp + (ikuQ — u' )(p' + (iapk 3 — uo — 2pk)cp —ik~ p<p -1 0 l 0 =0 . (5.55) N e x t we formulate the b o u n d a r y a n d interface conditions i n terms of (p. 1. N o - s l i p c o n d i t i o n : Mo) <i>'M = &(O) = 0' (O) 2 = = o (5-56) 0 (5.57) 2. C o n t i n u i t y of velocity: Mh) = 4>' + hi4>' \(h ) u 0 46 (5.58) Mho) 0 = 0 i + hi 2 (M (5-59) 3. C o n t i n u i t y of shear stress: m(</>ii(M + fri0o'i(M) = /-t2(0i2(M + hi(j}Q (ho)) (5.60) 2 4. Balance of n o r m a l stress: E q u a t i o n (5.53) gives Pi = u (f) + (i(op — k/j,) — u )(f> — ik~ p,(f) (5.61) l 0 0 Substitute into equation (5.36), get u ((t)i2-4>ii)+{i' p- o)<t>i2(ho)-(io-p-u )(j) (ho)-iiJ,2k~ <fr 2+ilJ.ik~ (/) -ak hi(h ) a u 1 0 0 u 1 1 =0 2 11 0 (5.62) 5. K i n e m a t i c free surface condition: ahi =—ik(f>i(ho) — ikhiuo(ho) (5.63) N o w we are ready to solve the eigenvalue p r o b l e m (5.55) n u m e r i c a l l y subject to b o u n d a r y a n d interface conditions (1 — 5) above. F i r s t we discretize (j> o n d o m a i n . L e t $ denote the discrete values. T h e discretized values i n c o r p o r a t i n g b o u n d a r y a n d interface conditions c a n then be w r i t t e n as: A a n d B are matrixes. T h i s is a generalized eigenvalue p r o b l e m w h i c h c a n be solved numerically. a is the eigenvalues of A. W e w i l l like to know i f any Re(a) > 0 for any wave numbers k. For the circular arc flow, the system is similar to the system above b u t w i t h m u c h more complexity since it is two dimensional. However, we s h o u l d also be able to formulate the corresponding eigenvalue p r o b l e m numerically. A real difficulty is how to handle perturbations to the contact line. M o v i n g contact line conditions are very p o o r l y u n d e r s t o o d i n literatures. 47 Chapter 6 Conclusions and Future Work T h e u n i d i r e c t i o n a l steady flow of a r i v u l e t driven by gravity a n d surface tension has been investigated b o t h a n a l y t i c a l l y a n d numerically. A n a l y t i c a l expressions for the cross-section, fluid velocity a n d flux have been obtained. T h e a p p r o x i m a t e s o l u t i o n gives good agreement w i t h the numerics for contact angles up to 3 0 ° . T h e r i v u l e t energy has also been considered to determine whether it is energetically favorable for r i v u l e t break-up. It is found that gravity d r i v e n rivulets are unstable to breakup by this mechanism for lengths of the order of millimeters for a range of contact angles but that shear d r i v e n rivulets are stable i n this way. A n interesting p o i n t is that the r e l a t i o n between fluid flux a n d r i v u l e t w i d t h is independent of surface tension i n the h i g h surface t e n s i o n / n a r r o w rivulet scenario. However, surface tension does appear i n the low surface t e n s i o n / w i d e rivulet result. T h i s s u r p r i s i n g result carries t h r o u g h to the non-zero inter-facial shear results. T h e a s y m p t o t i c study of r i v u l e t flow d o w n a n i n c l i n e d surface highlights a strong dependence on the l i q u i d - s o l i d contact angle shown i n the flux a n d energy vs. length plots. T h i s discrepancy may well arise from inaccurate measurements of the contact angle as i n d i c a t e d i n [1]. A n o t h e r source of error is that the a n a l y t i c a l results are s t r i c t l y only v a l i d for a s m a l l contact angle, hence a non-negligible correction error is i n t r o d u c e d for larger contact angles. T h e c o m b i n e d air a n d r i v u l e t steady flow d o w n a circular channel (ignore gravity) is studied using the annulus m o d e l a n d circular arc m o d e l . For the annulus r i v u l e t , the change i n pressure 48 gradient is p r o p o r t i o n a l to the square root of the v a r i a t i o n of water flux, for s m a l l water fluxes a n d fixed air flux. S i m i l a r result holds for the circular arc r i v u l e t m o d e l , b u t o n l y n u m e r i c a l results are presented due to algebraic c o m p l e x i t y of the p r o b l e m . For the circular arc m o d e l , we also investigated the inverse p r o b l e m , a n d found that pressure drop is more sensitive to change channel size t h a n to change i n contact angle w h i c h relates to the m a t e r i a l used to b u i l d the channel. F i n a l l y , linear s t a b i l i t y of the models are considered. W e found from literature that the annulus rivulet is unstable to b o t h a s y m m e t r i c a n d a x i s y m m e t r i c disturbances for long wave numbers. A l t h o u g h linear s t a b i l i t y analysis based on s m a l l p e r t u r b a t i o n theory predicts instability, more work s h o u l d be done to investigate the s t a b i l i t y p r o b l e m more thoroughly. T h e s t a b i l i t y of our c i r c u l a r arc rivulet flow has never been investigated before. O n e has to study the s t a b i l i t y of the c i r c u l a r arc r i v u l e t m o d e l numerically. Instead of s t u d i n g the stability of the c i r c u l a r arc r i v u l e t directly, I present the linearized s m a l l p e r t u r b a t i o n equations a n d interface/boundary conditions for a simpler rectangular geometry. T h e i m m e d i a t e next step w i l l be to discretize the d o m a i n , resulting i n a n eigenvalue p r o b l e m for the g r o w t h rates. M o r e work can be done to improve the a s y m p t o t i c analysis of r i v u l e t d o w n a n i n c l i n e d plane. W e could include the second order t e r m (e terms) i n the analysis, a n d expect that there w o u l d 2 be better agreements of the flux a n d energy vs. length plots between numerics a n d asymptotics. I n the end, the c o m b i n e d a i r / w a t e r flow m o d e l could be tested against e x p e r i m e n t a l d a t a to show its validity. 49 Appendix .1 Verification of femlab computation using a rectangular domain For a rectangular geometry shown below, we have Direchelet conditions o n b o t h the top a n d b o t o m , a n d N e u m a n b o u n d a r y conditions on the two sides. T h e r e are also continuity of velocity a n d stress at the interface. L e t p\ = 1 i n region 1 and ix = 2 i n region 2. After a p p l y i n g basic 2 u=0 region 1 [u]=0 o I ,|u.Un]=0 region 2 u=0 F i g u r e 1: Verification geometry, techniques, W e calculate the exact solution to the P D E . i.e. 2» l ^ 12 y 12 Integrate over the d o m a i n , we find that the flux is 0.05729. T h e n we used the P D E toolbox in F E M L A B to solve the same p r o b l e m numerically. A n d the flux we get is 0.05134 using 100 contour values. T h u s the relative error is less t h a n 0.1% . See figure 2 below for the solutions. Solution U vs. y analytical solution 1 numerical solution | i -''"""\ / / / / \ \ \ \ \ / x \ \ \\iX1 :\ numerical soluti 20 contour values / 0 0.1 0.2 0.3 0.4 0,5 y 0,6 0.7 0.8 ;\ ^ ; \ 0.9 1 F i g u r e 2: Analytical and numerical solutions 50 Bibliography [1] T . G . M y e r s Thin films with high surface tension, S i a m Review,40(3), 441-462, 1998. [2] T . G . M y e r s , M . H . X . L i a n g a n d B . W e t t o n , The stability and flow of a rivulet driven by interfacial shear and gravity, International J o u r n a l of N o n l i n e a r - M e c h a n i c s , s u b m i t t e d . [3] S h m u k i a n d Laso, On the stability of rivulet flow J F M 215, 125-143, 1990. [4] C h i a - S h u n Y i n , Wave formation o a liquid layer for de-icing airplane wingsJFM 212, 41-53, 1990. [5] T . G . M y e r s a n d C . P . T h o m p s o n , Modelling the flow of water on aircraft, A I A A J o u r n a l 36(6) 1010-1013, 1998. [6] G . K . Batchelor, An introduction to fluid dynamics, C a m b r i d g e U n i v e r s i t y press. [7] P . G . D r a z i n , W . H . R e i d Hydrodynamic stability, C a m b r i d g e U n i v e r s i t y press, 1981. [8] C . E . H i c k o x Instability due to viscosity and density stratification in axisymmetric pipe flow, Physics of F l u i d s , 14(2), 251-262, 1971. [9] D . J o s e p h , Y . R e n a r d y Fundamentals of Two-Fluid Dynamics, Interdisciplinary A p p l i e d M a t h e m a t i c s , V o l u m e 4, Springer-Verlag. [10] Handbook of physics and chemistry [11] M A T L A B / F E M L A B , T h e M a t h w o r k s Inc., http://www.mathworks.com. [12] E . B . D u s s a n V , On the spreading of liquids on solid surfaces: static and dynamic contact lines, A n n u a l R e v i e w of F l u i d Mechanics, 1 1 , 371-400, 1979. [13] A . C . F o w l e r Mathematical models in the applied sciences, C a m b r i d g e texts i n a p p l i e d m a t h ematics, 1997. 51
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Rivulet flow and stability
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Rivulet flow and stability Liang, Margaret Hongxia 2002
pdf
Page Metadata
Item Metadata
Title | Rivulet flow and stability |
Creator |
Liang, Margaret Hongxia |
Date Issued | 2002 |
Description | Two dimensional, steady-state rivulet flow down an inclined plane, and combined rivulet/air flow i n a circular pipe are studied. Both asymptotic analysis and numerical methods are used to solve the Navier-Stoke equations for rivulet flow down an inclined plane. Asymptotic approximation, which is valid only for rivulet of large cross-sectional aspect ratios, is substantiated by exact numerical solutions up to specified error tolerance. Standard shooting techniques are used to solve ODE for rivulet shape and finite element methods are used to approximate the flow velocity. The stability of the rivulet flow is also studied based on energy criterion. We find that pure gravity driven rivulets are subject to break up and pure shear driven rivulets are stable in that sense. For the combined air/rivulet flow, the relationship between pressure gradients and fluxes is investigated using the annulus rivulet and circular arc rivulet models. The annulus rivulet model is solved analytically, and the circular arc rivulet is computed using FEMLAB . It is found that pressure gradient needed to drive the rivulet flow is much more sensitive to change in channel size than to change in contact angle. Finally, linear stability analysis of the combined air/rivulet flow in a rectangular domain is formulated. |
Extent | 3549168 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-10-24 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080054 |
URI | http://hdl.handle.net/2429/14160 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2002-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2003-0263.pdf [ 3.38MB ]
- Metadata
- JSON: 831-1.0080054.json
- JSON-LD: 831-1.0080054-ld.json
- RDF/XML (Pretty): 831-1.0080054-rdf.xml
- RDF/JSON: 831-1.0080054-rdf.json
- Turtle: 831-1.0080054-turtle.txt
- N-Triples: 831-1.0080054-rdf-ntriples.txt
- Original Record: 831-1.0080054-source.json
- Full Text
- 831-1.0080054-fulltext.txt
- Citation
- 831-1.0080054.ris
Full Text
Cite
Citation Scheme:
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>
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.831.1-0080054/manifest