S L O W F L O W O F A B I N G H A M F L U I D IN A G A P O F S L O W L Y VARYING WIDTH by DANIEL PADRAIGH RYAN B.Sc. (Mathematics) University of Minnesota, 2000 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE D E G R E E OF M A S T E R OF SCIENCE in T H E FACULTY OF GRADUATE STUDIES Institute of Applied Mathematics, Department of Mathematics We accept this thesis as conforming to^the^required standard T H E 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 Chapter 1. vii Introduction 1 Chapter 2. Equations and Scalings 2.1 Dimensional Equations 2.1.1 The Equations of Motion 2.1.2 Domain and Constraints 2.2 Non-Dimensionalization 2.3 The Lubrication Approximation 3 3 3 4 5 7 Chapter 3. Numerical Solution 3.1 Weak Variational Formulation 3.1.1 Derivation of the Functional to be Minimized 3.1.2 Determining A P 3.1.3 The Minimization Problem 3.2 Solution of Augmented Lagrangian 3.2.1 The velocity subproblem 3.2.2 The q sub-problem 9 12 12 15 15 17 18 20 Chapter 4. Numerical Results 4.1 Notes on Numerical Implementation 4.2 Error in a Straight Channel 4.3 Results for Varying Gap Width 21 21 22 24 Chapter 5. Asymptotic Solution 5.1 Geometric Assumptions and Simplifications 5.2 Lubrication Rescalings 5.3 Regular Expansion for the Yielded Region 5.4 The Zeroth Order Solution 5.4.1 Zeroth Order Results 5.5 Corrected Solution for yr < y 5.5.1 Solution in the Transition Layer for yr < y 5.5.2 Matching for yr < y 5.5.3 Composite Solution: yr < Vy 33 33 34 35 37 40 44 45 47 50 y y y iii 5.6 Corrected Solution for y < yr 5.6.1 Composite Solution for yr > y • 5.7 Corrected Pressure Gradient 5.8 Plug Speed Determination 5.8.1 Derivation of Method to Determine the Plug Speed 5.8.2 Results of Plug Speed Determination Method y y 51 54 55 56 56 59 Chapter 6. Plug Breaking Criteria 6.1 A Solution for the Broken Plug 6.2 Derivation of Plug Breaking Criteria 6.2.1 The Necessary Condition for Plug Breaking 6.2.2 Sufficient Condition for a Broken Plug 6.3 Discussion on Results of Breaking Criteria 62 62 64 66 66 67 Chapter 7. 70 Summary Bibliography 71 iv List of Tables A table showing the relationship between grid resolution (k elements in each direction) and the L 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 4.1 2 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 A plot of the logarithm of the element size = h vs. the logarithm of the L 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 4.1 2 5.5 5.6 5.7 5.8 5.9 Plots of y and \po,x\ for B = 1,5 and 20 with h = .05 41. Plots of y and \po, \ for B = 1,5 and 20 with h = .15 41 Plots of uo(0,y),uo(j,y) and uo(^,y) for various amplitudes of yi(x) 42 A plot illustrating the relationship between the plug speed and the position of yx relative to y 43 A plot of UQ and u corrected to 0(5 ) on a cross-section where yr < y 49 A plot of uo and u corrected to 0(5 ) on a cross-section where yr > y 54 A plot of po, d the corrected pressure gradient for B = 5 and h = .05 55 A plot of the plug speed functional using B = 5 , h = .02 , 5 = .05 60 Plots of plug speeds for various amplitudes for B = 5 and B = 20 61 6.1 6.2 6.3 A sketch of the domain showing a broken plug Plots of the sufficient breaking criteria with B = 5 Plots of the sufficient breaking criteria with B = 20 5.1 5.2 5.3 5.4 y y x y 2 y 2 y a n x vi 64 67 68 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 Chapter 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 y=0 r-r<rTJJ~TTy-T i y : Figure 1.1: A sketch of the 2-D flow that will be examined. In this problem, we know that when the gap width variation is kept small enough there must be a rigid plug along the length of the gap. Herein, we will use a lubrication scaling, and then resolve the plug inconsistency without straying from the Bingham model. This will be done by constructing first and second order corrections to the pressure and velocity fields, resulting in an asymptotic solution that describes a rigid plug with uniform speed. Note that throughout this analysis, the words "gap" and "channel" are used interchangeably. An analytical technique for determining the uniform speed of the plug will also be developed at the end of Chapter 5. It will use an extremum principle for the solution stress field to determine the appropriate plug velocity. A criteria for determining when the gap width variation is large enough to break the plug will also be presented in Chapter 6. A finite element formulation of the problem will be derived and implemented in Chapters 3 and 4. It will use the method of augmented Lagrangian, which allows accurate computation of the unyielded plug region, wherein, the rate-of-strain will be computed to be zero. The agreement between the results of the numeric and analytic approaches will be discussed in Chapters 5 and 6. 2 Chapter 2 Equations and Scalings 2.1 Dimensional Equations We start with the following equations for the flow of an incompressible Bingham fluid through a 2-D gap. The flow is modelled by the Navier-Stokes Equations coupled with the Bingham fluid constitutive laws. 2.1.1 The Equations of Motion These equations relate the velocity u = (u , v), pressure p and the deviatoric stress-tensor f y . The fluid is assumed to have constant density p, and to be incompressible, and is governed by du ^du „du dp d „ dv „dv „dv. dp d „ du dx + d „ d „ dv dy (2.1) (2.2) (2.3) The rate-of-strain tensor , 7 y , and its second invariant, 7 , are defined as du dy y ixx - dv dx 2, g£ (2-4) (2.5) (2.6) 3 Chapter 2. Equations and Scalings Let f y be the yield stress of the fluid, jl be the plastic viscosity, f y be the deviatoric stress tensor, and f be the second invariant of f y , i.e. the stress. The Bingham fluid constitutive laws are Tij = (A + -^)7ij <=> 7 7 = 0 <S=> + > Ty, (2.8) f < fy, 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, f . y The areas where f < T have a zero rate-of-strain, hence they move like a rigid solid, and will V 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 = U D V f , 0 (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), u = v = 0 on y = yi(x). 4 the channel walls. This means (2-H) 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. = 0 on y = 0. r xy 2.2 (2.12) 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 V = UQV u = UQU i b p = ~ U 0 Uo. Ixy = 15 lyx UQ. 7xx — lyy = 7 p (IUQ Ty X lxy = jlUo -f7 T yx b lyx - b lxx Uo. b yy - T m UQ. ^b 0 b r yx fiUo Uo. l D p,U T = 1 5 ^ fiUo b Txx r vy Jj-Uo b T (2.13) Chapter 2. Equations and Scalings Substituting these non-dirnensional variables into the momentum equations yield du di du du. dp + v—) = dy dx dx dv di dv dx d dx dv. dp d dy dy dx~ du dv = 0. | dy dx Txy + d Txy dy (2.14) d Tyy dy (2.15) (2.16). i 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, du dy ixx dv dx = 2^, 2— dx' 2— dy' (2.19) (2.20) By introducing Bingham number,B, our constitutive relations become T« = (1 + ? ) 7 « r>B 7 = 0 <=^ T (2.21) < 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 = ^ . The condition of incompressibility gives j (2.22) = -j , xx yy which in turn gives T = -r . xx vy 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 = (r + r )l 2 2 xy 6 xx (2-24) 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 halfwidth, 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 = DUo —-—V L 7xy lyx 7xx xx T — - xx T 7KK T = ——T D Substituting these rescaled variables into the momentum equations yields T~xx ~l~ d_ n T y, X (2.25) (2.26) dx dy (2.27) The dimensionless small parameter 5 is defined by (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=(i 2 xy + S i )K 2 2 xx (2-29) and T = (T 2 XV + 5*T J. 2 (2.30) This scaling will let us use a perturbation expansion with 8 as the small parameter to analytically study the flow. 8 Chapter 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, Fig. 3.1. 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 (3.1) (3.2) 0 = (3-3) Vu, (3.4) dui Tij — dxj 9 duj dxi' (3.5). Chapter 3. Numerical Solution Figure 3.1: A sketch of the half-gap under consideration. F i , r 2 , r 3 and T4 are the contours making up the boundary of the domain Q. The boundary conditions imposed on u and T y ( u ) are the following: = o, T~xx Ul = 0, U ~-= U2 = 0, U2 = 0, = 0, on Ti = 0, on T = 0, on Ty -0, on 2 xx X X = L/2, y G (0, (L/2)), r = x e [-L/2.L/2], r = X = -L/2, ye r = X G {-L/2,L/2], 2 3 4 We note that the periodicity conditions on u W y = (x), yi (0, (-L/2)), yi y = 0. (3.6) (3.7) (3.8) (3.9) 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 f re u lect to this problem. However, since the geometry is exactly the 10 Chapter 3. Numerical Solution same as before, we know that reflect _ u u Equally, we can see that the problem for u is exactly that for u rejlect reverse but with the geometry and boundary conditions reversed. Thus, we have: u {x,y) = u[ (x,y) eflect x = -u r (-x,y) = r erse (-x,y), Ul and u (x,y) = u (x,y) = u^ (-x,y) r eflect 2 2 erse = -u (-x,y). 2 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 u = 0 does not necessarily imply that x T XX = 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 T = 0. XX 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^ + — dxdy xv dx Ul 1 Extending this to the Bingham fluid by assuming that A r dx ~ xy - o 1 ± - 0 dy ~ T yy 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 D e r i v a t i o n of the F u n c t i o n a l 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. I V2 2 \V\\HHSI) dxj Here U denotes the flow domain and V is defined by: V = {v e [ ^ ( f i ) ] : V • v = 0 e fi; w = 0 on r 2 2 b f , T ; v = 0 on r , }. 3 4 2 (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 <Sy = 0. The momentum equations, (3.1) and (3.2), written in terms of the stress, <jy, take on the simplified form (Recall that it is assumed that anytime we see a repeated subscript, summation over that subscript is implicit.) Typically the pressure is only defined uniquely up to a constant, say poIt follows from (3.10) that the pressure experiences a constant pressure drop along the channel, say A P . Thus, we write p as: AP p = Po--j-x + p, where we choose po such that the modified pressure field p satisfies: p(±L/2,y)=0. 12 (3.14) (3.15) Chapter 3. Numerical Solution Multiply (3.13) by — Vi and integrate over Q to get / (iH - Vi)-^-[-pSij JQ i + nj} + ^-{ui ax - vi) dx. = 0. (3.16). L Using the product rule for derivatives we get For now, focus on the first term in (3.17). Applying Green's Theorem results in a boundary integral on dfl: ~ In ~&x~j^ ~ ^~ Vi = Ui p5ij + f {(-pSi2 + Ti2)(vi Jan - IH) , -{-pSii + Tn)(vi - ut)) • dr = When the boundary integral, IQQ, is broken up into its 4 components (one from each edge) and the essential boundary conditions for u (and v) are applied we are left with Ian = {- f I Jo + rL/2 / ( - P + T H ) ( U I - ui)\ = x T (v! - xy 1/2 1 +ni)(u)(vi - dy - [ Jyi ui)\ _ dy x= wi)li,=o dx > . 1/2 (3.18) J J-L/2 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 AP f d Tij—-[vi Jn uxj / — -Ui]dx= h / (vi - wi) dx, Jn (3.19) and using the symmetry of r y we rewrite (3.19) as f 1 / Jn o7y( 2 v AP f - ) u * = ~r / ( i - i ) u r d L Jn y u d K - (3-20) 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 /Jry(u) 7 i j (v-u) = ^ ( l = / J ^ « ( v - u , + ^ ) I ^ ( „ ) ^ ( + - u ) v / B 57«(u)7«(v) « ^ - , ( u ) . n The Cauchy-Schwarz inequality tells us ^7ij (u)7y(v) < 7(u) (v). 7 We now have for r(u) > £ : / ^ j ( u ) 7 « ( v -u)< f 57y(u)7tf(v -u) (3.21) + B f (v) - 7(11). 7 For r < 5 we do not have an explicit expression for r like we did for T > B. However, we v shall see that inequality remains true. If r < B, we know 7 (3.21) 1 y (u) = 0 and 7(11) = 0, so 1 2 y ( )7ij(v - u) = 2 y ' ( ) ^ ' ( ) r u T u v Again, the Cauchy-Schwartz inequality can be used to get the result ^ T i j ( u ) „ - ( v ) < T(U) (V). 7 7 Therefore, for the case where T < B, we have / ^ « ( u ) 7 « ( v - u) Jn ' < / r(u) (v) Taking into consideration 7(11) = 0 for T 7 Jn z < J3 we see < B f 7 (v). Jn (3.22). 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 ( v ) - 7(11) >^f7 The solution u is the function u f v- . x Ul (3.23) £ 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 ) = / ^7 (v) + B f ( v ) -^fJn Jn 2 7 1 14 L fvi. Jn (3-24) 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 for UAP a n 2 AP2 > A P i . It is clear that UAP x is a test function d vice versa. Substituting these functions for v in (3.23) for the two solutions, and summing gives: 0< / Jn < 7 ( U A P I - UAP )dx 2 2 ( A p 2 L APi) t _ Jn U A P i i l dx. Thus, we see that the integral in (3.25), effectively the flow rate, increases monotonically with the pressure drop. Indeed, provided that AP2 > A P j and UAP X 7^ UAP 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 7x (v), X a n d the second component acting like ^ ( v ) . Consider the following two constraints: dvi dv = = 9 2 1 2 -di---dy- 2 q i = W v 1 l x x { v ) 'r - M (3 26) - ) ( 3 - 2 7 ) We then replace -y(v) and 7 (v) with |q| and | q | respectively. Then the minimization is taken 2 2 to be over v e V and q £ L (Q), but with (3.26) and (3.27) satisfied. Then the problem is 2 mm V q J ( v , q ) = / i | q | + 5 / |q| - ^ 2 JQ 2 JQ LJ / «!• (3-28) 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 vqA,p (l c max £(v,q,A,p) = „ AP ^ j - | q | + B|q| - — dvi dv . (3.29) 2 V l ,dvi 2 +1 ,dvi dv , ,l dv 2 2 dv 2 , 2 \2 dv _ p - v + ^ ( V - v ) J dx . 2 V where A, q € H = L (Q) x L (fl) and p G L (Q). k\ and k are constants. And now, 2 2 2 2 vG/C = {we {R\n)} 2 : w = 0 o n r i , T , r ;w = 0 o n r } . 2 3 4 2 The functional £(v, q, X,p) is referred to as the augmented Lagrangian functional. The variables 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\ = k = 0 this functional is referred to as a Lagrangian. 2 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 then, { q - \ X ,p ) n n n 1 GHxHx 1 L {9) given : 2 known, we define {q , A " , p n + 1 n + 1 } by: •£(u ,q - ,A ,f)<£(v,q - ,A ,p ) n n 1 B n 1 n £(u^q^A^^)<£(u ,z,A n r)ii hi) n X^ n n (3.30) Vz£H,q £H, (3.31) n p ) , l ) Vv € K,u e AC, n n = X?+ A ^ + — (3-33) P p»+! = p " - p ( V - u ) , (3.34) n 2 where pi,p 2 > Q are relaxation parameters associated with k\ and k respectively. For the 2 computations presented in Chapter 4, we use p\ = p =fci= k = 1, which seems to work well. 2 2 Glowinski [4] shows that this algorithm has the following convergence properties: u —> u strongly in 1C, n q™ —> q strongly in H, yi+l _ n x pn+l_pn y > Q g l s t r o n g S T R O N G L Y I y i n JJ f N A™ bounded in H, f 1 bounded in L {fl). 2 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 u 2 fixed, then u was solved for keeping w™ fixed. In essence, _1 2 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 subscripts. 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, 0 = dCR(u + th,v,q,X,p) dt j -ph + k (V • u)h x 2 + ki[u x x -v y qi]h x + ki[u y + v - q ]h, x •y 2 (3.35) where for ease of notation. a = (-p + k (u 2 x + v ) + h[u y x -v - qi] + Xi , h[u y y +v x q ] + A ). 2 2 After applying Green's Theorem, this becomes (3.36) 0 18 Chapter 3. Numerical Solution Thus, since h is any admissible variation, we know that, = V . a 0 — + = fciV u + k u 2 = AP kiqi + \i , A - hq ) + — + V • (k v -pAP + ku +V f + — , 2 fciVV xx 2 2 2 y 2 (3.37) xx 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£ (u,v + th,q,\,p) R 0 at = -ph + k (V • u)h - k\\u - Vy - q{\h + ki[u + v - jy y y 2 x y y x q ]h 2 x —\\hy + X h dx} 2 x [ Vh-bdx, = (3.38) Jn where b is given by - ki(u - v - 51) - A i ) . b = (ki(u + v - q ) + \ , -p+k (u +v ) y x 2 2 2 x y x y Again, we apply Green's Theorem. This time the homogeneous conditions on h leave = / /i(V-b)dx. Jn 0 (3.39) This is true for any h satisfying the boundary conditions, so we can conclude, 0 = k\V v + k v 2 2 = hS7 v 2 yy +kv 2 yy + V • (A - hq , k u -p - \\ + hqi) 2 +V • 2 2 x g, (3.40) everywhere in the interior of fl. As mentioned earlier, be solved using v" 1-1 (3.37) and (3.40) will be solved in a decoupled manner. First to evaluate f. Then (3.40) 19 will be solved using u to evaluate n (3.37) g. will Chapter 3. Numerical Solution 3.2.2 The q sub-problem and we only allow q to vary, we see our minimization If we have some fixed, u , A" and n problem becomes t A q ) = \ I |q| + B f |q| + £ / |q| - / C • q, 2 2 JQ Z JQ JQ 1 JQ (3.41) where C = (Ai + h(u - v ) , A + ki(u + v )). x y 2 y x 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 and p are piecewise constant on each element. The updates 2 (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. All 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 | d 2 2 + (B|C|-|C| )d. 2 (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 ' d q = = B ^ n F > 1+ h dC. 20 0 ' (3.43) Chapter 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 = where Nj is piecewise bilinear. To obtain the element equations, we J2VJNJ 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 / where f = (k v 2 y u (k V N N j 1 2 j — p — kiqi + k N , N )dx=- i 2 + Xi , \ 2 j xx [ i - V-fJVidx-/ dx, (4.1) k\q ). 2 The form of this equation suggests that the appropriate basis functions for p, Ai, X , q\ and q 2 2 are constant on each element. This is because we want V • f ~ V i V j , and Nj is bilinear. Now 2 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 + k Nj N ) dx = 2 X x - / f • VN dx + / t ^ 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 + k Nj yN ) 2 t 21 iiV dx = - / g • ViVj dx, (4.3) Chapter 4. Numerical Results where g = (A - hq 2 2 , k u - p-Xi 2 x + 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 L (error) 10 .00581103 20 .00146366 40 .00040736 60 .000245451 80 .000131523 100 .0000631973 120 .0000401938 2 Table 4.1: A table showing the relationship between grid resolution (k elements in each direction) and the L velocity error norm. 2 In(element size) Figure 4.1: A plot of the logarithm of the element size = h vs. the logarithm of the L 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. 2 23 Chapter 4. Numerical Results -p ) L (u -u - ) L (A? 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 2 n n 1 2 + 1 - A") L (p n L (A" 2 + 1 - A?) 2 n+1 n 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 L 2 norms of {u -u - ), n n 1 (A? -A?), + 1 (A^-A^and (p -p ) n+1 n all converge to zero as n —• oo. Because of how we update Ai, A and p, this means that 2 « + < - « # ) , 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 \T \. XX When the slope of the gap wall is increased, the size of \T \ should increase as well. So, when we decrease the size XX 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 \T \ in this region. Looking at Fig. 4.6 we see that this is what XX 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. is shown in black and the yielded region is shaded gray. The plug region 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 , -1.75 -1.50 , -1.25 -1.00 , -0.75 I -0.50 I -0.25 I 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 - 1 . 7 5 I - 1 . 5 0 I - 1 . 2 5 — r— 1 - 1 . 0 0 X - 0 . 7 5 1 - 0 . 5 0 Figure 4.6: A two color plot of the plug region with B = 5,h= shown in black and the yielded region is shaded gray. 28 I 1 - 0 . 2 5 0 . 0 0 .12, L = 4. The plug region is 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 i g . 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 Chapter 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 r xy = 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(5 ) and since we 3 only use terms up to 0(S ), they will not play a role in this analysis. The assumption of slow 2 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 5 R e „ o „ 6 R e { u ^x + V dv dx du. d-y dp - / x = + P d o V x 5 dv. dp rf = ~i + . + ) r 9 T d - d V ^ + ,_ ,. ( 5 o d V ' ' „. r + W ^ + ^ = 0. dx dy 1 } ( 5 " 2 ) (5.3) Using the same rescaling on the constitutive laws for the Bingham fluid yields the following reduced constitutive laws: du dy Ixy — 7xx = 7OT = 7 ^ dv dx' (5-4). 2 (5-5) dx' dv dy (5.6) (~2 2-2 i5. V ixy ~ " ixx = +6 (5.7) The relation between the deviatoric stress and the rate-of-strain is now: i + ?)-yy T>B, (5.8). =0 T < B, (5.9) 7 where B is the Bingham number given by B = (5.10) fiUo and r = (r + 6 rl)h. 2 2 xy 34 (5.11). 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 + 6 p2 + ... , u = UQ + dui + 6 U2 V = Vo + Sv\ + 5 V2 + ... . 2 + ... , 2 2 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: 7 2 = "0,2/ + Sui, + S (U2, Ixy = ixx = 2uo, + 62u x ul + 52uo,yU (5.12) + 0(6 ), 3 + Vo ) V tX (5.13) + S 2u , + 0{S ). 2 ltX 3 2 x 5 [2uo, {u2,y + hy iV 2 y 2 y + v ) + u\ + 4ulj + 0{5 ) 3 0tX y (5.14) Taking the square root of both sides, we find 7 = Itio ,| (1 + ^ + , [ 2 ^ K , + ^ ) + ^ + 4< ] 2 x (5.15) + Now we use a Taylor expansion in 5 to get an expansion for 7-1 7-1 U ,y 2 = \U0,y\ — I U ,y I 0 V + VQ,x UQ,y 0,V U (l- 5^- - s2 °> ( ' u UQA y U2 y l,y 2 m U + 2ul 0,y '^ + Vo x 35 + 0,X 0,y '0,2/ u ^ ~ ^y + o(5 2u x U 3 (5.16) Chapter 5. Asymptotic Solution Using our expressions for the rate-of-strain and (5.8), we now find an expansion for the shear stress, r xy = {uo, + Bsign(u ,3/) + Sui 0 y tV + S U2,y + V0,x ~ Bsigli(uo )—if- 2 + For T r 8 T <v y > 1 3 (5.17) + 0(<5 ). + 5T 3 2 x 0(c5 ) XVI2 Tyy we shall only need the leading order term, 2 B 1+ Substitute the expansions for u, v, p, and r uo,x + 0(5). (5.18) into the Momentum Equations,(5.1) and (5.2), y 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= (5.19) -po,x + uo, , yy (5.20) 0 = PQ,y 0(S) : 0 = ~Pl,x 0= Pl,y. + (5-21) Ul yy, t (5.22) d 0{5 ) : 2 0= -u uo,x 0 = -P2, 0 y + - d v u y - p ,x + Qy~ xy,2 + d ~ <° T 0 0t faTxyfi 2 + Txx x (5.23) ' (5.24) Q^ yyfiT 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 u ~ 0(5) and by integrating the incompressibility condition out from y = 0 we x get vo — 0. This means the inertial terms in (5.23) should actually be in the 0(<5 ) problem. 3 36 Chapter 5. Asymptotic Solution We also get as a result of h ~ 0(5): J~ *2/>° T ~ 0(5), Txxfi — 0, Tyyfl = 0, T y,2 = X U-2,y Therefore, from the assumption of h ~ 0(5), we get the following problem at second order: 0 = - 2,x + 0(6 ): 2 P 0= 5.4 (5.25) u ,yy, 2 (5.26) p ,y. 2 T h e Zeroth Order Solution The leading order equations give: Po,y = 0 = > p o =Po(x), Po,x = J ~ ( o , i / + sign(u , )B) u 0 Y <=> 14(3,2/ = 0 To leading order, TO = |T X2/] O|. (5.27) T>B, (5.28) T < B. The continuity equation at leading order is «o,i + v , = 0. 0 y (5.29) The leading order solution should also satisfy the imposed flow constraint, i.e. mix) / Jo 10 u dy = l V x . 0 (5.30) In the region where r > B, we have at leading order P0,X = g^Txyfi- 37 (5.31) 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, y = u x So, for T 0tV + sign(u ,j/)£. (5.32) 0 > B, uo, = Po,xy - sign(u ,jy)B. y (5.33) 0 Now, when r = B there is a yield surface, call it y . In the unyielded region, uo, = 0, so for y y continuity of derivatives, and hence continuity of shear stress, we need uo —> 0 as y —» y from tV v 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 y = sign(u ,y) y • 0 (5-34) P0,x But, we know 0 < y < y%, so, y • y y = ^ , • (5.35) where r(y )o = B, and sign(uo ) = sign(po.x)- This means (5.33) can be written as y >y uo,y = Po,x(v ~ Vy)- ( -36) 5 Recalling that u(yi) = 0, we again integrate with respect to y, this time in from yi. ry / u dy= Jvi 0iy ry / PQ, {y-y )dy, Jyi (5.37) ~(y-y ) 2 (5.38) x y so, ivi-Vy) 2 2 uo(Vi(x),y) =Po,x v 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 - y ) 2 y - (v - y f] y • (5-39) Now, recall r = (4,,o + 0(5))$ = 2/bo,x|(l + 38 0(5)) 2. l (5.40). Chapter 5. Asymptotic Solution Clearly r is monotonically increasing with y. This means r is largest near the wall, and r > B on y € (y ,yi\y This means B uo = 7T- {(Vi ~ Vy? ~ (V ~ Vyf] , Vy V £ (Vy,yi]- z (- ) 5 41 For the unyielded region, we have u constant. For continuity of u at y , we need y B uo = ^(yi-yy) ' 2/e[0, ]. 2 (5.42) % Note, according to (5.35), at leading order, there is always a plug region. However, if y > yt, y 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 ir(Vi - Vy? + / z B 7T-[(Vi ~ Vy? - (v - Vy?] V d Jy Vy z (5.43). = 1 v Integrating, expanding terms, and multiplying both sides by y gives y !*2-(^ + D* + ^=0. (5-44) If we let Vv = j , (5-45) then %>-(?£ y )e + Divide (5.46) by Bv i + ?f = 0. (5.46) z and let B = By\ to get 2 ^ - ( 3 + ^ ) ^ + l =0. 3 (5.47) 2 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 < y < yi which means that £ > 1, leaving exactly one choice for £ t = £(B). y r00 39 Chapter 5. Asymptotic Solution 5.4.1 Zeroth Order Results The root £ oot(B) is easy to find numerically using e.g. bisection. In F i g . (5.1) we plot y r y and yi vs. x for various values B. The same is done in F i g . (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 y (x). It is y (x) y y that characterizes UQ, through (5.41), (5.42) and defines po,x through (5.35). It is easy to see from F i g . 5.1 that y y closely follows the shape of yi- It is difficult to tell from these pictures, but in fact, the value yi(x) — y (x) y is n o t constant in x. A s B increases, y y shifts closer to the mud interface. Furthermore, from our plots of |po,x| in F i g . 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=.05 B=1 B=5 1 B=20 X X Figure 5.1: Plots of y and |po,x| for B = 1,5 and 20 with h = .05. y h = .15 -0.5 0 0.5 -0.5 X 0 0.5 X Figure 5.2: Plots of y and |po,x| for B = 1,5 and 20 with h = .15. y 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 o,x 7^ 0 for y < y , and therefore UQ does not have real plug regions at all. For larger h we u y 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 o(y ) = u (x). u pp y 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 u < u^. This figure illustrates why we must have y < yr when u < in order to maintain a unit flow rate across the half-gap width. p y p 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 = u , yet to be determined. In order to capture this behavior, we will assume that p there is a real yield surface, yr, inside of which u = (u ,0). p 0(6 ) 2 We will then construct 0(5) and 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(y ) = u (x). y pv If u (x) then to preserve a flow-rate of 1, would require yr < y . However, to have u y require yx > y - If u y w < u, vp pp p > u would p = u , then we will have u = UQ. See Fig. (5.4) to get a visual p demonstration of this. Note that although u is constant, it has not yet been determined. This is dealt with in §5.8. p Without knowing u , we can assume that it is in the range: p (5.48) < u < p x=0 since these are the plug velocities of the outer solution. 43 Chapter 5. Asymptotic Solution We will treat the cases yr < y and yr > y separately as they will require slightly different y y solutions. For the case yr <y , we will have an inner transition layer of thickness < 5 connecting y the plug to the fully plastic region. For the case yr > y , we will match the plug and the plastic y region at a single point, yrLater, 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 y < y T y In this section we assume that we are in a part of the channel with u < u and hence yr < y pp p y The velocity profile will have the form: u = Up for y G [0, y ] , u = Up + 6 u (x,y) + ... u = (5.49) T 2 2 for y G (y ,y ) , T (5.50) y uPP(x) + 5ui(x,y) + 5 U2(x,y) + ... for y G [y , yt] , 2 y (5.51) and we shall see y — yr ~ 0(5). We will begin by constructing ui and u . In the next section, y 2 we will construct the inner solution, u\, and match the three layers. Only after that will the case yr > y will be treated. y The 0(S) problem in the yielded region is: Ul,yy-Pl,x = 0, (5.52) = 0, (5.53) ui, + vi y = 0, (5.54) Pi, y x t ui(yi) = vi(yi) = 0. (5.55) We make the following definitions: ui(y ) = u\, (5.56) ui = <?i, (5.57) y rVi 'Vv 44 Chapter 5. Asymptotic Solution and note that for given uj and q\, we can solve for u\ and by integrating (5.52) with respect pi tX to y. We find that: (V-Vy)(v-Vi) n , .,• ( 2 y V-Vi (5.58) \ y - 2/i y PX,X + « I — « i= - T _ W - W (5.59) 91, and Pi,x = 12gi - 6M? - %) (5.60) (Vi ~ Vv) ' 3 2 Substituting (5.60) into (5.58) gives the following formula for u\ depending on u\ and q\ only: Ml (y - vi) (Vi - %) Ui(3w - 2y - yi) - 6qi y 2 (y - y ) y '(Vi-Vy)] ' (5.61) The values u* and q\ will be found in the next section by matching u\ to an inner expansion in a transition layer between the plug and and the fully yielded flow. Under the assumption that h ~ 0(S), the 0(5 ) problem looks exactly like the 0(5) problem. 2 Therefore, we have u = 2 (y - yi) {yi - yy) u (3y 2 ~1y - yi) -6q y 2 2 (yi-yy)l (5.62) Again, u and q are defined by 2 2 Mvv) rut Jyy u 2 "2, (5.63) Q2- (5.64) 5.5.1 Solution in the Transition Layer for y < y T y In the transition layer, we take the expansion v} = u + 5 u (x, p 2 2 y) + 5 u\(x, 3 y) + .. (5.65) As we will see, this will result in a layer where 7 ~ 0(5) and r xy = B + 0(5). This layer has thickness 0(5) so we rescale y like y = y (x) + 5£, T 45 (5.66) Chapter 5. Asymptotic Solution and we define Vy(x) -yr(x) (5.67) W i t h this expansion for u, we can find the size of v through integrating the continuity equation out from the true yield surface. «*= fuidy^ f 8 u\ ~0{5 ). V 'VT 2 (5.68) 3 x VT J We find that v ~ 0(5 ) because v = 0 in the plug and we have only integrated over an interval l 3 of size 5. This makes v such small magnitude that it can be neglected i n the expansions for % rate-of-strain and hence does not play a role in the transition layer. Taking our expansion for u and our knowledge that v <~ 0(<5 ), we get the following expressions l 3 l for our rate-of-strain components: (5.69) (5.70) ixx = S 2u\ + ... , 2 tX = 5 (uy 2 r = Y 1 ' + 5 2u\^ui 2 + ... , 3 4 (\ 1 -^ 6\u%J u, 2 , £ l \\ "2,£ 1 2 • 0{5 /n/*2 )] (5.71). . + (5.72) / x We now use equation (5.72) to get an expansion for the shear stress in the inner region. iy X 1+ B S\u:'2,51 l-5^ + 0(6' ) z 2,C 2,£ 2 iC 2,4 (5.73) B sign(t4 ) + 5u\ + 0(S ). A The leading order equations for the transition layer are now: d i (5.74) (5.75) 0 = 0 = -PO,X+Q^U , U u\(y ) 46 T = i4,»(«r)- (5.76) Chapter 5. Asymptotic Solution Because po does not depend on £, integrating equation (5.74) from £ = 0 out into the transition layer twice and applying the boundary conditions (5.76) gives the following expression for u\: 4=Po,x^ = - ^ < 0 . (5.77) 2 5.5.2 M a t c h i n g for y < y T y In this section, we will determine q\, q , u{, u\ and yr for the case yr < y - This will lead 2 y to a description of the solution that depends only on the remaining unknown quantity u . A p method for determining u is presented in §5.8. p To avoid confusion in notation, the outer solution in the plastic region is denoted with a superscript o, e.g. u°, and the inner solution with a superscript i, e.g. u . 1 The overall aim now is to match the velocity and stress at each order, within the inner layer. fVi At the same time, we must ensure that / u dy = 1. Jo Recall from the previous section that at y we have, in terms of the inner variables: y u\y ) = -5 ^-& 2 y = -B-6^ T (Vy) xy + 0(5 ), (5.78) + 0(6 ). (5.79) z Up 2 T Vy Recall from our outer solution that at y we have y y = u (x) + 5u* +5 u* + 0(S ), r° (y ) = -B + Sul (y ) u°{y ) xy y 2 pp y 2 + S ul (y ) + 2 y y (5.80)" 3 1 y 0(6 ). (5.81) 3 Continuity of Velocity when yx < y y It is evident that for h ^ 0 we have u p(x) ^ u . However, if h ~ 0(6), then the total variation p p in u^a:) is also 0(6). To see this we could take a perturbation expansion of the solution &oot(B) about the value £, roo i.e. for t(B), yi = 1. Therefore, we have that u - u (x) ~ 0(5), pp p provided that equation (5.48) holds. In the inner layer, u ~ u + 0(5 ), so we must match the % 2 p discrepancy between u and UQ(y ) p y with uj, so u * = \[u -u (x)]. pp 1 o p 47 (5.82) Chapter 5. Asymptotic Solution Then to match the 0(S ) terms in (5.78) and (5.80) we must have 2 ul = (5-83) 2y y Finding qi and q2 when y r < y y To find q\ and q we consider the total flow across the channel. For the outer solution, u°, we 2 have r ° u u° + 5ul + 6 u° dy = 2 0 = 2 l-yyU (x) + 5q +5 q pp 2 1 (5.84) + 0(S ). i 2 We also know rvv rvr rvv . / udy= u dy+ u + 5 u dy = y u Jo Jo Jyr 2 p p 2 y (5.85) + 0(5 ), J p since y — yx ~ 0(6). y Combining this with the fact that the total flow across the gap half-width is 1, we have 0 = y [u - PP{X)} + 8 y p (5.86) + 5 q + 0(S ). 2 qi U 3 2 We now recognize the expression in brackets as Su* and rewrite the equation as 0 = 5(y ul+q ) y 1 (5.87) +5q. 2 2 Matching orders gives qi = -y u\, (5.88) = 0. (5.89) y q 2 Matching the Shear Stress for yx < y y In order to match T ° y ( y y ) from (5.81) with find an explicit expression for u\ y T l x y (y ) from (5.79), we need to know u\ . y We can by differentiating equation (5.61), and using the fact that Qi = -y u* y y v 48 Chapter 5. Asymptotic Solution u„ and Corrected Solution when y < y T 0 0.1 0.2 0.3 0.4 0.5 y 0.6 y T 0.7 0.6 0.9 y Figure 5.5: A plot of UQ and u corrected to 0(5 ) on a cross-section of the channel half-width where yr < Vy The parameters used were B = 5 , /i = .04 , 5 = .05. The plug speed was determined by the method presented in §5.8. 2 2u\ (5.90) Substitute (5.90) into equation (5.81) and we have T°{y ) y = 2u\ -B-6 (yi - (y + 2^) + 0{5 ). 2 y yy) (5.91) 2 For (5.81) and (5.91) to agree at 0(5) we must have 2y u\ y B{vi (yy + = {y v - y) 2 (y y + 2y») y 2yi)H + 2y») Notice that £ r ~ 0(1) as required. 49 -1 (5.92) Chapter 5. Asymptotic Solution 5.5.3 Composite Solution: yr < yy The solution for yr < y is now complete and is given by: y u~ y >y y y£(yr,y ) UQ + 5u° + S u , 2 2 u ~ u + 6u, 2 y p y € [0, y ] u~ T 2 Up, where _p yr = y - Sfr = y - (y + 2 ) y B u° 0 y [(yi - y f 2y y y yi i (5.93) •o (5.94) - (y - y ) . 2 y y (y - yi) (yi-y ) l [ 2 y u° (y-u. u 2 2 (yi - yy)' (3y - 2y - yi) + 6y y y (y - Vy) (yi-y )\' (5.95) y (5.96) (3y-2y -yi), y u - u° {y ) p 0 (5.97) y S (5.98)- u\ = B 2y (y-yr) 2 s (5.99) 2 y and B Px = | f Qui 5 y (yi - y ) 2 y y 2y y + i)+5 2 \ (yi - y ) y * &u + o(s ). 2 (yi - yy) s 2 (5.100) The only thing not yet determined is the actual plug speed, u , which will be found i n §5.8. p F i g . (5.5) shows a plot of UQ and u corrected to 0(5 ) vs. y on a cross-section taken i n a part 2 of the channel where yr < y - To get the plug speed for this plot, the method presented i n §5.8 y was used. 50 Chapter 5. Asymptotic Solution 5.6 Corrected Solution for y < yr y To construct a solution when u < vP (y ), we assume the following structure: p 0 y u = Up for y £ [0, yx], u = u (yi,y) 0 (5.101) + Sui(x,y) + 8 u (x,y) 2 for y £ (yx,yi\- 2 (5.102) We must choose u% and u so that velocity and stress are continuous at the point yr, and the 2 unit flow-rate is maintained. This differs from our approach in the case y > yx where we had y an inner transition layer, but the results will look very similar. In order to have continuity of velocity at yx we need Up = u (yx) + Sui(y ) 0 Recall, because h <~ 0(6), we have \u - uo(y )\ p implies that \uo(y ) — uo(yx)\ y 0(8) ~ 0(8 ). difference between u and uo(y ), over the interval (y ,yx), 8 u (y ). 2 2 T ~ 0(8). Furthermore, \yx — y \ ~ 0(8) which y y This suggests that we let u\(yx) 2 p + T and let u (yx) y 2 compensate for the compensate for the 0(5 ) change in un 2 i-e- y 1 (up-uPP(x)) 6 ui(y )= T MVT) 1 = =u* (5.103) = u. (5.104) 1: -p (MVy) ~ MVT)) 2 In order to maintain a unit flow-rate, we need / Up + / u + Sui + 8 u 2 0 = 1. 2 (5.105) 'VT However, we have already constructed UQ in such a way that / u+ / u = l. 0 JO 10 (5.106) 0 Jy„ Jy y Subtracting (5.106) from (5.105) we see rvv / Jo [u -u (y )}+ p 0 y rvr / [u -uo(y)} Jy p y 51 +8 rVi ryi u + 5\ u = 0. Jy Jy l T 2 2 T (5.107) 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 5y u\ y + SUKVT - y f + T « i + 5 /""" u* = 0. - yy) + ^-(y 2 v T Oy Jy V Note the first two terms can be combined to make 5yru\. (5.108) 2 Jy T T Balancing orders in (5.108) tells us about qi and q . 2 wi = j -yru*i = 9i, (5.109) 0 =92. (5.110) Jy T u = 2 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-y )(y- ) = T yi / + \yT-yiJ 2 ox The definition of q\ has also been altered to use yr instead of y in the lower limit of integration. y Integrating the new expression for u\ relates i P Jy Solving for pi tX Ul ~h = ' { V i and q\. tX ' VT)3PX X + ^-Y V L = 9 1 • ( 5 - U 2 ) gives an expression depending on yr, u\ and q\. _ 12<?i ' ^ {yi - VT) Pl x 6u* (yi - 3 ^ yr) ' 2 Substituting (5.113) into (5.111) and also using (5.109) to eliminate qi, we arrive at an expression for u\ that only depends on yr and u*. Ul = Ui (y ~ yi) (yi - yr) 3y - 2y - yi + 6yr ^ ^ (Vi-yr). T 2 V VT (5.114) The same procedure can be used to solve the 0(5 ) problem with the new boundary conditions 2 u (yr) 2 = u 2 The solution is: U 2 = {y - yT) U2 V i 2 ^y-^-y^- 52 ( 5 - 1 1 5 ) Chapter 5. Asymptotic Solution Now we must find yr to have a complete description of the solution. This is done by matching the shear stress to 0(5) at yx- This requires y(yr) « uo,y(yr) + 5u^ (y ) = 0. (5.116) m, {yT) = — { V T - Vy) ~ (5.117) u y T Recall, 0{5). y Vy Therefore, in order to have a balance of terms in (5.116) we need , ^ U s B (yr-yy) ~y~y 6 ' = (5.118) If we differentiate (5.114), we will then be able to use (5.118) to get an expression for yr in terms of u*. Ul,y {Vi - yr) 2 2y-yr3y - 2y -yi + 3(y - y ) + 6y T t Vi T (5.119) Which means, -2wT hy(yr) = , _ y ( y (5.120) ( l/< + VT) • u 2 r)2 Combining (5.118) with (5.120) gives the following expression for yr' yr5 y y _ ^M^ B{yi - y ) { 2 y i + y T ) 2 T ^ ^ ^ _ B(yt { 2 m + y y ) . (5.121) Note that this means yr is given approximately by the same equation as the case yr < y y yr~y y {tyi + y ) v '0 Now that yr has been found, we know u\ from (5.114) and u from (5.115). 2 Note that the shear stress has not been matched at 0(6 ) at the plug interface. This discon2 tinuity is also present in the solution for the case yr < y , but in that case it occurred at y , y y where the transition layer meets the outer solution. As a result, 5 and h must be quite small (no more than about .08) for the solution to keep a smooth appearance. 53 Chapter 5. Asymptotic Solution u and the Corrected Solution vs. y when y < y Q y T r- 1 Corrected Soln 0 0.1 0.2 0.3 0.4 0.S y 0.6 y y 0.7 0.8 0.9 1 y, T Figure 5.6: A plot of uo and u corrected to 0(5 ) on a cross-section of the channel half-width where yr > y . The parameters used were B = 5 , h = .05 , 5 = .075. The plug speed was determined by the method presented in §5.8. 2 y 5.6.1 Composite Solution for yr > yy The complete description of the solution for yr > y is: y y > yr V u~ UQ + 8u° + u= u, e [0, yr] 6 U2, 2 p where p_ pp L yr = y - (y + tyi) y y i1 (5.122) 0 B (5.123) Vy z f (y-yj) 1 U° 2 = Uo (Vi - yr) 3y - 2y -y% + Qyr T 2 (y - y%) [y% - yr) [3y - 2y - y ], T Un v 2 (5.124) (5.125) 2 U u { (y - yr) {yi-yr)!' = UQ -u (y ) > 0 5 54 (5.126) 2 T (5.127) 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 (Vi - VT) Vv 5.7 2 + V (Vi - l ) + S VT) > ^ L ^ (Vi - yr) + 2 (5*). 0 (5.128) Corrected Pressure Gradient To get the corrected pressure gradient, (5.100) must be used for regions where yr < y and y (5.128) must be used for regions where yx > y . Fig. (5.7) shows a plot of po,x and also p y x corrected to 0(5 ). Notice that the correction terms give p a significantly greater amount of 2 x variation in x than is seen in £>o,xIt is worth noting that the value of 5 does not directly affect p , but the amplitude of the wall x variation, h, does. This is because of p i , having a factor of u* and p_,x having a factor of u^. x For illustration, assume we are in a region where yr > Vy, then if we were to use the definitions of u* and u in (5.128) we would get, 2 Px B y Q[u -ulf] p (yi-yr) 2 y ( T \(yi-yr) , o[-u-pp- -uayy )\ , ,\ 2y + 1 | A similar thing happens in regions where yr < y y 55 , « l f - " o y (yi - yr) 0 + T 2 + 0 ( ^ 3 i » ) , (5.129) 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 velocityfield,then the solution stressfieldminimizes the following functional(5.130) Summation is implicit in the functional above with u\ = u and u = v. Different plug speeds 2 give different stress fields, which means we should choose u so that it minimizes K. p Using the fact that u = 0 throughout all of dfl lets us write the second term as follows, 2 -2 o-ijUiUj / ds = - 2 / {u(-p + T) xx , uT ) xy •n ds. Many more terms fall out when the rest of the boundary conditions are applied to u and r y leaving, 2 / (u(-p + T XX ) , UT ) XY •n ds = 2 / [up] i x= dy-2 [up] .i 1 = dy. However, we have shown that p ~ 0(5 ), so if we keep only the 0(5 ) terms and larger, then 3 2 y 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 •1/2 2[p(l/2) - p ( - l / 2 ) ] = 2 / ri/2 rl/2 p dx = 2 / p dx = 2 1/2 7-1/2 7-1/2 7-1/2 x po,* + < W + <5 p ,x + 0(<* ) dx. 3 2 x 2 (5.131) Note that the 0(1) term does not depend on the choice of u , so it does not need to be considered p in the final functional. Now focus on the integral over the interior domain. We have different solution forms depending on whether u < u (x) or u > u (x). PV pp p p 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 x such that: p u > u (x) for |a;| > x , pp p p Up = u (x) for | j = x , pp p Up < u (x) for |a;| < x . pp p 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 Jo / [\T — B \ + T — B] dy dx = Jo B\ I Jo Jy 'VT ,1/2 + p B] 2 dy dx T / Jx + T° - (5.132) [\T -B\+T -B] dy % % + 2 Jyr [\T° — B\+T° / — B\ dy dx Jyy Let us first assume that we are in a region where x > x which means yr < Vy Recall that i n p the inner region T = - B + 5u £ + 0(5 ), 2 XY 2 and r ' = \ri \ + 0{5 ) =B2 y Su + 0{5 ). 2 24 Therefore rVy fVy / [\ -B\+T -B] dy= Jyr i i / AS u Jyr 2 2 T 57 2 + 0(S ) dy. 3 u Chapter 5. Asymptotic Solution Remember u\ ^ = y y and £ = V^ML _ u j g s n j t u information, we get an expression in terms of s and y : T f V -B\WVV Bf dy, 4 [ (y 4 jy y T yf - VV 4 B dy = T 2 Jy y ~ ^ l 0(6 ). 3 y 6 T y For the integral in the fully yielded outer region, we know T ~ 0(5) because h ~ 0(5). This XX means = \r \ + 0(5 ) T = 4 xy \T \ 5r y i + XYFI X sign^^o) + t sign(r 5T , 2 XY 2 xy>0 0(5 ). )+ 3 Note that sign(r j o) = — 1- Therefore X /] f \\r-B\+r-B} dy - V 2 ^ Vy J f Vi = / Jy [\r ,o\ - 5r y ' {(\Txyfi\ ~ B) xy X o - 5T 2 tl 2 ~ B + 0(5 )} 3 xy>2 dy y = / + 5 [r 2 - 2 W r ^ i d T ^ o | - 2r , (\r \ 2 ytl xy 2 - B) - B)} + 0(5 )} dy. 3 xyfi (5.133) Note that the 0(1) term in (5.133) above does not depend on u . For this reason, it can be p 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 / -25p \p \y(y 1:X - y ) + 5 \p\ y 0iX 2 y x - 2p , \p , \y(y 2 2 x - y )] + 0(5 ) dy. 0 x 3 y Jyy (5.134) Now, if we were in a region where «r > yy, the integral would have a lower limit of yx instead of y . However, we shall see that this only makes a 0(5 ) difference to (5.134). You can see 3 y this by noting that, / Jy -2<5pi, |p ,x|y(y - y ) + 5 [p y I 2 y 0 2 hx 2 - 2p , \po, \y(y 2 x x - y )] + 0(5 ) dy 3 v T = \ j [Jyy -i8p\, \po,x\y(y x - yy) + 52 [pi,xy 2 - ^P2, \po,x\y(y x - yy)} + o(6 ) dy 3 1 rvr - J'Vv -26pi, \po, \y(y x x - y ) + 5 [p\ y y 2 x - 2p , \po,x\y{y 2 2 x 58 ~ Vy)] + °(^) v t; d < Chapter 5. Asymptotic Solution and, the integral from y to yr above is 0(8 ) because the interval of integration is 0(5) and 3 y the term (y — y ) is 0(5) on this interval. Therefore, the expression in (5.134) is an accurate y approximation for all x. However, the definitions of p i and p2, vary depending on whether x x x < x or x > Xp. p Now, by doing the integration and expanding the result it can be shown that, rVi \Po,x\ / Jy rVi y(y - y ) dy = Jo y y u dy = 1. (5.135) 0 This means that (5.134) can be written as -26 + 5 [pl y - 2p 2 , x ] + 0{5 ) dy. 2 2 PUx (5.136) 3 x When we include the boundary integral given by (5.131), the 0(5) terms cancel and the resulting functional dependent on u to be minimized is p £(«p) = [ Jo r J 1/2 PlxV + 0(5 ) dy dx. 2 3 yy The integration in y can be done exactly. The higher order terms are then dropped, giving the final form of the functional, fC{u ) ~ J p 5.8.2 o ^ (Vi - Vy) dx. (5.137) 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 x . One plot has B = 5 and the other has B = 20. Two things are p evident from these plots. The first is that a higher Bingham number results in x being closer p to the symmetry line x = 0. The other is that x < .25 for h > 0. This has been validated for p 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 x —> .25. p Another interesting feature of the functional given in (5.8) is that it is independent of 5. However, we shall see in Chapter 6 that the existence of an intact plug is highly dependent on 5. 59 Chapter 5. Asymptotic Solution T h e Plug S p e e d Functional 15001 "0 1 0.05 1 0.1 1 0.15 1 0.2 1 1 0.25 0.3 X 0.35 0.4 0.45 0.5 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 u . This comes from the fact that p y, y po,x,Pi,x and p , only depend on 2 x yi, and u . p 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 x < .25 and the trend that x decreases as h increases. We also see that for all h shown, p p x is smaller for B = 20 than it is when B = 5. p 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. A l 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 + S p2 + ... u = u (y (x)) v = vo + 5vi + 8 V2 + ... , 2 0 + 5u (x,y) i + 6 u {x,y) 2 1 2 +... , (6.1) . 2 This means that u ~ 0(5) in the pseudo-plug. The result of this expansion is a solution that v 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 = 2u j = 8ui y + 8 (u2, xy 0jX + 52u + ltX 2 t y 62 O(6 ), 2 +uo,x) + 0(5 ). 3 Chapter 6. Plug Breaking Criteria The second invariant of strain looks like + 4ul + 5[2u y(u ,y 7 = 6 (uly lt x 0{5 )Y ~ 2 uo,*) + 8u ,x«i,x] + + 2 0 <D(S). (6.2) For regions where r > B the constitutive relations to leading order become - " B 2 U O * - ~ T + ... ~0{h, Bu 0(1). ixy Kj/ + (6-3) 4 u 0,x) (6.4) 5 The leading order equations are still P0,y =0=>p = Po(x), 0 d_ P0,x = But now the expansion for r xy Q^{Txyfi)- dy has changed, which when substituted gives P 0 X = J L — ^ _ _ y (6 .5) We know po is independent of y, so integrating both sides of (6.5) gives VPo,x = l,y Bu KJ/ + 0,X) 4u i • 5 Finally, solving for u\, we find y 2ypo uo x = 7r72 ± ,„ , t iX ' Ul y R 2~^vi- ( 6 ' 6 ) B Equation (6.6) can be simplified if written in terms of y = : y P0,x «i,_=T-?pU. y i ' (6.7)' 1 fe Because urj is also independent of y, we can integrate (6.6) in from y and have an explicit y expression for u\ depending on u*: I ui(x,y) y 2 = ±2uo,_2M/l V 2 + u i for V yy < (-) 6 8 ^y This is where we will leave the analysis. In order to complete the solution, we would have a matching layer at y that connects the pseudo-plug solution to the fully yielded outer solution. y 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 , r , rs and 2 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 r component of the line integral (see Fig. 6.1). The integrand 2 for the r component of the line integral involves T . When the plug is broken, we know that 2 T XX XX 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 r and which can 2 be directly evaluated to determine the order of r and hence whether or not the plug is broken. 2 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 2 U 3 U 4' r r r s e e Fig- 6-1). We then use the fact that r xy = 0 on r\. (6.9) 64 Chapter 6. Plug Breaking Criteria Now, define I(x )==-5 b rVT&b) rVT(x ) b T (x ,y) xx io Jo (6.10) dy. b Prom (6.9) we have ! / rO ryrixb) -pdy+ Jy (0) -pdy JO T + j -T (x,y )+y' [-p(x,y ) xy T T + S T (x,y )} dx\ . 2 T xx T Jx (6.11) > b 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(x ) can be evaluated once we assume b a form for T (x ,y). xx b To develop the necessary condition for plug breaking, we assume that T (x ,y) takes on the xx b form given by the pseudo-plug solution, i.e. T is given by (6.3). Then, if there exists an XX x € [0,1/2] satisfying (6.11), we say that the necessary condition is satisfied. 0 To develop the sufficient condition for plug breaking, we assume that T (x ,y) = — along all xx b of r . This will mean that r(x ,y) > B along all of r , and hence the plug must be broken. 2 b 2 Then, if there exists an x £ [0,1/2] satisfying (6.11), we say that the sufficient condition is b satisfied. First we will simplify the RHS of (6.11). Ignoring terms of 0(5 ) and higher, we see 2 rxb 8I(x ) = \po + <5PI] =O2/T(0) - [po + Spi} X b y {x ) - Bx + / Jo x=Xb T b b p (a:)^(a:) dx. 0 (6.12) 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 p = 0. So, we can choose p = 0 on y 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(x ) = - [ p + <5pi]i= 2/ (2;b) - Bx + [pour]^ - / Jo b 0 X6 b T 65 VO,XVT dx. (6.13) Chapter 6. Plug Breaking Criteria B Now we use the fact that Pox — to get Vv + B [ — - I dx. Jo y • SI(x ) = -6 (x )y {x ) b Pl b T (6.14) b b y 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) r?/l \ 2 I Using the substitution sin(0) = I(x ) « By (x ) b y b 2 rVy( b) x I 1 2\ 2 dy. (6.15) the above equation becomes [' cos (0)d0 = Jo 2 J 2 / * 1 - cos(20)d0 = (6.16). 4 0 With this methodology, we would satisfy the necessary condition for a broken plug when there exists an x G (0,1/2) satisfying F(x) = 0, where b F(x) = -5 (x)y (x) Pl T + B f* — - 1 dx - 5^^-. Jo Vy (6.17) 4 When there is no x that can satisfy this equation, the plug must remain intact. b 6.2.2 Sufficient Condition for a Broken Plug To develop the sufficient condition for a broken plug, we assume T (x ,y) xx b = BS. With this assumption (6.10) becomes I(x )= b rVT&b) fVr(xb) / Jo 10 Bdy = By {x ). T (6.18) b This means that we would satisfy the sufficient condition for plug breaking when there exists an x G (0,1/2) satisfying G(x) = 0, where b G{x) = -Sjn (x)y {x) + B T — - 1 dx - 6By . Jo Vy T When this condition is satisfied, we know that the plug should be broken. 66 (6.19) 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 C r i t e r i a 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 p to break order x which results in a poor approximation for u . The poor approximation for u then causes a p p 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 0.1 0.2 0.3 0.4 0.5 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 \y — VT\)y for the underestimation of the sufficient condition for plug breaking. 69 This could account 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 [4] R. Glowinski. Numerical New York, 1984. to Fluid Methods Dynamics. for Nonlinear Cambridge University Press, 1967. Variational Problems. Springer-Verlag, [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., 4 (%), 5 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
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Slow flow of a Bingham fluid in a gap of slowly varying...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Slow flow of a Bingham fluid in a gap of slowly varying width Ryan, Daniel Padraigh 2003
pdf
Page Metadata
Item Metadata
Title | Slow flow of a Bingham fluid in a gap of slowly varying width |
Creator |
Ryan, Daniel Padraigh |
Date Issued | 2003 |
Description | 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. |
Extent | 3207439 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-03 |
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.0080067 |
URI | http://hdl.handle.net/2429/14634 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2003-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2003-0602.pdf [ 3.06MB ]
- Metadata
- JSON: 831-1.0080067.json
- JSON-LD: 831-1.0080067-ld.json
- RDF/XML (Pretty): 831-1.0080067-rdf.xml
- RDF/JSON: 831-1.0080067-rdf.json
- Turtle: 831-1.0080067-turtle.txt
- N-Triples: 831-1.0080067-rdf-ntriples.txt
- Original Record: 831-1.0080067-source.json
- Full Text
- 831-1.0080067-fulltext.txt
- Citation
- 831-1.0080067.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-0080067/manifest