S L O W F L O W O F A B I N G H A M FLUID IN A G A P O F S L O W L Y V A R Y I N G W I D T H by DANIEL PADRAIGH RYAN B.Sc. (Mathematics) University of Minnesota, 2000 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES Institute of Applied Mathematics, Department of Mathematics We accept this thesis as conforming to^the^required standard THE UNIVERSITY OF BRITISH COLUMBIA September 2003 © Daniel Padraigh Ryan, 2003 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Abstract The slow Poiseuille flow of a Bingham fluid in an infinitely long gap of slowly varying width is considered. A consistent analytical theory is developed using a lubrication scaling without deviating from the Bingham model. It is shown that for sufficiently small gap width variations there is a rigid plug zone around the center of the gap. A numerical treatment of the problem is also presented using the method of augmented Lagrangian. ii Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgement vii Chapter 1. Introduction 1 Chapter 2. Equations and Scalings 3 2.1 Dimensional Equations 3 2.1.1 The Equations of Motion 3 2.1.2 Domain and Constraints 4 2.2 Non-Dimensionalization 5 2.3 The Lubrication Approximation 7 Chapter 3. Numerical Solution 9 3.1 Weak Variational Formulation 12 3.1.1 Derivation of the Functional to be Minimized 12 3.1.2 Determining A P 15 3.1.3 The Minimization Problem 15 3.2 Solution of Augmented Lagrangian 17 3.2.1 The velocity subproblem 18 3.2.2 The q sub-problem 20 Chapter 4. Numerical Results 21 4.1 Notes on Numerical Implementation 21 4.2 Error in a Straight Channel 22 4.3 Results for Varying Gap Width 24 Chapter 5. Asymptotic Solution 33 5.1 Geometric Assumptions and Simplifications 33 5.2 Lubrication Rescalings 34 5.3 Regular Expansion for the Yielded Region 35 5.4 The Zeroth Order Solution 37 5.4.1 Zeroth Order Results 40 5.5 Corrected Solution for yr < yy 44 5.5.1 Solution in the Transition Layer for yr < yy 45 5.5.2 Matching for yr < yy 47 5.5.3 Composite Solution: yr < Vy 50 iii 5.6 Corrected Solution for yy < yr 51 5.6.1 Composite Solution for yr > yy • 54 5.7 Corrected Pressure Gradient 55 5.8 Plug Speed Determination 56 5.8.1 Derivation of Method to Determine the Plug Speed 56 5.8.2 Results of Plug Speed Determination Method 59 Chapter 6. Plug Breaking Criteria 62 6.1 A Solution for the Broken Plug 62 6.2 Derivation of Plug Breaking Criteria 64 6.2.1 The Necessary Condition for Plug Breaking 66 6.2.2 Sufficient Condition for a Broken Plug 66 6.3 Discussion on Results of Breaking Criteria 67 Chapter 7. Summary 70 Bibliography 71 iv List of Tables 4.1 A table showing the relationship between grid resolution (k elements in each direction) and the L2 velocity error norm 23 4.2 A table illustrating the convergence of the augmented Lagrangian. The parameters B = 5, h = .05 and L = 4 were used to generate these results .' 24 List of Figures 1.1 A sketch of the flow 2 2.1 A sketch of the domain before scaling 5 3.1 A sketch of the half-gap under consideration 10 4.1 A plot of the logarithm of the element size = h vs. the logarithm of the L 2 velocity error norm 23 4.2 A two color plot of the plug region with B = 5,h = .085, L = 4 26 4.3 A two color plot of the plug region with B — 1, h = .085, L — 4 27 4.4 A two color plot of the plug region with B = 20, h = .085, L = 4 27 4.5 A two color plot of the plug region with B = 5, h = .05, L = 4 28 4.6 A two color plot of the plug region with B — 5,h = .12, L = 4 28 4.7 A two color plot of the plug region with B = 5, h = .085, L = 10 29 4.8 A two color plot of the plug region with B = 5,h = .085, L = 20 29 4.9 Cross-sections of the gap half-width with B = 5, h = .05, L = 4 31 4.10 Cross-sections of the gap half-width with B = 5,h = .085, L = 20 31 4.11 A colormap plot of the pressure field for B = 5, h = .05, L = 4 32 5.1 Plots of yy and \po,x\ for B = 1,5 and 20 with h = .05 41. 5.2 Plots of yy and \po,x\ for B = 1,5 and 20 with h = .15 41 5.3 Plots of uo(0,y),uo(j,y) and uo(^,y) for various amplitudes of yi(x) 42 5.4 A plot illustrating the relationship between the plug speed and the position of yx relative to yy 43 5.5 A plot of UQ and u corrected to 0(52) on a cross-section where yr < yy 49 5.6 A plot of uo and u corrected to 0(52) on a cross-section where yr > yy 54 5.7 A plot of po,x a n d the corrected pressure gradient for B = 5 and h = .05 55 5.8 A plot of the plug speed functional using B = 5 , h = .02 , 5 = .05 60 5.9 Plots of plug speeds for various amplitudes for B = 5 and B = 20 61 6.1 A sketch of the domain showing a broken plug 64 6.2 Plots of the sufficient breaking criteria with B = 5 67 6.3 Plots of the sufficient breaking criteria with B = 20 68 vi Acknowledgement This research was conducted under the supervision of Dr. Ian Frigaard. I would like to thank him for for his guidance, advice, encouragement, and extreme patience throughout. I would also like to thank NSERC and Schlumberger Cambridge Research for funding the research in this thesis, via collaborative research grant CRD 245434. The Pacific Institute for Mathematical Sciences is also thanked for their additional funding. vii C h a p t e r 1 Introduction A Bingham fluid is an incompressible visco-plastic yield stress fluid. It is characterized by the fact that when the stress is below the yield stress, the rate-of-strain is zero, and the fluid moves as a rigid solid. When the stress is above this yield value, the rate-of-strain is in a linear relationship with the stress and the fluid flows in a viscous manner. Here, we consider the 2-D Poiseuille flow of a Bingham fluid in a gap of slowly varying width. The gap is assumed to be symmetric about the centerline and to have periodic width in the flow direction. This results in a flow with zero shear stress on the centerline. As a result, if the variation of the gap width is kept small enough, the stress on the centerline will be below the yield value, and a solid plug will exist along the infinite length of the gap. According to the Bingham model, this solid plug will move with uniform speed. See Fig. 1.1 for a sketch of the flow and domain. It is known that when a lubrication scaling is used with a Bingham fluid, there is an inconsistency at leading order. This inconsistency is that the lubrication model shows a plug existing, yet shows a non-uniform speed for the plug. This is not allowed by the Bingham model which requires the plug to behave like a rigid solid. In some instances people have shown that this is because there is not actually a plug region because the stress is, in fact, slightly above the yield value. Others have proposed alternate constitutive models for the fluid in order to resolve this inconsistency. Some such models use methods of viscosity regularization in which the solid is replaced by a fluid of large viscosity (e.g. [9] & [10]) and the visco-elastic model in which the rigid solid is replaced with an elastic. 1 r-r + > Ty, (2.8) 7 7 = 0 f < f y , where It should be noted that it is not possible to explicitly express the deviatoric stress in terms of the rate-of-strain when in a region where the stress is below the yield value, fy. The areas where f < TV have a zero rate-of-strain, hence they move like a rigid solid, and will be hereinafter referred to as "plug regions". 2.1.2 Domain and Constraints Fig. 2.1 is a sketch of the domain. In order to simplify the problem, it is assumed that the gap walls are symmetric about the centerline, y = 0. It is also assumed that the gap walls are periodic in the x spatial dimension with periodicity L. Because of the symmetry about y = 0, we can restrict our attention to only positive y. We will assume that there is a fixed average flow rate, UQ, being pumped into the gap infinitely downstream. This can be expressed using the mean flow rate UQ and mean half-gap width D as / udy = U0D V f , (2.10) where yi(x) is the local half-width of the slot. There is a no-slip condition on the velocity field at ±yi(x), the channel walls. This means u = v = 0 on y = yi(x). (2-H) 4 Chapter 2. Equations and Scalings Figure 2.1: A sketch of the domain before scaling. As a result of the symmetry of the flow about y = 0, we know there is no shear stress on the centerline, i.e. rxy = 0 on y = 0. (2.12) 2.2 Non-Dimensionalization We now use the mean half-gap width, D, to scale the space variables, and the mean flow rate, UQ, to scale the velocity components. X = Dx y = Dy u = UQU V = UQV i b ~ U0 p = jlUo -f7p Ixy = Uo. 15lxy TXy (IUQ D lyx = UQ. blyx Tyx p,U0 - b ryx 7xx — Uo. blxx fiUo - ^ Txx lyy = Uo. b l m Tyy fiUo - b rvy 7 UQ. ^ b 1 T = Jj-Uo b T (2.13) 5 du di du dx du. + v—) dy = -dp dx d dx d dy dv dv dx dv. dy dp dy d dx~Txy + d di dy du dv = 0. i | dx dy Chapter 2. Equations and Scalings Substituting these non-dirnensional variables into the momentum equations yield Txy (2.14) Tyy (2.15) (2.16). In (2.14) and (2.15), Re is the Reynold's number, a dimensionless quantity measuring the ratio of inertial effects to viscous effects which is given by Re = (2.17) The non-dimensional rate-of-strain tensor is now, ixx = 2 ^ , (2.19) (2.20) By introducing Bingham number,B, our constitutive relations become du dv dy dx 2— dx' 2— dy' T « = (1 + ? ) 7 « r>B (2.21) 7 = 0 <=^ T < B. The Bingham number is a dimensionless parameter that represents a ratio between the yield strength of the fluid and the size of the viscous terms, and is given by B = ^ . (2.22) The condition of incompressibility gives jxx = -jyy, which in turn gives Txx = -rvy. This allows us to write the second invariants of the rate-of-strain tensor and the deviatoric stress tensor as 7 = (7^ + 7**)i (2-23) and r = (r2xy + r2xx)l (2-24) 6 Chapter 2. Equations and Scalings 2.3 The Lubrication Approximation The approach is numerical in Chapters 3 and 4. Afterwards, we use a lubrication approximation which we derive below to study the flow analytically. For analytical purposes, we make use of the fact that the gap has a slowly varying width along the length of the channel, so the characteristic length, L, is much larger than the mean half-width, D. Taking this into account, we choose to scale x and y differently for the analytic approach. This type of scaling is referred to as a lubrication scaling. x = Lx y = Dy V = —-—V L DUo 7xy lyx 7xx Txx — - Txx 7KK T = ——T D Substituting these rescaled variables into the momentum equations yields d_ T~xx ~l~ n TXy, (2.25) dx dy The dimensionless small parameter 5 is defined by (2.26) (2.27) (2.28) 7 Chapter 2. Equations and Scalings and Re is given by (2.17). S will play the roll of an asymptotically small parameter in the analysis to come. The stress/strain relationship, rate-of-strain tensor and the deviatoric stress tensor are still given by (2.18) - (2.22), but the definitions of their second invariants are now i=(i2xy + S2i2xx)K (2-29) and T = (T2XV + 5*T2J. (2.30) This scaling will let us use a perturbation expansion with 8 as the small parameter to analytically study the flow. 8 C h a p t e r 3 Numerical Solution We describe here the background to the algorithm used to generate the results in Chapter 4. Since the focus is on the effects of the geometry on the flow, and in particular on the formation of the plug region, we consider only a Stokes flow (i.e. the inertial terms are neglected in the momentum equations). Furthermore, we assume that the channel is periodic in the z-direction and consider only one period of the channel: x G [-L/2, L/2]. The flow domain is shown in, For. this chapter it is very convenient to write the velocity components in a subscripted manner, i.e. u = (ui,U2). The spatial components are also written using subscripts x = {x\,X2)- This notation is only used in this chapter. We begin the analysis with the classical formulation of the problem. We take the non-dimensional 2-D Stokes equations, and the Bingham constitutive model detailed in Chapter 2 which were Fig. 3.1. (3.1) (3.2) 0 = V u , (3-3) (3.4) Tij — dui duj dxj dxi' (3.5). 9 Chapter 3. Numerical Solution Figure 3.1: A sketch of the half-gap under consideration. Fi ,r2 ,r3 and T4 are the contours making up the boundary of the domain Q. The boundary conditions imposed on u and Ty ( u ) are the following: = o, T~xx = 0, on Ti = X = L/2, y G (0 , W (L /2)) , (3.6) Ul = 0, U2 ~-= 0, on r2 = x e [ -L /2 .L /2 ] , y = yi(x), (3.7) U2 = 0, Txx = 0, on r 3 = X = -L/2, ye (0,yi(-L/2)), (3.8) U2 = 0, TXy - 0 , on r4 = X G {-L/2,L/2], y = 0. (3.9) We note that the periodicity conditions on u at x = ±L/2 have been replaced by (3.6) & (3.8), and this requires some explanation. In fact, for a Newtonian fluid, conditions such as (3.6) & (3.8) result purely from symmetry and can also be imposed at any of the positions x = nL/2, n an integer, as we argue below as long as y%(x) has maxima at x = nL/2. First, (for either Bingham or Newtonian fluid), we assume that there exists a unique solution u to the classical problem with periodicity conditions. Then we may replace the boundary conditions on u by their negative and we find that ^reverse ^ is the unique solution to this reversed problem, (note that the pressure gradients are reversed). This is the classical result of Stokes flow, e.g. see [3]. Now, since the geometry of the domain is symmetric, we can simply reflect the domain about x = 0, impose the same boundary conditions and we find a new solution ureflect to this problem. However, since the geometry is exactly the 10 Chapter 3. Numerical Solution same as before, we know that ureflect _ u Equally, we can see that the problem for urejlect is exactly that for ureverse but with the geometry and boundary conditions reversed. Thus, we have: ux{x,y) = u[eflect(x,y) = -urrerse(-x,y) = Ul(-x,y), and u2(x,y) = ur2eflect(x,y) = u^erse(-x,y) = -u2(-x,y). From this we see that «2 (0 , y) = 0 and -^-ui(0, y) = 0. Now, simply by shifting the coordinates to x = nL/2, n an integer, we can repeat the above argument. For a Newtonian fluid this suffices, but for a Bingham fluid we need note that ux = 0 does not necessarily imply that TXX = 0, except where yielded. However, if we suppose for example that the model is regularized, (e.g. by a biviscosity model), then it is clear that TXX = 0. Now we note that the symmetry argument is still stronger (at least for a Newtonian fluid). Taking higher derivatives, we see that at any point x = nL/2, n an integer, we have that: — ' ^ + — dx^xv dxdyUl dx1 Extending this to the Bingham fluid by assuming that A r - o ± T - 0 dxxy~ 1 dyyy~ at x = nL/2, n an integer, we can evaluate the y-momentum equation and have that §?(±L/2,j/) = 0. (3.10) oy Thus, at the ends of the domain we will also assume that p is independent of y. This will be used below in deriving our variational formulation. 11 Chapter 3. Numerical Solution 3.1 Weak Variational Formulation 3.1.1 Der iva t ion of the Funct ional to be M i n i m i z e d We assume that there exists a unique solution u to the classical problem derived above, in a weak sense which we define below. The solution u £ V, where V is the closure of V with respect to the norm, i.e. 2 I V2 \V\\HHSI) dxj Here U denotes the flow domain and V is defined by: V = {v e [ ^ ( f i ) ] 2 : V • v = 0 e fi; w2 = 0 on r b f 3 , T4; v = 0 on r2, }. (3.11) In deriving the variational formulation from the classical formulation, we assume that u , v G V and whatever other regularity is implicit (e.g. for pressure and stress fields). The stress is defined as <7y = —pSij + Tij, (3.12) where <5y = 1 if i = j , otherwise . (3.18) J-L/2 J Now, we see that these remaining terms vanish due to the non-essential conditions in (3.6), (3.8) & (3.9) and due to the non-essential condition (3.15) on p at the ends of the channel section. Thus, returning to (3.17) and noting that functions in V are divergence free, we have f d AP f / Tij—-[vi -Ui]dx= — / (vi - wi) dx, (3.19) Jn uxj h Jn and using the symmetry of r y we rewrite (3.19) as f 1 AP f / o 7 y ( v - u ) r u d* = ~r / ( y i - u i ) d K - (3-20) Jn 2 L Jn We now derive an inequality by considering the first term in (3.20). We consider separately and pointwise, regions in which T > B and those in which r < B. First, using the definition of T y 13 Chapter 3. Numerical Solution for r > B gives / J r y ( u ) 7 i j ( v - u ) = ^ ( l + ^ ) I ^ ( „ ) ^ ( v - u ) 5 7 « ( u ) 7 « ( v ) = / J ^ « ( v - u , + B / n « ^ - , ( u ) . The Cauchy-Schwarz inequality tells us ^7ij (u)7y(v) < 7(u) 7(v). We now have for r(u) > £ : / ^ j ( u ) 7 « ( v -u)< f 57y(u)7tf(v -u) + B f 7(v) - 7(11). (3.21) For r < 5 we do not have an explicit expression for r v like we did for T > B. However, we shall see that inequality (3.21) remains true. If r < B, we know 7 y(u) = 0 and 7 ( 1 1 ) = 0, so 1 1 2 r y (u)7ij(v - u) = 2 T y ' ( u ) ^ ' ( v ) -Again, the Cauchy-Schwartz inequality can be used to get the result ^ T i j ( u ) 7 „ - ( v ) < T(U) 7(V). Therefore, for the case where T < B, we have / ^ « ( u ) 7 « ( v - u) < / r(u) 7(v) < B f 7(v). (3.22). Jnz ' Jn Jn Taking into consideration 7(11) = 0 for T < J3 we see that in fact, equation (3.21) holds for all values of r . Combining with (3.20), this leads to the following weak variational formulation of the problem: I ^7«(u)7«(v -u) + B f 7 ( v ) - 7(11) >^f- f vx-Ul. (3.23) The solution u is the function u £ V for which inequality (3.23) holds for all v G V. We assume that such a function can be found. It follows that u will also minimize the functional J (v ) = / ^72(v) + B f 7 ( v ) -^f- fvi. (3-24) Jn1 Jn L Jn 14 Chapter 3. Numerical Solution 3.1.2 Determining A P Note that the solution that we require is in fact the pair (u, A P ) . We have scaled the velocity with the mean velocity through the channel and therefore should have that at any point x. By differentiating this integral with respect to x and using the incompressibility condition, we see that there is no x variation in the above integral. Thus, our flow rate constraint above can be written as: - f uidx = 1, (3.25) L Jn which corresponds to the linear functional in (3.23). The question now is whether we can find a constant A P such that (3.25) is always satisfied. We can see this is the case by using (3.23). Denote by U A P the solution of (3.23) corresponding to pressure drop A P , and suppose we have A P 2 > A P i . It is clear that U A P x is a test function for U A P 2 a n d vice versa. Substituting these functions for v in (3.23) for the two solutions, and summing gives: 0 < / 7 2 ( U A P I - UAP 2 )dx < ( A p 2 AP i ) t _ U A P i i l d x . Jn L Jn Thus, we see that the integral in (3.25), effectively the flow rate, increases monotonically with the pressure drop. Indeed, provided that A P 2 > APj and UAP X 7 ^ U A P 2 , the increase is strictly monotone. It appears therefore that we can decouple the problem for (3.25), by first solving for fixed A P and then iterating with respect to A P to satisfy (3.25). 3.1.3 The Minimization Problem It is worth noting that the term in (3.24) containing 7(11) is non-differentiable, which causes problems for any conventional minimization. This term will require special treatment and here 15 Chapter 3. Numerical Solution we follow the approach taken by Glowinski [4]. First we relax the minimization problem, by introducing an auxiliary variable, q. This variable will have two components, the first of which will behave like 7xX(v), a n d the second component acting like ^ ( v ) . Consider the following two constraints: dvi dv2 1 1 q i = -di---dy-= 2 l x x { v ) ' r M (3-26) 9 2 = W v ) - ( 3 - 2 7 ) We then replace -y(v) and 7 2(v) with |q| and |q | 2 respectively. Then the minimization is taken to be over v e V and q £ L2(Q), but with (3.26) and (3.27) satisfied. Then the problem is mm V q J ( v , q ) = / i | q | 2 + 5 / |q| - ^ / «!• (3-28) JQ 2 JQ LJ JQ The next step is critical. In order to relate the auxiliary variable q to the velocity field, penalty terms are needed. The reason for doing so is that (3.26) and (3.27) are difficult constraints to satisfy directly in a numerical discretisation. We thus relax (3.26) and (3.27) but add linear and quadratic penalty terms for these constraints. The other difficulty numerically is the incompressibility condition. This we also remove from the problem via the introduction of linear and quadratic penalty terms for the incompressibility. Thus, we consider: m i n max c (l „ A P v q A , p £(v,q,A,p) = ^ j - | q | 2 + B |q | - — V l (3.29) dvi dv2 . ,dvi dv2 ,dvi dv2 2 , ,dvl , dv2 \2 +1 _ p V - v + ^ ( V - v ) 2 J dx . where A, q € H = L2(Q) x L2(fl) and p G L2(Q). k\ and k2 are constants. And now, v G / C = {we {R\n)}2 : w 2 = 0 o n r i , T 3 , r 4;w = 0onr 2 } . The functional £(v, q, X,p) is referred to as the augmented Lagrangian functional. The vari-ables Ai , \ 2 and p play the role of Lagrange multipliers for the constraints (3.26), (3.27) and 16 Chapter 3. Numerical Solution incompressibility, respectively. With k\ = k2 = 0 this functional is referred to as a Lagrangian. The Lagrangian is then augmented with the quadratic penalty terms for stability, leading to (3.30). Note that p will in fact also play the role of the modified pressure p introduced earlier. 3.2 Solution of Augmented Lagrangian A solution to the augmented Lagrangian is attained by an iterative procedure of Uzawa type. This algorithm is presented in [4] as A L G 2 in §6.3. The algorithm goes as follows, {q 0,* 1,? 1} GHxHx L2{9) given : then, { q n - \ Xn,pn) known, we define {q n , A " + 1 , p n + 1 } by: • £ ( u n , q n - 1 , A B , f ) < £ ( v , q n - 1 , A n , p n ) Vv € K,un e AC, (3.30) £ ( u ^ q ^ A ^ ^ ) < £ ( u n , z , A n ) p , l ) Vz£H,qn£H, (3.31) r)iin hi)n X^ = X?+ P A ^ + — (3-33) p»+! = p " - p 2 ( V - u n ) , (3.34) where pi,p2 > Q are relaxation parameters associated with k\ and k2 respectively. For the computations presented in Chapter 4, we use p\ = p2 = fci = k2 = 1, which seems to work well. Glowinski [4] shows that this algorithm has the following convergence properties: u n —> u strongly in 1C, q™ —> q strongly in H, yi+l _xn > g s t r o n g l y i n JJf pn+l_pn y Q S T R O N G L Y I N A™ bounded in H, f1 bounded in L2{fl). To simplify the solution of (3.30), in the implementation that produced the results in Chapter 4, the individual velocity components were solved for separately in a decoupled fashion. First 17 Chapter 3. Numerical Solution u™ was solved for keeping u2 _1 fixed, then u2 was solved for keeping w™ fixed. In essence, this simplification amounts to weakening the quadratic penalty term for the incompressibility condition, as a result, it needed more iterations to converge. However, the linear systems being solved are only half the size, which results in much faster iterations. It is also much simpler to implement. 3.2.1 The velocity subproblem The subscript notation for velocity and spatial directions are no longer useful. The notation that we use now will be the same notation used for the remainder of the paper. It is u = (u, v) and x = (x,y). This notation is useful because it allows us to write partial derivatives as We now consider problem (3.30) for the ^-component of velocity, to see what the differential form of this minimization is. Let t be a real number and h G K be a variation of the z-component of the velocity field. Then, if u is the solution to (3.30) we know that, subscripts. 0 = dCR(u + th,v,q,X,p) dt j -phx + k2(V • u)hx + ki[ux -vy- qi]hx + ki[uy + vx - q2]h, •y (3.35) where for ease of notation. a = (-p + k2(ux + vy) + h[ux -vy- qi] + Xi , h[uy +vx- q2] + A 2 ) . After applying Green's Theorem, this becomes 0 (3.36) 18 Chapter 3. Numerical Solution Thus, since h is any admissible variation, we know that, 0 = V . a + — A P = fciV2u + k2uxx + V • (k2vy -p- kiqi + \i , A 2 - hq2) + — A P = fciVV + k2uxx + V f + — , (3.37) everywhere in the interior domain of Q. In the final formulation of the finite element equations du for u, there will be the natural conditions — = 0 on the ends of fl and the centerline, y = 0. on We now use the same method to derive an explicit differential form for v. Begin by taking a variation h about the second velocity component, v, using homogeneous boundary conditions for h on all boundaries: d£R(u,v + th,q,\,p) 0 at = jy -phy + k2(V • u)hy - k\\ux - Vy - q{\hy + ki[uy + vx- q2]hx —\\hy + X2hx dx} = [ Vh-bdx, (3.38) Jn where b is given by b = (ki(uy + vx - q2) + \ 2 , -p+k2(ux+vy) - ki(ux - vy - 51) - A i ) . Again, we apply Green's Theorem. This time the homogeneous conditions on h leave 0 = / / i (V-b)dx. (3.39) Jn This is true for any h satisfying the boundary conditions, so we can conclude, 0 = k\V2v + k2vyy + V • (A 2 - hq2 , k2ux -p - \\ + hqi) = hS72v + k2vyy + V • g , (3.40) everywhere in the interior of fl. As mentioned earlier, (3.37) and (3.40) will be solved in a decoupled manner. First (3.37) will be solved using v"1-1 to evaluate f. Then (3.40) will be solved using un to evaluate g . 19 Chapter 3. Numerical Solution 3.2.2 The q sub-problem If we have some fixed, u n , A" and and we only allow q to vary, we see our minimization problem becomes t A q ) = \ I | q | 2 + B f | q | + £ / | q | 2 - / C • q , (3.41) Z JQ JQ 1 JQ JQ where C = (Ai + h(ux - vy) , A 2 + ki(uy + vx)). We will see in §1 of Chapter 4 that when we choose bilinear elements for u and v, that the proper choice of basis for Ai , A 2 and p are piecewise constant on each element. The updates (3.32)-(3.34) are then performed in an average sense over the element, which is equivalent to evaluating the velocity derivatives at the centroid of the element. Al l this implies that q\ and qi should also be piecewise constant on each element, and the velocity derivatives in C should be evaluated at the centroid of the element. With this scheme, the relationships (3.26)-(3.27) between q and jij hold in an average sense over each element, and are satisfied exactly at the centroid of each element. This means that the minimization for q discussed here can be done separately over each individual element, which makes the solution of (3.41) simple. It is obvious from looking at (3.41) that for any fixed length of q , C will be minimized when q points in the same direction as C, i.e. q = dC, where d is a non-negative real. Therefore, we can reduce this minimization to the following 1-D problem, (where it is implicit that we minimize on each element): d>0 £^(d) = ^ ^ | C | 2 d 2 + ( B | C | - | C | 2 ) d . (3.42) Immediately we see that if B > |C | , then d = 0 and q = 0. However, for the case B < |C | , we differentiate the integrand and set the result equal to zero to find 1 - B ' d = ^ n F > 0 ' 1 + h q = dC. (3.43) 20 C h a p t e r 4 Numerical Results 4.1 Notes on Numerical Implementation For implementation, we represent u and v with bilinear quadrilateral basis functions, i.e. u = ^UjNj and v = J2VJNJ where Nj is piecewise bilinear. To obtain the element equations, we multiply both sides of (3.37) by an arbitrary basis function, Nj, integrate over a single element and write u and v as a sum of basis functions, giving V / uj(k1V2NjNi + k2Nj,xxNi)dx=- [ V - f J V i d x - / dx, (4.1) where f = (k2vy — p — kiqi + Xi , \ 2 - k\q2). The form of this equation suggests that the appropriate basis functions for p, Ai , X2, q\ and q2 are constant on each element. This is because we want V • f ~ V 2 iV j , and Nj is bilinear. Now we use Green's Theorem on both sides to get the final form of the element matrix equations arising from the variation of u: I V UJ (fciVJV,- • VN + k2Nj XN x) dx = - / f • VNt dx + / ^ dx. (4.2) The natural condition — arising from symmetry is used to eliminate the boundary integral on arising from this use of Green's Theorem. The same procedure is used to arrive at the following element equation for v: f V VJ(fciVJVj • VJVi + k2NjtyNiiV) dx = - / g • ViVj dx, (4.3) 21 Chapter 4. Numerical Results where g = (A2 - hq2 , k2ux - p-Xi + hqi). As mentioned in Chapter 3, to have the correct numerical solution we need to choose a constant mean pressure drop, A P , that yields a unit flow rate across the gap half-width. To find this A P , we implement an outer-loop around the solution of the augmented Lagrangian, i.e. first we solve the augmented Lagrangian with an approximate A P " , and then we choose A P n + 1 according to the flow rate given by the velocity field (recall the flow rate is monotonically related to the mean pressure drop). This process is repeated until the resulting flow rate across the gap half-width is within an acceptable tolerance range. To ensure this procedure works well, the outer-loop should only occur once the augmented Lagrangian has iterated enough to give a good approximation of the converged flow rate. In practice, it only took 4 or 5 iterations of A P to get an acceptable solution. 4.2 Error in a Straight Channel When there is no variation in the gap width, we know the analytical solution to the problem (it is developed in Chapter 5). We use this for a test problem to see how the element size affects the 1? norm of the velocity error. We calculate on a computational domain that is a 1 unit by 1 unit square (this corresponds to a channel with periodic length 2 and gap half-width 1) and let k denote the number of elements in either direction (we use a uniform square mesh). It is worth noting that when there is no variation in gap width that the solution is fully converged after only two iterations of the augmented Lagrangian. 22 Chapter 4. Numerical Results k L2 (error) 10 .00581103 20 .00146366 40 .00040736 60 .000245451 80 .000131523 100 .0000631973 120 .0000401938 Table 4.1: A table showing the relationship between grid resolution (k elements in each direc-tion) and the L2 velocity error norm. In(element size) Figure 4.1: A plot of the logarithm of the element size = h vs. the logarithm of the L2 velocity error norm. The line shown is found by a least-squares approximation, and has slope « 1.93. This suggests that 1?(velocity error) ~ 0(h ) for the straight channel problem. 23 Chapter 4. Numerical Results n L2(un-un-1) L 2 ( A ? + 1 - A") L 2 ( A " + 1 - A?) L2(pn+1 -pn) 50 3.95292 e-4 .00524789 .00163119 .00737470 100 1.13832 e-4 .00159296 5.63465 e-4 .00218611 150 3.09011 e-5 5.15219 e-4 1.95149 e-4 6.75485 e-4 200 9.55098 e-6 2.11226 e-4 8.15549 e-5 2.58076 e-4 Table 4.2: A table illustrating the convergence of the augmented Lagrangian. The parameters B = 5, h = .05 and L = 4 were used to generate these results. 4.3 Results for Varying Gap Width Table 4.2 gives an idea of how well the augmented Lagrangian converges. We see that the L2 norms of {un-un-1), ( A ? + 1 - A ? ) , ( A ^ - A ^ a n d (pn+1-pn) all converge to zero as n —• oo. Because of how we update Ai , A 2 and p, this means that « + < - « # ) , and V - u converge to zero as n —> oo. Thus we see that u is indeed converging to the solution of the original problem. If the velocity components were solved in a coupled fashion (as shown in (3.30)) then it would surely take fewer iterations to get convergence of the augmented Lagrangian. However, by decoupling the velocity problem, the size of the linear system is halved, which results in a much faster solution of the system. Experimentally, it has been the case that a higher Bingham number results in slower convergence of the algorithm, i.e. more iterations are needed. Figs. 4.2 - 4.8 are two color plots of the plug regions for various choices of B, h and L. The grey regions indicate where the fluid is yielded and the black regions indicate the unyielded 24 Chapter 4. Numerical Results plug regions. Fig. 4.2 can be considered a base case. Then Figs. 4.3-4.4 vary B, Figs. 4.5-4.6 vary h, and Figs. 4.7-4.8 vary L. We see that in the base case, where B = 5, an unbroken plug exists along the length of the gap. When we leave h and L unchanged and choose B = 1 (see Fig. 4.3), we see that the plug has broken into two plug regions separated by a stretch of fully yielded fluid. This seems intuitive as lowering B amounts to reducing the effective yield strength of the fluid. When B is increased to 20 (see Fig. 4.4), we see that there is an unbroken plug and in fact the plug region has increased in size compared to the base case. Again this seems intuitive as the effective yield strength of the fluid has been increased. This implies that the plug is more likely to remain intact for larger values of B. Also, for larger values of B we expect the yield surface to move out further towards the wall. Now, to investigate the effects of varying the amplitude of the wall variation, we refer to Figs. 4.5 and 4.6. In the base case, h = .085 is used. In Fig. 4.5 we decrease this value to h = .05. We see that when compared to the base case, the yield surface takes a different shape. The yield surface in Fig. 4.5 has less of a dip around x = —1/4 when compared to the yield surface in the base case. This could be understood by considering the size of \TXX\. When the slope of the gap wall is increased, the size of \TXX\ should increase as well. So, when we decrease the size of h, we have decreased the slope of the gap wall, especially in the region around x = —1/4. As a result the plug has increased in size around this region. When we increase the size of h, we would be increase the size of \TXX\ in this region. Looking at Fig. 4.6 we see that this is what appears to have happened, as now the plug region has broken into two separate plugs. This means that |r | has now exceeded the yield value around the region x = —1/4. This implies that the plug is more likely to remain intact for lower values of h. Lastly we consider the effect of the parameter L. In the base case we have L = 4. In Figs. 4.7 and 4.8 we use L = 10 and L = 20 respectively. We see that in both cases the plug has now broken into two separate plug regions. It does not appear that there is much of a difference between Fig. 4.7 and Fig. 4.8, however, looking closely we see that the size of the plug regions 25 Chapter 4. Numerical Results i. o 0 . 8 0 . 6 0'. 4 0 . 2 -1 .75 -1 .50 -1 .25 -1 .00 -0.75 -0 .50 -0.25 0.00 X Figure 4.2: A two color plot of the plug region with B = 5, h = .085, L = 4. The plug region is shown in black and the yielded region is shaded gray. differ very slightly. The most we can really conclude from these runs is that the plug is more likely to break when L is a larger value. We will see in Chapters 5 and 6 that all of the general conclusions reached here agree with the analytical results (when the Bingham number is not large, as this causes the asymptotic results to break down). 26 Chapter 4. Numerical Results 0 . 8 0 . 6 0 . 4 0 . 2 , , , I I I -1 .75 -1 .50 -1.25 -1 .00 -0.75 -0 .50 -0.25 0.00 X Figure 4.3: A two color plot of the plug region with B = 1, h = .085, L = 4. The plug region is shown in black and the yielded region is shaded gray. 1 . 0 0 . 8 0 . 6 0 . 4 0.2 -1 .75 -1 .50 -1 .25 -1 .00 -0.75 -0 .50 -0 .25 0.00 X Figure 4.4: A two color plot of the plug region with B = 20, h — .085, L = 4. The plug region is shown in black and the yielded region is shaded gray. 27 Chapter 4. Numerical Results I I I 1 r—— 1 1 I - 1 . 7 5 - 1 . 5 0 - 1 . 2 5 - 1 . 0 0 - 0 . 7 5 - 0 . 5 0 - 0 . 2 5 0 . 0 0 X Figure 4.6: A two color plot of the plug region with B = 5,h= .12, L = 4. The plug region is shown in black and the yielded region is shaded gray. 28 Chapter 4. Numerical Results i . o Figure 4.7: A two color plot of the plug region with B = 5,h= .085, L = 10. The plug region is shown in black and the yielded region is shaded gray. 1 . 0 0 . 8 0 . 6 0 . 4 0.2 Figure 4.8: A two color plot of the plug region with B = 5, h = .085, L = 20. The plug region is shown in black and the yielded region is shaded gray. 29 Chapter 4. Numerical Results Taking cross-sections of the u\ velocity component at various points along the channel will let us see what is happening more clearly. Fig. 4.9 shows cross-sections at x = —L/2,x = —L/4, and x = 0 for the case B = 5, h = .05 and L = 4 (note the plug region for this set of parameters was shown in Fig. 4.5). We see that the plug is indeed moving at a uniform speed. Looking at Fig. 4.10 we see that the cross-section taken at x = —L/4 has no plug region, and the cross-sections taken at x = 0 and x = —L/2 both have plug regions, but they have different plug speeds. This is consistent with Fig. 4.8 which shows that the plug has broken into two separate plugs and there is a region around x = —L/4 where no plug exists. In the region where no plug exists, the solution appears to be only very slightly yielded up to a certain value, after which it takes on a fully yielded profile. In fact, we will derive an asymptotic solution that has this property in Chapter 6. We call this type of solution a "pseudo plug solution" which is a term taken from Balmforth and Craster [2]. However, unlike the problem in [2], our solution can contain both regions where a real plug exists and regions with pseudo plug behavior. Other problems that exhibit pseudo plug behavior can be .found in Mei & Yuhi [5], and Piau [6]. 30 Chapter 4. Numerical Results Figure 4.9: Cross-sections of the gap half-width. This is for B = 5, h = .05, L = 4, which yields the plug region shown in Fig. 4.5. Figure 4.10: Cross-sections of the gap half-width. This is for B = 5,h = .085, L = 20, which yields the plug regions shown in F ig. 4.8. 31 Chapter 4. Numerical Results Figure 4.11: A colormap plot of the pressure field for B = 5, h = .05, L = 4. Looking at Fig. 4.11, we see there is little change in pressure in the y direction. However, when looking closely, there is a small pressure change right on the interface between the plug region and the yielded flow (the plug region is show in Fig. 4.5). 32 C h a p t e r 5 Asymptotic Solution Is this chapter, we use the lubrication rescaling of the Navier-Stokes equations presented in Chapter 2 to build an asymptotic solution. The purpose of this analysis is to present a consistent solution for the Bingham model, i.e. a solution that exhibits a rigid plug region. Recall that the notation in this section for velocity is u = (it, v), which is not to be confused with the notation in Chapter 3. 5 . 1 Geometric Assumptions and Simplifications For simplicity, it is assumed that the gap is symmetric about y = 0. This means that we can restrict our attention to only the upper half of the gap, y > 0. By assuming the gap is symmetric about y = 0, we are given by symmetry that rxy = 0 and v = 0 on this line. Also for simplicity, we take a sinusoidal variation of the gap half-width, i.e. yt(x) = 1 — hcos(2irx). This assumption allows easy control over both the size of the variation, and the size of its derivative through a single parameter, h. However, we shall see that this analysis would remain true for any smooth wall variation with the size of both its amplitude and slope bounded by h. We will consider only one period of the channel in x. In this analysis we assume a slow, steady flow and a small h. This means that we have Re ~ 0(6) and both the amplitude and slope of the wall variation are 0(6) as well. As we shall see, this causes the inertial terms in the momentum equations to be no larger than 0(53) and since we only use terms up to 0(S2), they will not play a role in this analysis. The assumption of slow 33 Chapter 5. Asymptotic Solution flow is also essential to our method developed to determine the speed of the rigid plug which uses a result concerning slow visco-plastic flow presented by Prager [7]. 5.2 Lubrication Rescalings Recall from Chapter 2, that if we scale x with L, y with D, and introduce the small parameter 5 = 5, then the Navier-Stokes equations for steady flow become L , du du. dp P o d d ,_ ,. 5 R e ^ x + V d - y ) = - / x + 5 V x T - + ^ ( 5 ' 1 } „ o „ . dv dv. dp r 9 d r o d „. 6 R e { u d x + Vrf = ~i + + V W ' ( 5 " 2 ) ^ + ^ = 0. (5.3) dx dy Using the same rescaling on the constitutive laws for the Bingham fluid yields the following reduced constitutive laws: (5-4). (5-5) (5.6) i5. (5.7) The relation between the deviatoric stress and the rate-of-strain is now: i + ? ) - y y T>B, (5.8). 7 = 0 T < B, (5.9) where B is the Bingham number given by B = (5.10) fiUo and r = (r2xy + 62rl)h. (5.11). du ^2 dv Ixy — dy dx' 7xx = dx' dv 7OT = dy 7 = (~2 +62-2 V ixy ~ " ixx 34 Chapter 5. Asymptotic Solution We now will construct an asymptotic solution for this flow, under the assumption of a small 5. First we will construct a solution for the yielded region. Then we will attempt to match this expansion to the plug solution. Lastly we will consider the problem of how to determine the speed of the plug (which is not given by the initial analysis). 5.3 Regular Expansion for the Yielded Region We begin our general analysis of the flow described by the lubrication scaling by constructing an outer expansion in the fully yielded region. We take regular expansions in 5 for u, v and p: p = po + Spi + 62p2 + ... , u = UQ + dui + 62U2 + ... , V = Vo + Sv\ + 52V2 + ... . For ease of notation, partial derivatives will now be written as a subscript x or y. The subscripts from expansion indices and subscripts from partial derivatives will be separated by a comma. Putting these regular expansions into the equations for the rate-of-strain, (5.4) - (5.7), gives the following expansions in 6: Ixy = "0,2/ + Sui,y + S2(U2,V + VotX) + 0(63), ixx = 2uo,x + 62ultX + S22u2,x + 0{S3). 7 2 = uliV + 52uo,yUhy + 52[2uo,y{u2,y + v0tX) + u\y + 4ulj + 0{53) Taking the square root of both sides, we find 7 = Itio ,| (1 + ^ + , 2 [ 2 ^ K , + ^ ) + ^ + 4< x ] + (5.12) (5.13) (5.14) (5.15) Now we use a Taylor expansion in 5 to get an expansion for 7 -1 7-1 = \U0,y\ V U0,V U2,y + VQ,x Ul,y 2 m 0 , X UQ,y + 2ul 0,y u0,y '0,2/ — (l- 5^- - s2u°>y(U2'y + Vo' x^ + 2u^ x ~ U^ y + o(53 I U0,y I UQA (5.16) 35 Chapter 5. Asymptotic Solution Using our expressions for the rate-of-strain and (5.8), we now find an expansion for the shear stress, rxy = {uo,y + Bsign(u0,3/) + SuitV + S2 U2,y + V0,x ~ Bsigli(uo 1 + 52TXVI2 + 0(<53). 0(c53) (5.17) For Tr Tyy we shall only need the leading order term, 2 1 + B uo,x + 0(5). (5.18) Substitute the expansions for u, v, p, and r y into the Momentum Equations,(5.1) and (5.2), and balance orders of 5 to get the zeroth, first and second order problems. Recall that we have Re ~ 0{5) because we have assumed a slow flow. 0(1) : 0 = -po,x + uo,yy, 0= PQ,y (5.19) (5.20) 0(S) : 0 = ~Pl,x + Ultyy, 0 = Pl,y. (5-21) (5.22) d d 0{52) : 0 = -u 0uo,x - v0u0ty - p 2,x + Qy~Txy,2 + dx~Txx<° ' 0 = -P2,y + faTxyfi + Q^Tyyfi-(5.23) (5.24) Recall we have the channel half-width being Ui(x) = 1 — hcos(2irx) and we have assumed h ~ 0(6). We shall see in the following section that uo = uo(yi(x), y) and therefore if h ~ O(S), we have that ux ~ 0(5) and by integrating the incompressibility condition out from y = 0 we get vo — 0. This means the inertial terms in (5.23) should actually be in the 0(<53) problem. 36 Chapter 5. Asymptotic Solution We also get as a result of h ~ 0(5): J~T*2/>° ~ 0(5), Txxfi — 0, Tyyfl = 0, TXy,2 = U-2,y Therefore, from the assumption of h ~ 0(5), we get the following problem at second order: 0(62): 0 = -P2,x + u2,yy, (5.25) 0 = p2,y. (5.26) 5.4 The Zero th Order Solut ion The leading order equations give: Po,y = 0 = > p o =Po(x), (5.27) Po,x = J~( u o,i/ + sign(u0,Y)B) <=> T>B, (5.28) 14(3,2/ = 0 T < B. To leading order, TO = | T X 2 / ] O | . The continuity equation at leading order is «o, i + v0,y = 0. (5.29) The leading order solution should also satisfy the imposed flow constraint, i.e. 10 mix) / u0dy = l Vx. (5.30) JoIn the region where r > B, we have at leading order P0,X = g^Txyfi- (5.31) 37 Chapter 5. Asymptotic Solution Remembering that there is no shear stress on y = 0 and, po = po(x), we integrate (5.31) with respect to y out from the centerline to get Po,xy = u0tV + sign(u0,j/)£. (5.32) So, for T > B, uo,y = Po,xy - sign(u0,jy)B. (5.33) Now, when r = B there is a yield surface, call it yy. In the unyielded region, uo,y = 0, so for continuity of derivatives, and hence continuity of shear stress, we need uotV —> 0 as y —» yv from the yielded region. Using this in (5.33) tells us about the position of the yield surface. (Note: This "yield surface" is only the leading order yield surface, and the actual yield surface may be different, or may not exist at all.) B yy = sign(u0,y) • (5-34) P0,x But, we know 0 < yy < y%, so, • yy = ^ , • (5.35) where r(yy)o = B, and sign(uo>y) = sign(po.x)- This means (5.33) can be written as uo,y = Po,x(v ~ Vy)- (5-36) Recalling that u(yi) = 0, we again integrate with respect to y, this time in from yi. ry ry / u0iydy= / PQ,x{y-yy)dy, (5.37) Jvi Jyi so, ~(y-yv)2 ivi-Vy) uo(Vi(x),y) =Po,x (5.38) 2 2 We have imposed a pressure drop to cause a flow in the positive u direction, i.e. UQ > 0. This means our leading order velocity is Mm, y) = ^ [(vi - yy)2 - (v - yyf] • (5-39) Now, recall r = (4,,o + 0(5))$ = 2/bo,x|(l + 0(5))l2. (5.40). 38 Chapter 5. Asymptotic Solution Clearly r is monotonically increasing with y. This means r is largest near the wall, and r > B on y € (yy,yi\- This means B uo = 7T- {(Vi ~ Vy? ~ (V ~ Vyf] , V £ (Vy,yi]- (5-41) zVy For the unyielded region, we have u constant. For continuity of u at yy, we need B uo = ^(yi-yy)2' 2/e[0,%]. (5.42) Note, according to (5.35), at leading order, there is always a plug region. However, if yy > yt, then the pressure drop is not enough to force movement and UQ = 0. Now, we use the imposed flow rate across the gap, (5.30), and the expressions for UQ, (5.41), and (5.42) to get B B ir(Vi - Vy? + / 7T-[(Vi ~ Vy? - (v - Vy?] dV = 1-z Jyv zVy Integrating, expanding terms, and multiplying both sides by yy gives !*2-(^ + D* + ^ = 0 . (5-44) (5.43). If we let then Vv = j , (5-45) %>-(?£+yi)e + ?f = 0. (5.46) Bvz Divide (5.46) by and let B = By\ to get 2^ 3 - (3 + ^ ) ^ 2 + l =0. (5.47) B Call the above polynomial P(£) and remember that B > 0. Note that, P(-oo) < 0 < P(0), P(0) > 0 > P( l ) , P( l ) < 0 < P(oo). This means that there is exactly one zero of P(£) in each of the intervals: (—oo,0), (0,1), and (1, oo). Now, 0 < yy < yi which means that £ > 1, leaving exactly one choice for £ r 0 0 t = £(B). 39 Chapter 5. Asymptotic Solution 5.4.1 Zeroth Order Resul ts The root £root(B) is easy to find numerically using e.g. bisection. In F ig. (5.1) we plot yy and yi vs. x for various values B. The same is done in Fig. (5.2), but with an increased value of h. We note that with a fixed B, for each x we will have a different value of B and hence a different yy(x). It is yy(x) that characterizes UQ, through (5.41), (5.42) and defines po,x through (5.35). It is easy to see from Fig. 5.1 that yy closely follows the shape of yi- It is difficult to tell from these pictures, but in fact, the value yi(x) — yy(x) is n o t constant in x. As B increases, yy shifts closer to the mud interface. Furthermore, from our plots of |po,x| in F ig. 5.1 we can see that for higher Bingham numbers, a larger pressure drop is needed to drive the flow. Also, we notice that |po,x| is larger at the narrower parts of the channel. 40 Chapter 5. Asymptotic Solution h = . 0 5 B=1 B = 5 1 B = 2 0 X X Figure 5.1: Plots of yy and |po,x| for B = 1,5 and 20 with h = .05. h = . 1 5 -0.5 0 0.5 -0.5 0 0.5 X X Figure 5.2: Plots of yy and |po,x| for B = 1,5 and 20 with h = .15. 41 Chapter 5. Asymptotic Solution 1.5F 0 0.2 0.4 0.6 0.8 1 Figure 5.3: Plots of un(0,y),uo(\,y) and uo(^,y) for various amplitudes of y%{x). B = 5 was used in all 3 plots. To examine uo, we look at plots of uo(y) for x = 0, \ and \ superimposed on the same axes in Fig. 5.3. Note that uo has different "plug speeds" for different yi(x). Hence if y'^x) ^ 0, then uo,x 7^ 0 for y < yy, and therefore UQ does not have real plug regions at all. For larger h we see that the difference in "plug speeds" is greater over the length of the gap. This leading order inconsistency arising from the lubrication scaling has been known and studied for some time now. See [1], [2], [6] and [5] for some discussions. For the remainder of this paper,we will use the notation uo(yy) = upp(x). This notation emphasizes the fact that the "plug speed" in the zeroth order solution varies in x. The "pp" superscript stands for "pseudo plug". 42 Chapter 5. Asymptotic Solution Figure 5.4: This is an example of a plot of UQ and the asymptotically corrected solution for an instance when up < u^. This figure illustrates why we must have yy < yr when up < in order to maintain a unit flow rate across the half-gap width. In our problem we know that if the amplitude of the sinusoidal wall variation, h, is small enough, there will be an intact plug for the entire length of the slot. This plug will have some uniform velocity u = up, yet to be determined. In order to capture this behavior, we will assume that there is a real yield surface, yr, inside of which u = (up,0). We will then construct 0(5) and 0(62) corrections to the outer solution that connect to the true plug. For any given x, the outer solution has a "pseudo-plug" speed of uo(yy) = upv(x). If uvp(x) < up, then to preserve a flow-rate of 1, would require yr < yy. However, to have upp > up would require yx > yy- If uw = up, then we will have u = UQ. See Fig. (5.4) to get a visual demonstration of this. Note that although up is constant, it has not yet been determined. This is dealt with in §5.8. Without knowing up, we can assume that it is in the range: < up < (5.48) x=0 since these are the plug velocities of the outer solution. 43 Chapter 5. Asymptotic Solution We will treat the cases yr < yy and yr > yy separately as they will require slightly different solutions. For the case yr yy, we will match the plug and the plastic region at a single point, yr-Later, in §6.1 we will construct a solution for the case when the wall amplitude is large enough to make the plug actually break. The criteria for judging when the plug will actually break is derived in §6.2. 5.5 Corrected Solution for yT < yy In this section we assume that we are in a part of the channel with upp < up and hence yr < yy-The velocity profile will have the form: u = Up for y G [0, yT] , (5.49) u = Up + 62u2(x,y) + ... for y G (yT,yy) , (5.50) u = uPP(x) + 5ui(x,y) + 52U2(x,y) + ... for y G [yy, yt] , (5.51) and we shall see yy — yr ~ 0(5). We will begin by constructing ui and u2. In the next section, we will construct the inner solution, u\, and match the three layers. Only after that will the case yr > yy will be treated. The 0(S) problem in the yielded region is: Ul,yy-Pl,x = 0, (5.52) Pi,y = 0, (5.53) ui,x + vity = 0, (5.54) ui(yi) = vi(yi) = 0. (5.55) We make the following definitions: ui(yy) = u\, (5.56) rVi ui = yy u ~ UQ + 5u° + S2u2, y£(yr,yy) u ~ up + 62u2, y € [0, yT] u ~ Up, where u°0 u°2 yr = y y - Sfr = y y - (yy + 2yi) _p i •o B 2yy [(yi - yyf - (y - yy)2. u2 (y - yi) l(yi-yy)2 [ (y-u. (3y - 2yy - yi) + 6yy (y - Vy) (yi - yy)' (3y-2yy-yi), (yi-yy)\' (5.93) (5.94) (5.95) (5.96) and u\ = up - u°0{yy) S B (y-yr)2 Px = 2yy B | 5 Qui f 2yy yy (yi - yy)2 \ (yi - yy) s2 + i)+52 &u*2 + o(ss). (yi - yy)2 (5.97) (5.98)-(5.99) (5.100) The only thing not yet determined is the actual plug speed, up, which will be found in §5.8. F ig . (5.5) shows a plot of UQ and u corrected to 0(52) vs. y on a cross-section taken in a part of the channel where yr < yy- To get the plug speed for this plot, the method presented in §5.8 was used. 50 Chapter 5. Asymptotic Solution 5.6 Corrected Solution for yy < yr To construct a solution when up < vP0(yy), we assume the following structure: u = Up for y £ [0, yx], (5.101) u = u0(yi,y) + Sui(x,y) + 82u2(x,y) for y £ (yx,yi\- (5.102) We must choose u% and u2 so that velocity and stress are continuous at the point yr, and the unit flow-rate is maintained. This differs from our approach in the case yy > yx where we had an inner transition layer, but the results will look very similar. In order to have continuity of velocity at yx we need Up = u0(yx) + Sui(yT) + 82u2(yT). Recall, because h <~ 0(6), we have \up - uo(yy)\ ~ 0(8). Furthermore, \yx — yy\ ~ 0(8) which implies that \uo(yy) — uo(yx)\ ~ 0(82). This suggests that we let u\(yx) compensate for the 0(8) difference between up and uo(yy), and let u2(yx) compensate for the 0(52) change in un over the interval (yy,yx), i-e-ui(yT)= 1(up-uPP(x)) =u*1: (5.103) 1 6 MVT) = -p (MVy) ~ MVT)) = u2. (5.104) In order to maintain a unit flow-rate, we need / Up + / u0 + Sui + 82u2 = 1. (5.105) 'VT However, we have already constructed UQ in such a way that 10 Jyy Subtracting (5.106) from (5.105) we see / u0+ / u0 = l. (5.106) JO J „ rvv rvr rVi ryi / [up-u0(yy)}+ / [up-uo(y)} + 8 ul + 52\ u2 = 0. (5.107) Jo Jyy JyT JyT 51 Chapter 5. Asymptotic Solution Furthermore, we can use (5.103) to write the first term in (5.107) in terms of u*. Doing this we get 5yyu\ + SUKVT - yy) + ^-(yT - yvf + T « i + 52 /""" u*2 = 0. (5.108) OyV JyT JyT Note the first two terms can be combined to make 5yru\. Balancing orders in (5.108) tells us about qi and q2. j wi = -yru*i = 9 i , (5.109) JyT u2 = 0 =92. (5.110) JVT The 0(5) problem presented at the beginning of §5.5 still applies here, except now u* = UI(J /T)-The solution is now = dj±(y-yT)(y-yi) + / ox 2 \yT-yiJ The definition of q\ has also been altered to use yr instead of yy in the lower limit of integration. Integrating the new expression for u\ relates PitX and q\. Jy Ul = ~ h { V i ' VT)3PX'X + V^-YL = 9 1 • ( 5 - U 2 ) Solving for pitX gives an expression depending on yr, u\ and q\. _ 12 yy. The parameters used were B = 5 , h = .05 , 5 = .075. The plug speed was determined by the method presented in §5.8. 5.6.1 Composite Solution for yr > yy The complete description of the solution for yr > yy is: y > yr u~ UQ + 8u° + 62U2, V e [0, yr] u = up, where U°2 yr = yy- (yy + tyi) Lp_ i pp 1 0 B zVy = Uo f (y-yj) 1 (Vi - yr)2 (y - y%) [y% - yr)2 3y - 2yT -y% + Qyr [3y - 2yT - y{], (y - yr) {yi-yr)!' (5.122) (5.123) (5.124) (5.125) u2 = Uv Un UQ > -u0(yT) 52 (5.126) (5.127) 54 Chapter 5. Asymptotic Solution p„ and p Corrected to Order 5 -0.5 -0.4 -0.3 -0.2 -0.1 Figure 5.7: A plot of po,x and the corrected pressure gradient for B = 5 and h = .05. The plug speed was determined by the method presented in §5.8. The value of S does not affect this plot. and B 6ut Vv (Vi - VT)2 V (Vi - VT) 5.7 C o r r e c t e d P r e s s u r e G r a d i e n t + l ) + S > ^ L ^ + 0(5*). (Vi - yr)2 (5.128) To get the corrected pressure gradient, (5.100) must be used for regions where yr < yy and (5.128) must be used for regions where yx > yy. Fig. (5.7) shows a plot of po,x and also px corrected to 0(52). Notice that the correction terms give px a significantly greater amount of variation in x than is seen in £>o,x-It is worth noting that the value of 5 does not directly affect px, but the amplitude of the wall variation, h, does. This is because of p i , x having a factor of u* and p_,x having a factor of u^. For illustration, assume we are in a region where yr > Vy, then if we were to use the definitions of u* and u2 in (5.128) we would get, pp (yi - yr)2 Px B Q[up-ulf] ( 2yT , ,\ , o[-u-0- -uayyT)\ , ^ 3 yy (yi-yr)2 \(yi-yr) A similar thing happens in regions where yr < yy-+ 1 | + « l f - " o y + 0 ( i » ) , (5.129) 55 Chapter 5. Asymptotic Solution 5 . 8 Plug Speed Determination 5.8.1 Derivation of Method to Determine the P l u g Speed Now, we devise a method to determine the plug velocity. Once the plug speed is determined, we will have a totally complete solution for the flow when the plug is intact. We assume that there exists an intact plug for the calculations in this section. We will see in Chapter 6 that this assumption is not true if h is too large. First, we refer to [7]. In §6, Prager develops the following minimum principle for stress of a slow visco-plastic flow. He states: / / u is the solution velocity field, then the solution stress field minimizes the following functional-Summation is implicit in the functional above with u\ = u and u2 = v. Different plug speeds give different stress fields, which means we should choose up so that it minimizes K. Using the fact that u2 = 0 throughout all of dfl lets us write the second term as follows, Many more terms fall out when the rest of the boundary conditions are applied to u and r y leaving, (5.130) - 2 / o-ijUiUj ds = - 2 / {u(-p + Txx) , uTxy) • n ds. 2 / (u(-p + TXX) , UTXY) • n ds = 2 / [up]x=i dy-2 [ u p ] 1 = . i dy. However, we have shown that py ~ 0(53), so if we keep only the 0(52) terms and larger, then we can pull the pressure out of the integral to get But we also know: 56 Chapter 5. Asymptotic Solution That means the entire boundary integral reduces to r-i/2 ri/2 px dx = 2 / 1/2 7-1/2 Note that the 0(1) term does not depend on the choice of up, so it does not need to be considered in the final functional. /•1/ rl/2 2[p(l/2) - p ( - l / 2 ) ] = 2 / pxdx = 2 po,* + upp(x). Furthermore, the flow is non-inertial, and we have taken yi(x) = 1 — hcos(2irx), so our solution is symmetric about x = 0. Therefore, there must be some xp such that: up > upp(x) for |a;| > xp, Up = upp(x) for | j = xp, Up < upp(x) for |a;| < xp. From this point forward, we will restrict our attention to positive values of x only. We can do this because our solution is even about x = 0. From our asymptotic solution, we have /•1/2 rVi / / [\T — B \ + T — B] dy dx = (5.132) Jo Jo B\ + T° - B]2 dy dx I Jo JyT 'VT ,1/2 + / [\T%-B\+T%-B]2dy + / [\T° — B\+T° — B\ dy Jxp Jyr Jyy dx Let us first assume that we are in a region where x > xp which means yr < Vy Recall that in the inner region TXY = - B + 5u2£ + 0(52), and Therefore r ' = \riy\ + 0{52) =B- Su24 + 0{52). rVy fVy / [\Ti-B\+Ti-B]2dy= / AS2u2u + 0(S3) dy. Jyr Jyr 57 Chapter 5. Asymptotic Solution Remember u\ ^ = and £ = V^ML _ u s j n g t u j s information, we get an expression in terms of yy and yT: fVVV -B\W- Bf dy,44 [VV(y - yTf dy = 4 B 2 ^ l ~ 0(63). jyT yy JyT 6yy For the integral in the fully yielded outer region, we know TXX ~ 0(5) because h ~ 0(5). This means T = \rxy\ + 0(54) = \TXYFI\ + 5rXyti sign^^o) + 52TXY,2 sign(r x y > 0) + 0(53). Note that sign(rXj/]o) = — 1- Therefore - fV\\r-B\+r-B}2dy ^ JVy fVi o 2 = / [\rxy,o\ - 5rXytl - 52Txy>2 ~ B + 0(53)} dy Jyy = / ' {(\Txyfi\ ~ B)2 - W r ^ i d T ^ o | - B) + 52 [r2ytl - 2rxy,2(\rxyfi\ - B)} + 0(53)} dy. (5.133) Note that the 0(1) term in (5.133) above does not depend on up. For this reason, it can be neglected in the final functional we are considering. Dropping the 0(1) term and writing the shear stresses in terms of the pressure gradients, the RHS of (5.133) becomes / -25p1:X\p0iX\y(y - yy) + 52 \p\xy2 - 2p2,x\p0,x\y(y - yy)] + 0(53) dy. Jyy (5.134) Now, if we were in a region where «r > yy, the integral would have a lower limit of yx instead of yy. However, we shall see that this only makes a 0(53) difference to (5.134). You can see this by noting that, / -2<5pi, I|p0,x|y(y - yy) + 52 [p2hxy2 - 2p2,x\po,x\y(y - yv)] + 0(53) dy JyT = \ j -i8p\,x\po,x\y(y - yy) +52 [pi,xy2 - ^P2,x\po,x\y(y - yy)} + o(63) dy [Jyy rvr 1 - J -26pi,x\po,x\y(y - yy) + 52 [p\xy2 - 2p2,x\po,x\y{y ~ Vy)] + °(^) dv t; 'Vv 58 < Chapter 5. Asymptotic Solution and, the integral from yy to yr above is 0(83) because the interval of integration is 0(5) and the term (y — yy) is 0(5) on this interval. Therefore, the expression in (5.134) is an accurate approximation for all x. However, the definitions of p i x and p2,x vary depending on whether x < xp or x > Xp. Now, by doing the integration and expanding the result it can be shown that, rVi rVi \Po,x\ / y(y - yy) dy = u0dy = 1. (5.135) Jyy Jo This means that (5.134) can be written as -26PUx + 52 [plxy2 - 2p2,x] + 0{53) dy. (5.136) When we include the boundary integral given by (5.131), the 0(5) terms cancel and the resulting functional dependent on up to be minimized is £ ( « p ) = [1/2 r PlxV2 + 0(53) dy dx. Jo Jyy The integration in y can be done exactly. The higher order terms are then dropped, giving the final form of the functional, fC{up) ~ Jo ^ (Vi - Vy) dx. (5.137) 5.8.2 Results of Plug Speed Determination Method A plot of fC(xp) is shown in Fig. 5.8. There is a clear minimum at x « .2165. Fig. 5.9 shows the relationship between h and xp. One plot has B = 5 and the other has B = 20. Two things are evident from these plots. The first is that a higher Bingham number results in xp being closer to the symmetry line x = 0. The other is that xp < .25 for h > 0. This has been validated for very small values of both B and h that are not shown on in Fig. 5.9. For any Bingham number, as h —• 0, we have xp —> .25. Another interesting feature of the functional given in (5.8) is that it is independent of 5. How-ever, we shall see in Chapter 6 that the existence of an intact plug is highly dependent on 5. 59 Chapter 5. Asymptotic Solution The Plug Speed Functional 15001 1 1 1 1 1 1 "0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 X p Figure 5.8: A plot of the plug speed functional using B = 5 , h = .02 , 5 = .05. This means that S strongly affects whether or not (5.8) is valid, but does not affect the value that (5.8) chooses for up. This comes from the fact that po,x,Pi,x and p2,x only depend on yi, yy, and up. The values obtained form the numerical simulations of Chapter 4 are shown by diamond shaped markers in Fig. 5.9. We see that the values from the numerical results are higher than the values from the analytic method derived in this section. We see that in the numerical results, we still have xp < .25 and the trend that xp decreases as h increases. We also see that for all h shown, xp is smaller for B = 20 than it is when B = 5. The relative difference between the analytic and numerical method is about 10 percent for B = 5 and h = .05. The relative difference seems to be less when either B or h is decreased. 60 Chapter 5. Asymptotic Solution |x | vs. Amplitude of y. Variation 0.251 1 E 1 1 ' 1— 0.16r I i i i i i I 0 0.01 0.02 0.03 0.04 0.05 0.06 h Figure 5.9: Plots of plug speeds for various amplitudes for B = 5 and B = 20. These plots are not affected by the value of 5. Some values determined from the numerical simulations of Chapter 4 are marked with diamonds of the appropriate color. 61 Chapter 6 Plug Breaking Criteria 6.1 A Solution for the Broken Plug The term Pseudo-Plug Region characterizes a region in our domain where the stress on the centerline is slightly above the yield stress of the fluid. The result is a flow structure that resembles that of the plug region, but is actually yielded. We will not present the complete solution here because an inner layer is needed to match this solution to the fully yielded region, and a complete solution is not necessary to develop our plug breaking criteria. Al l that is needed is an expansion for the normal stress under the pseudo-plug flow structure. To capture the behavior of the pseudo-plug region, we take expansions of the following form: p = po + Spi + S2p2 + ... , u = u0(yi(x)) + 5u1(x,y) + 62u2{x,y) +... , (6.1) v = vo + 5vi + 82V2 + ... . This means that uv ~ 0(5) in the pseudo-plug. The result of this expansion is a solution that is yielded at 0(5) in both x and y (recall y'iix) ~ h ~ 0(5)). Craster [2] used an expansion of this type to develop a solution for the flow of a visco-plastic fluid down an inclined plane. Using the new expansions for velocity, we obtain new expansions for the strain components: ixx = 2u0jX + 52ultX + O(62), jxy = 8uity + 82(u2,y + u o , x ) + 0(53). 62 Chapter 6. Plug Breaking Criteria The second invariant of strain looks like 7 = 6 (uly + 4ulx + 5[2ulty(u2,y + uo,*) + 8u0,x«i,x] + 0{52)Y ~ B the constitutive relations to leading order become - " B 2 U O * - ~ T + ... ~0{h, (6-3) Bu ixy K j / + 4 u 0 ,x ) 5 The leading order equations are still 0(1). (6.4) P0,y =0=>p0= Po(x), d_ dy But now the expansion for rxy has changed, which when substituted gives P0,x = Q^{Txyfi)-P 0 X = J L — ^ _ _ y ( 6 .5) We know po is independent of y, so integrating both sides of (6.5) gives Bul,y VPo,x = KJ/ + 4u0,X)5 Finally, solving for u\,y we find i • 2ypoiXuotx ,„ R, Ul'y = ±7r72 2 ~ ^ v i - ( 6 ' 6 ) B Equation (6.6) can be simplified if written in terms of yy = : P0,x « i , _ = T - ? p U . (6.7)' y i 1 ' fe Because urj is also independent of y, we can integrate (6.6) in from yy and have an explicit expression for u\ depending on u*: I y2 ui(x,y) = ±2uo,_2M/l 2 + u i for V < y y (6-8) V ^y This is where we will leave the analysis. In order to complete the solution, we would have a matching layer at yy that connects the pseudo-plug solution to the fully yielded outer solution. 63 Chapter 6. Plug Breaking Criteria Plug Regions y = Yi Figure 6.1: A sketch of the domain showing a broken plug. The integration uses a closed curve around the central plug region and is broken into four separate line integrals over r i , r2, rs and n. 6.2 Derivation of Plug Breaking Criteria The general idea used to develop the criteria for when the plug breaks comes from Szabo and Hassager in their study of visco-plastic flows in eccentric annular geometries [8]. The idea is that when the plug is broken we could integrate the x-momentum equation around the central plug region and then use Green's Theorem to evaluate this integral. All of the terms can be directly evaluated except for the r2 component of the line integral (see Fig. 6.1). The integrand for the r2 component of the line integral involves TXX. When the plug is broken, we know that TXX changes orders compared to when the plug is not broken. Therefore, we can use the size of the other components of the line integral, n , r 3 and r±, which balance with r2 and which can be directly evaluated to determine the order of r2 and hence whether or not the plug is broken. This is explained in more detail below. First we integrate the x-momentum equation around the plug region (following the path r = r i U r 2 U r 3 U r 4 ' s e e Fig- 6-1). We then use the fact that rxy = 0 on r\. (6.9) 64 Chapter 6. Plug Breaking Criteria Now, define rVT&b) io Prom (6.9) we have rVT(xb) I(xb)==-5 Txx(xb,y) dy. (6.10) Jo !rO ryrixb) / -pdy+ -pdy JyT(0) JO + j -Txy(x,yT)+y'T[-p(x,yT) + S2Txx(x,yT)} dx\ . (6.11) Jxb > We will now use the force balance idea briefly discussed earlier to develop both necessary and sufficient conditions for the breaking of the plug. The RHS of (6.11) will be evaluated up through 0(6). I(xb) can be evaluated once we assume a form for Txx(xb,y). To develop the necessary condition for plug breaking, we assume that Txx(xb,y) takes on the form given by the pseudo-plug solution, i.e. TXX is given by (6.3). Then, if there exists an x0 € [0,1/2] satisfying (6.11), we say that the necessary condition is satisfied. To develop the sufficient condition for plug breaking, we assume that Txx(xb,y) = — along all of r 2 . This will mean that r(xb,y) > B along all of r 2 , and hence the plug must be broken. Then, if there exists an xb £ [0,1/2] satisfying (6.11), we say that the sufficient condition is satisfied. First we will simplify the RHS of (6.11). Ignoring terms of 0(52) and higher, we see rxb 8I(xb) = \po + <5PI]X=O2/T(0) - [po + Spi}x=XbyT{xb) - Bxb + / p0(a:)^(a:) dx. (6.12) Jo Remember that the solution pressure field is unique only up to an arbitrary constant, p. We have shown that at the maxima and minima of yi(x), that py = 0. So, we can choose p = 0 on the line x = 0 to set the arbitrary constant. We would then have po = pi = 0 on the line x = 0. This means the first term on the RHS of (6.12) drops out. Integrating the last term by parts gives fXb SI(xb) = -[p 0 + <5pi]i=X62/T(2;b) - Bxb + [pour]^ - / VO,XVT dx. (6.13) Jo 65 Chapter 6. Plug Breaking Criteria B Now we use the fact that Pox — to get Vv SI(xb) = -6Pl(xb)yT{xb) + B [ b — - I dx. (6.14) Jo yy • 6.2.1 The Necessary Condition for Plug Breaking To develop the necessary condition for plug breaking we use (6.3) to evaluate the RHS of (6.10). This gives rvr&b) I r ? / l 2 \ 2 rVy(xb) I Using the substitution sin(0) = the above equation becomes I(xb) « Byy(xb) [' cos2(0)d0 = / * 1 - cos(20)d0 = (6.16). Jo 2 J 0 4 With this methodology, we would satisfy the necessary condition for a broken plug when there exists an xb G (0,1/2) satisfying F(x) = 0, where F(x) = -5Pl(x)yT(x) + B f* — - 1 dx - 5^^-. (6.17) Jo Vy 4 When there is no xb that can satisfy this equation, the plug must remain intact. 6.2.2 Sufficient Condition for a Broken Plug To develop the sufficient condition for a broken plug, we assume Txx(xb,y) = BS. With this assumption (6.10) becomes fVr(xb) 10 1 2\ 2 dy. (6.15) r T&b) I(xb)= / Bdy = ByT{xb). (6.18) Jo This means that we would satisfy the sufficient condition for plug breaking when there exists an xb G (0,1/2) satisfying G(x) = 0, where G{x) = -Sjn (x)yT{x) + B — - 1 dx - 6ByT. (6.19) Jo Vy When this condition is satisfied, we know that the plug should be broken. 66 Chapter 6. Plug Breaking Criteria Figure 6.2: Plots of the sufficient breaking criteria with B = 5 , he [.02 , .025] and 5 = .1. The sufficient criteria is not satisfied for h — .02 or .0225, but it is satisfied for h = .025. 6.3 Discussion on Results of Breaking Criteria Both F(x) and G(x) appear to always be concave down on [0,1/2] for any choice of parameters. We always have F(0) < 0 and G(0) < 0. As h increases, the apex of F(x) and G(x) increase. For any set of parameters there exists h € (0,1) that will give an apex of both F(x) and G(x) greater than zero. Figs. 6.2 and 6.3 each show G(x) plotted for h = .02, .0225 and .025. In Fig. 6.2 we have a Bingham number of 5. We can see that the plug is not broken when h = .0225, but is broken when h = .025. Notice how in Fig. 6.3 we have increased the Bingham number to 20 and as a result the plug is not broken for h = .025. In fact, the apexes of all the curves have shifted down, away from the breaking line. The trend of the Bingham number making the curves shift downward does not continue for large Bingham numbers. It seems that when B becomes large, it will cause px to break order which results in a poor approximation for up. The poor approximation for up then causes a poor approximation for yx and the breaking criteria becomes inaccurate. 67 Chapter 6. Plug Breaking Criteria 0 -0.2 -0.4 -0.6 O -0-8 -1 -1.2 -1.4 -1.6 0.3 0.4 0.5 0 0.1 0.2 X Figure 6.3: Plots of the sufficient breaking criteria with B = 20 , h G [.02 , .025] and 5 = .1. The sufficient criteria is not satisfied for any h G [.02, .025]. The most dominant parameter in both the necessary and sufficient plug breaking conditions is S. For example, the plots in Figs. 6.2 and 6.3 have 8 = .1. For 5 = .1 , B = 20 the sufficient criteria for the plug breaking is satisfied at h « .0235. When 8 = .25 it requires h « .091 to satisfy the sufficient condition. The necessary condition with 8 = .1 , B = 20 is satisfied for h « .019, but when we increase 8 to .25 the necessary condition is now satisfied a t / i « .065. When we compare the results of this section with our numerical simulations from Chapter 4, we find that the sufficient condition for plug breaking underestimates the amplitude required to break the plug in the numerical simulation. For example, when using B = 5 and L = 10, the numerical simulations tell that the plug breaks for h somewhere between .04 and .05 (numerical results not shown). However, we can see from Fig. 6.2 that the sufficient condition for plug breaking is satisfied for some h between .0225 and .025. This trend of the sufficient condition of plug breaking being satisfied for a lower amplitude than the numerical results show a broken plug seems to hold for all parameter choices tried. However, since the necessary condition for plug breaking yields an even smaller value than 68 Chapter 6. Plug Breaking Criteria the sufficient condition, we know that the necessary condition never overestimates the proper amplitude, so in this respect, we can trust it. There is some agreement between the numerical simulations and the analytical results. Both approaches agree on the fact that a higher Bingham number requires a larger amplitude to break the plug (as long as B does not become very large). They also agree that a smaller 6 requires a smaller amplitude to break the plug. One observation that should be made is that when the amplitude is large enough to satisfy the sufficient condition for plug breaking, we see a change in appearance of the yield surface in the numerical simulations. For example, 5 = .25, B = 5 and h = .085 satisfies the sufficient condition for plug breaking. These are the parameters used in Fig. 4.2. When we look at this figure, we see that the yield surface near x = —1/4 (for x scaled) is curved in towards the centerline. However, delta = .25, B = 5, h = .05 are the parameters used in Fig. 4.5. These parameter choices does not satisfy the sufficient condition and we can see from the figure that there is no dip in the yield surface near x = —1/4 (for x scaled). This suggests that there is probably some transitional change in the flow structure between the "intact plug" solution and the "pseudo plug" solution (perhaps characterized by the size of \yy — VT\)- This could account for the underestimation of the sufficient condition for plug breaking. 69 Chapter 7 Summary Both numerical and analytical techniques are used to solve the unmodified Bingham model for slow flow in a long-thin gap of slowly varying width. We have developed an analytical method to determine when an unbroken plug exists throughout the length of the gap as well as a complete description of the flow when the plug is intact. This flow description uses an asymptotic scheme based on the lubrication type scaling of the dimensional equations. This solution does not exhibit the classical inconsistencies that are usually present for a visco-plastic fluid scaled in this manner. The inconsistencies are resolved by assuming that there is an intact plug and constructing first and second order corrections to the flow that result in a fully rigid region. In order to attain the speed of this rigid solid a minimization principle for slow visco-plastic flow is used. A method for numerical simulation based on finite elements is also presented, as well as results. This augmented Lagrangian method is ideal for a Bingham fluid because it permits the use of the full Bingham model and can compute the fully rigid regions. 70 Bibliography [1] N.J. Balmforth and R.V. Craster. A consistent thin-layer theory for Bingham plastics. J. Non-Newtonian Fluid Mech., 84, 1999. [2] N.J. Balmforth, R.V. Craster and R. Sassi. Shallow viscoplastic flow on an inclined plane. J. Fluid Mech., 470, 2002. [3] G.K. Batchelor. An Introduction to Fluid Dynamics. Cambridge University Press, 1967. [4] R. Glowinski. Numerical Methods for Nonlinear Variational Problems. Springer-Verlag, New York, 1984. [5] C C . Mei and M . Yuhi. Slow flow of a Bingham fluid in a shallow channel of finite width. J. Fluid Mech., 431, 2001. [6] J .M. Piau. Flow of a yield stress fluid in a long domain: Application to flow on an inclined plane. J. Rheol., 40(4), 1996. [7] W. Prager. On Slow Visco-Plastic Flow. Office of Naval Research Contract N7onr-35801. Taken from: Studies in mathematics and mechanics, presented to Richard von Mises by friends, colleagues, and pupils., 1954. [8] P. Szabo and O. Hassager. Flow of viscoplastic fluids in eccentric annular geometries. J. Non-Newtonian Fluid Mech., 45(%), 1992. [9] S.D.R. Wilson and A.J . Taylor. The channel entry problem for a yield stress fluid. J. Non-Newtonian Fluid Mech., 65, 1996. [10] Th. Zisis and E. Mitsoulis. Viscoplastic flow around a cylinder kept between parallel plates. Non-Newtonian Fluid Mech., 105, 2002. 71