Theoretical and numerical study of free-surface flowof viscoplastic fluids: 2D dambreaks, axisymmetricalslumps and surges down an inclined slopebyYe LiuB.Sc, Mathematics, Peking University, 2012M.Sc, Mathematics, University of British Columbia, 2014A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Mathematics)The University of British Columbia(Vancouver)May 2019c© Ye Liu, 2019The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled:Theoretical and numerical study of free-surface flow of viscoplas-tic fluids: 2D dambreaks, axisymmetrical slumps and surgesdown an inclined slopesubmitted by Ye Liu in partial fulfillment of the requirements for the degree ofDoctor of Philosophy inMathematics.Examining Committee:Neil Balmforth, MathematicsSupervisorSarah Hormozi, Mechanical EngineeringCo-supervisorAnthony Wachs, Mathematics and Chemical & Biological EngineeringSupervisory Committee MemberJames J. Feng, Chemical and Biological EngineeringAdditional ExaminerChristian Schoof, Earth, Ocean and Atmospheric SciencesAdditional ExaminerAdditional Supervisory Committee Members:Anthony WachsSupervisory Committee MemberBrian WettonSupervisory Committee MemberiiAbstractThe dynamics of free surface flow of yield stress fluid under gravity has been anopen problem, both in theory and in computation. The contribution of this thesiscomes in three parts.First, we report the results of computations for two dimensional dambreaksof viscoplastic fluid, focusing on the phenomenology of the collapse, the modeof initial failure, and the final shape of the slump. The volume-of-fluid method(VOF) is used to evolve the surface of the viscoplastic fluid, and its rheology iscaptured by either regularizing the viscosity or using an augmented-Lagrangianscheme. The interface is tracked by the Piecewise-Linear-Interface-Calculation(PLIC) scheme, modified in order to avoid resolution issues associated with theover-ridden finger of ambient fluid that results from the no slip condition and theresulting inability to move the contact line. We establish that the regularized andaugmented-Lagrangian methods yield comparable results. The numerical resultsare compared with asymptotic theories valid for relatively shallow or verticallyslender flow, with a series of previously reported experiments, and with predic-tions based on plasticity theory.Second, we report computations of the axisymmetric slump of viscoplasticfluid, with the PLIC scheme improved for mass conservation. The critical yieldstress for failure is computed and bounded analytically using plasticity methods.The simulations are compared with experiments either taken from existing liter-ature or performed using Carbopol. The comparison is satisfying for lower yieldstresses; discrepancies for larger yield stresses suggest that the mechanism of re-iiilease may affect the experiments.Finally, we report asymptotic analyses and numerical computations for surgesof viscoplastic fluid down an incline with low inertia. The asymptotic theoryapplies for relatively shallow gravity currents. The anatomy of the surge consistsof an upstream region that converges to a uniform sheet flow, and over which atruly rigid plug sheaths the surge. The plug breaks further downstream due to thebuild up of the extensional stress acting upon it, leaving instead a weakly yieldedsuperficial layer, or pseudo-plug. Finally, the surge ends in a steep flow front thatlies beyond the validity of shallow asymptotics.ivLay SummaryYield-stress fluids are common in our daily life: toothpaste, ketchup, polymergels, fresh concrete, etc. The common feature of these fluids is that they all tendto form certain shapes, when poured on a plane, or squeezed out of a tube. Unlikewater, which behaves totally like liquid and deforms with even a tiny bit of force,yield-stress fluid behaves like solid, as long as the force applied on it is withincertain threshold, called the yield-stress.So what is the final shape of a can of ketchup when it is poured on a plane andspreads out? What is the stress state of that slump? In which condition will theketchup not deform at all? How it is related to practical application? In this study,we will answer these questions in a systematic manner.vPrefaceIn this section, we briefly explain the contents of the journal papers that are pub-lished or submitted for publication from this thesis and clarify the contributionsof co-authors in the papers.• [Liu, Y.], Balmforth, N.J. & Hormozi, S. & Hewitt D.R. (2016) Two-dimensional viscoplastic dambreaks. J. Non-Newtonian Fluid Mech.238, 65-79.This publication has mostly focused at trying to understand the flow dy-namics of yield stress fluid in 2D geometry. Chapter 2 includes the contentsof this publication. The author of this thesis was the principal contributorto this publication. Professor Sarah Hormozi assisted with code develop-ment. Dr. Hewitt assisted with the numerical computation. Professor NeilBalmforth supervised the research and assisted with writing the paper.• [Liu, Y.], Balmforth, N.J. &Hormozi, S. (2018) Axisymmetric viscoplas-tic dambreaks and the slump test. J. Non-Newtonian Fluid Mech. 258,45-57This publication has mostly focused at trying to understand the fluid dy-namics of yield stress fluid in axisymmetric geometry. Chapter 3 includesthe contents of this publication. The author of this thesis was the princi-pal contributor to this publication. Professor Sarah Hormozi assisted withcode development. Professor Neil Balmforth supervised the research andassisted with writing the paper.vi• [Liu, Y.], Balmforth, N.J. & Hormozi, S. (2018) Viscoplastic surgesdown an incline. J. Non-Newtonian Fluid Mech. Submitted for pub-lication, under review.This publication has mostly focused at trying to understand the flow dy-namics on a slope in 2D geometry. Chapter 4 includes the contents of thispublication. The author of this thesis was the principal contributor to thispublication. Professor Sarah Hormozi assisted with code development. Pro-fessor Neil Balmforth supervised the research and assisted with writing thepaper.viiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Objectives of the thesis . . . . . . . . . . . . . . . . . . . . . . . 21.2 Viscoplastic fluids . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Mathematical Tools . . . . . . . . . . . . . . . . . . . . . . . . . 41.3.1 Lubrication theory . . . . . . . . . . . . . . . . . . . . . 41.3.2 Lubrication paradox . . . . . . . . . . . . . . . . . . . . . 51.3.3 Slipline model . . . . . . . . . . . . . . . . . . . . . . . . 61.3.4 Numerical tools . . . . . . . . . . . . . . . . . . . . . . . 7viii1.3.5 Volume of fluid . . . . . . . . . . . . . . . . . . . . . . . 81.3.6 Regularization method . . . . . . . . . . . . . . . . . . . 91.3.7 Augmented-Lagrangian method . . . . . . . . . . . . . . 121.3.8 PLIC advection scheme . . . . . . . . . . . . . . . . . . . 131.3.9 Treatment of moving contact lines . . . . . . . . . . . . . 151.4 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . 161.4.1 Two-dimensional viscoplastic dambreaks . . . . . . . . . 161.4.2 Axisymmetric viscoplastic dambreaks and the slump test . 181.4.3 Viscoplastic surges down an incline . . . . . . . . . . . . 181.5 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . 192 Two dimensional viscoplastic dambreak . . . . . . . . . . . . . . . . 212.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.1.1 Dambreak arrangement and solution strategy . . . . . . . 222.1.2 Model equations . . . . . . . . . . . . . . . . . . . . . . 232.2 Newtonian benchmark . . . . . . . . . . . . . . . . . . . . . . . 262.3 Bingham slumps . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.3.1 Slump and plug phenomenology . . . . . . . . . . . . . . 302.3.2 Shallow flow . . . . . . . . . . . . . . . . . . . . . . . . 352.3.3 Slender slumps . . . . . . . . . . . . . . . . . . . . . . . 402.3.4 Failure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 412.3.5 The final shape and slump statistics . . . . . . . . . . . . 482.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . 533 Axisymmetric viscoplastic dambreaks and the slump test . . . . . . 553.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 563.1.1 Problem set-up and solution strategy . . . . . . . . . . . . 563.1.2 Dimensionless model equations . . . . . . . . . . . . . . 573.1.3 PLIC scheme with interface correction . . . . . . . . . . . 593.2 Bingham slumps . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.2.1 Shallow flow . . . . . . . . . . . . . . . . . . . . . . . . 63ix3.2.2 Slender columns . . . . . . . . . . . . . . . . . . . . . . 673.2.3 Failure mode . . . . . . . . . . . . . . . . . . . . . . . . 693.3 Comparison with experiments . . . . . . . . . . . . . . . . . . . 743.3.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 743.3.2 Cylindrical dambreaks . . . . . . . . . . . . . . . . . . . 753.3.3 Extrusions . . . . . . . . . . . . . . . . . . . . . . . . . . 853.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . 894 Surges of viscoplastic fluids on slopes . . . . . . . . . . . . . . . . . 924.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 934.1.1 Model equations . . . . . . . . . . . . . . . . . . . . . . 934.1.2 The sheet-flow solution . . . . . . . . . . . . . . . . . . . 954.1.3 Volume-of-fluid computations . . . . . . . . . . . . . . . 954.2 Asymptotic analysis . . . . . . . . . . . . . . . . . . . . . . . . . 974.2.1 Plugged Flow . . . . . . . . . . . . . . . . . . . . . . . . 984.2.2 Breaking the plug . . . . . . . . . . . . . . . . . . . . . . 1004.2.3 Lubrication Zone . . . . . . . . . . . . . . . . . . . . . . 1014.2.4 Improving the lubrication theory . . . . . . . . . . . . . . 1024.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . 1054.3.1 Comparison with asymptotics . . . . . . . . . . . . . . . 1054.3.2 Comparison with experiments . . . . . . . . . . . . . . . 1114.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . 1165 Summary and future research directions . . . . . . . . . . . . . . . 1185.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1185.2 Limitations of the study . . . . . . . . . . . . . . . . . . . . . . . 1195.3 Summary of contributions . . . . . . . . . . . . . . . . . . . . . . 1205.4 Future research directions . . . . . . . . . . . . . . . . . . . . . . 121Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123xA Two-dimensional viscoplastic dam break . . . . . . . . . . . . . . . 131A.1 Further numerical notes . . . . . . . . . . . . . . . . . . . . . . . 131A.1.1 Parameter settings and other details . . . . . . . . . . . . 131A.1.2 The failure computation for Re = t = 0 . . . . . . . . . . . 132A.1.3 Thickness of the over-ridden finger . . . . . . . . . . . . . 133A.2 Shallow flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134A.3 Slender columns . . . . . . . . . . . . . . . . . . . . . . . . . . . 136A.4 Plasticity solutions . . . . . . . . . . . . . . . . . . . . . . . . . 138A.4.1 Slipline fields . . . . . . . . . . . . . . . . . . . . . . . . 138A.4.2 Simple failure modes . . . . . . . . . . . . . . . . . . . . 141B Axisymmetric viscoplastic dambreaks and the slump test . . . . . . 144B.1 Numerical convergence study . . . . . . . . . . . . . . . . . . . . 144B.2 Higher-order shallow asymptotics . . . . . . . . . . . . . . . . . 147B.3 Shallow slippy flows . . . . . . . . . . . . . . . . . . . . . . . . 149B.4 Tall slender cones . . . . . . . . . . . . . . . . . . . . . . . . . . 150C Surges of viscoplastic fluids on slopes . . . . . . . . . . . . . . . . . 152C.1 Additional numerical details . . . . . . . . . . . . . . . . . . . . 152C.1.1 Resolution study . . . . . . . . . . . . . . . . . . . . . . 153C.1.2 Inertial effects . . . . . . . . . . . . . . . . . . . . . . . . 154C.1.3 The effect of the back wall . . . . . . . . . . . . . . . . . 154C.2 Improved lubrication solution . . . . . . . . . . . . . . . . . . . . 157C.3 Improved slump profiles . . . . . . . . . . . . . . . . . . . . . . 160xiList of TablesTable 4.1 Experimental parameters from [36]. . . . . . . . . . . . . . . . 110xiiList of FiguresFigure 1.1 A sketch of the slipline field . . . . . . . . . . . . . . . . . . 7Figure 2.1 A sketch of the geometry for the case of a rectangular initialblock. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Figure 2.2 Snapshots of the evolving interface for a Newtonian fluid . . . 27Figure 2.3 Flow front X(t) plotted against time and interface shape forcomputations with Newtonian fluids . . . . . . . . . . . . . . 28Figure 2.4 Front position and central depth for Bingham dambreaks witha square initial block . . . . . . . . . . . . . . . . . . . . . . 31Figure 2.5 Snapshots of a collapsing square . . . . . . . . . . . . . . . . 32Figure 2.6 Snapshots of the evolving interface for a triangular initial con-dition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33Figure 2.7 Snapshots of the evolving interface for an incomplete slump . 35Figure 2.8 Slumps of slender columns . . . . . . . . . . . . . . . . . . . 36Figure 2.9 Comparison of shallow-layer theory with numerical results fora collapsing square block . . . . . . . . . . . . . . . . . . . . 38Figure 2.10 Critical yield stress, Bc, plotted against initial width X0 . . . . 42Figure 2.11 Trial velocity fields to compute lower bounds on Bc . . . . . . 44Figure 2.12 Strain-rate invariant plotted logarithmically as a density on the(x,z)−plane for blocks . . . . . . . . . . . . . . . . . . . . . 46Figure 2.13 Strain-rate invariant plotted logarithmically as a density on the(x,z)−plane for triangles . . . . . . . . . . . . . . . . . . . . 47xiiiFigure 2.14 Profiles of the final deposit . . . . . . . . . . . . . . . . . . . 49Figure 2.15 Final numerical solution for the slump of an initial square andsliplines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49Figure 2.16 Scaled final runout and central depth as a function of B/√X0for Bingham slumps from square initial conditions . . . . . . 51Figure 2.17 Scaled final runout and central depth as a function of B/√X0for Bingham slumps from rectangular and triangular initialconditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52Figure 3.1 A sketch of the geometry for the case of a cylindrical block . . 56Figure 3.2 Time series of flow height and front position for Binghamslumps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 3.3 Profiles of the final deposit of cylindrical slumps . . . . . . . 62Figure 3.4 Evolution of the interface for a slumping cylinder with (R,B)=(1,0.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63Figure 3.5 Final profile for a cylindrical dambreak with B = 0.0074 andR = 0.2546 . . . . . . . . . . . . . . . . . . . . . . . . . . . 66Figure 3.6 Collapse of a cylindrical column . . . . . . . . . . . . . . . . 69Figure 3.7 Critical yield stress for collapse Bc as a function of initial ra-dius R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70Figure 3.8 Strain-rate invariant plotted logarithmically as a density on the(r,z) plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure 3.9 Comparison of experimental slump height with previously pub-lished data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 3.10 Experimental and numerical slump heights against B . . . . . 77Figure 3.11 The central depth and radius of the final deposit against Bing-ham number . . . . . . . . . . . . . . . . . . . . . . . . . . . 80Figure 3.12 Side images from the slump experiments and theoretical finalprofiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82Figure 3.13 The dimensionless slump height for cylindrical dambreaks withR = 1, varying B and the Reynolds numbers indicated . . . . . 84xivFigure 3.14 Maximum radius versus depth for experimental extrusions overthe rough and smooth plexiglass . . . . . . . . . . . . . . . . 86Figure 3.15 Images and profiles from the extrusion experiments . . . . . . 88Figure 3.16 Slump heights for Bingham dambreaks with an initially coni-cal shape corresponding to the ASTM standard . . . . . . . . 91Figure 4.1 Sketch of the geometry of the surge in the frame of referencein which there is no net flux and the flow is steady. The freesurface is located at z = h; the level z = Y divides a fullyyield region underneath from either a true plug or a weaklyyielded pseudo-plug. The flow divides into three regions: aplugged flow region (PF) where the superficial layer of fluidis not yielded, a lubrication zone (LZ) where the pseudo-plugarises, and the flow front (FF) where the dynamics is not shallow. 93Figure 4.2 Asymptotic surge solutions for (a) the plastic limitY → 0, and(b) Y∞ = 1−S−1B = 0.2. The solid lines show the improvedlubrication solutions for h(x) plotted against S(x− x f ) with12εpiB = 0.2 (as given by either equations (4.41) and (4.42));the dotted lines indicate the predictions of the leading-ordertheory. In (b), the corresponding fake yield surfaces Y (x) arealso plotted. The insets show magnifications near the flow front.104Figure 4.3 Numerical solution for θ = 10◦ and Y∞ = 0.8 showing (a)τI, (b) τxz and (c) τxx. The darker (blue) lines show samplestreamlines and the dotted white line is the true yield surfacewhere τI = B. The dashed line shows the contour level whereτxz = B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106xvFigure 4.4 Flow profiles and fake yield surfaces for (a) numerical simu-lations, (b) leading-order lubrication theory and (c) improvedasymptotic theory (using (4.42)), with θ = 10◦ and Y∞ = 0.4,0.6 and 0.8. The solid and dashed lines show h and Y , with Ydefined as the contour level where τxz = B for the simulations.The flow profiles deepen with increasing Y∞. The shaded re-gions in (a) show the true plugs (which extend up to the freesurface in each case, and with the shading darkening with in-creasing Y∞). . . . . . . . . . . . . . . . . . . . . . . . . . . 107Figure 4.5 Numerical solution for Y∞ = 0.8 and θ = 5◦ (B = 0.0175),showing a density plot of τI with z = h(x) and the fake yieldsurface superposed; the true plug is shaded black. The dashedlines show the prediction of the leading-order lubrication the-ory (with the flow front aligned). The dotted lines show thepredictions of the plugged-flow solution in (4.26), with 1−hmatched to the simulation at x = 5. . . . . . . . . . . . . . . . 109Figure 4.6 Experiment C5 of [36], showing (a) u(x,z) and (b) w(x,z) asdensities over the (x,z)−plane, for a qualitative comparisonwith their figure 6. A selection of streamlines is also shown(thinner blue lines). The dashed lines shows the levels whereu = 0 and z =Y (τxz = B). In leading-order lubrication theory,u = 0 along z =Y −Y [nY/h(2n+1)]n/(n+1); this prediction isalso drawn as the lighter (pink) solid line. . . . . . . . . . . . 110Figure 4.7 Experiment C2 of [36]: shown are (a) log10(γ˙/0.26), (b) thesurface velocity, and (c)–(h) velocity profiles at the x−positionsindicated in (b), for a comparison with figure 10 in [36]. The(red) dashed lines indicate the predictions of the leading-orderlubrication theory. In (a) the solid white line shows the con-tour τxz = B, whereas the dashed white line is the fake yieldsurface z = Y (x) of the leading-order lubrication theory. . . . 113xviFigure 4.8 A similar picture to figure 4.7, but for Experiment C5 of [36]and for comparison with their figure 9 (with log10(γ˙/0.38)plotted in (a)) . . . . . . . . . . . . . . . . . . . . . . . . . . 114Figure 4.9 Flow profiles for a simulation with θ = 14.6◦, τY = 6 Pa, κ =6.65 Pa sn and n = 0.405 (the experiment C PIV of [23]), for(a) ub = 0.2, 0.4, 0.8, 1.2 and 1.4 m/s, in the steady state,and (b) ub = 1.6 m/s at a succession of times, starting fromthe initial profile shown by the dashed line (t = 0, 0.06, 0.08,0.11, 0.17 and 0.19s). In (a), the profile for ub = 1.4 m/s isshown by the solid line and several streamlines are also plotted(dotted lines). . . . . . . . . . . . . . . . . . . . . . . . . . . 116Figure A.1 Slender asymptotic solutions . . . . . . . . . . . . . . . . . . 139Figure A.2 Slipline solutions for a rectangle . . . . . . . . . . . . . . . . 141Figure B.1 Computations with varying grid size . . . . . . . . . . . . . . 145Figure B.2 Resolution studies for various numerical solvers . . . . . . . . 146Figure B.3 Scaled stress components and radial speed . . . . . . . . . . . 148xviiFigure C.1 Simulations of a Bingham surge profile with θ = 10◦ andB = 0.0522 (Y∞ = 0.7, fluid area of 24, ub = 0.041 and Re =1) with the grid sizes indicated. (a) shows h and Y ; (b) amagnification near the flow front; (c) contours of constantτxz/B, as indicated (with the surge profile of the finest res-olution case shown by the lighter grey line); (d) the shearstress along z = 0.025 (the lowest grid cell of the coarsestcomputation). In (a), the inset shows the flow front X f , meanshear stress, τxz/B =∫ ∫(τxz/B) c dx dz, and τxz at the point(x,z) = (18,0.4), all plotted against the vertical grid spacing∆z (∆x = 4∆z). The convergence of τxz is impeded by theneed to resolve the sharply localized region of high stress un-derneath the flow front (cf. figure 4.3 and panels (c) and (d)).. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155Figure C.2 Simulations of Bingham surges for θ = 10◦, ub = 0.041 andB = 0.0522 (Y∞ = 0.7) for the Reynolds numbers indicated(the total area of fluid is 24). Plotted are h, Y (τxz = B) and thestress levels τI/B = 1, 2 and 3. . . . . . . . . . . . . . . . . . 156Figure C.3 Simulations for θ = 10◦ and B = 0.0354 (Y∞ = 0.8) with dif-ferent left-hand boundary conditions: (a) (u,w) = (usheet ,0)and (b) (u,w) = (0,0). Shown is the second invariant τI asa density on the (x,z)−plane. The (true) yield surfaces areindicated by the solid (green) lines. The simulation in (a) cor-responds to the truncated solution shown in figure 4.3. Theprofiles and the true and fake yield surfaces are compared in(c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162xviiiFigure C.4 Slump profiles for (a)-(b) S = 0 and (b)-(c) S 6= 0. For (a)-(b),scaled variables are plotted, which eliminates any free param-eters; in (c)-(d) the profiles are shown for 12εpiB = 0.2. Thedark solid lines show the solutions to (4.38), whereas the dot-ted lines show the solutions to the iterated version in(4.40); thedashed lines show the leading-order approximation. In (a)–(b)the lighter (red and blue) lines show the results of a series ofsimulations from [53] with ε = 1, B = 0.02, 0.03, ..., 0.1 andthe flow fronts aligned. In (c)–(d) the lighter (red) lines showadditional simulations of the slump of a rectangular block onan incline with B = (0.1,0.2,0.3,0.4)/pi and θ = 5◦; the insetshows the aligned, but unscaled profiles. . . . . . . . . . . . 163xixAcknowledgmentsI would like to convey my gratitude to all people who gave me the possibility tocomplete this dissertation.In the first place, I would like to express my sincere gratitude to my supervi-sor, Professor Neil Balmforth for the continuous support of my Ph.D study andresearch, for his patience, motivation, enthusiasm, and immense knowledge. Hiswide knowledge and his logical way of thinking have been of great value to me.I am deeply grateful to my co-supervisor Sarah Hormozi for her guidance andassistance in my research. I would not be able to accomplish the development ofthe computational tool without her teaching and help.I gratefully acknowledge Professor Anthony Wachs, Brian Wetton and JamesFeng, for their supervision, advice, and guidance from the early stage of this re-search. Above all and the most needed, they provided me support and friendlyhelp in various ways.I would like to thank all my colleagues and schoolmates who have supportedme all along my study and research.Most importantly, I would like to thank my entire family. My immediate fam-ily to whom this dissertation is dedicated, has been a constant source of love,concern, support and strength all these years. I would like to express my heart-feltgratitude to all my friends who have supported and encouraged me throughout thisendeavor.xxDedicationTo my familyxxiChapter 1IntroductionViscoplastic fluids are commonly encountered in natural settings in geophysics(e.g. mud and and lava) and biology (mucus and blood clots), and feature inmany engineering processes in, for example, the food (fruit pulp, dairy productsand chocolate confections) and petroleum industries (drilling mud, cement andwaxy crude oil). These materials flow like viscous fluid once the stress exceedsa certain threshold (the yield stress), while remaining solid-like, with any defor-mation often assumed small and discarded. In that situation, which underpinspopular constitutive models such as the Bingham and Herschel-Bulkley laws, thestress state of the material becomes formally indeterminate [14]. Together withthe need to track the boundaries of the yielded regions, this complicates signifi-cantly efforts aimed at theoretical modelling. Numerical strategies to overcomesuch difficulties have been developed in recent years and here we apply them tothe particular problem of the free-surface flow of viscoplastic fluid under gravity,including two dimensional dambreaks, axisymmetric slumps and surges onan inclined slope, which feature in a diverse range of problems from geophysicsto engineering. These flows can constitute natural or manmade hazards, as in thedisasters caused by mud surges and the collapse of mine tailing deposits. In anindustrial setting, the controlled release of a reservoir in a simple dambreak exper-iment forms the basis of a number of practical rheometers, including the slump1test for concrete [61, 68] and the Bostwick consistometer of food science [12].The Bostwick consistometer, in which materials such as ketchup are released in arectangular channel, is conveniently idealized by a two-dimensional flow (Chap-ter 2), while the slump test features the release of a cylinder of yield-stress fluid,which is widely used as a measurement of concrete rheology (Chapter 3). Finally,viscoplastic surges on an inclined plane (Chapter 4) have been studied experimen-tally, which provides useful results for theoretical comparisons.1.1 Objectives of the thesisThis thesis aims at the three topics described above. Overall our research seeks tofully understand the dynamics of fluid in each context. Different objectives of thisstudy include the following.1. We want to explore or improve the theoretical model in these problems.The lubrication model, widely used in simulating those processes, could beexpanded to higher order to provide more accurate results.2. We want to provide robust and accurate numerical simulations of theseproblems. The difficulties here mainly include the treatment of indeter-minate stress state in the plug and the advection of fluid interface.3. We want to understand if there is a true plug or a plug-like region held closeto, but above, the yield stress, in each case, which has been confusing andcontroversial in previous studies. We want to show it both asymptoticallyand numerically.4. By the theoretical and numerical tools we have developed, we want to per-form various numerical simulations of existing experiments in literature.We have also conducted some experiments ourselves, in order to supportthe improved model and numerical strategies.21.2 Viscoplastic fluidsBefore introducing the concept of viscoplastic fluid, we first review the rheologyof Newtonian fluid, where the shear stress tensor τi j(u) and shear rate tensorγ˙i j(u)are proportional,τi j(u) = µγ˙i j(u)where µ is the viscosity of the material(constant for Newtonian fluids). This equa-tion is called the constitutive law. Many fluids in reality, however, tend to followmore complicated constitutive laws, and are called non-Newtonian fluids.Viscoplastic fluid is a specific type of non-Newtonian fluid, which remainsundeformed until the stress exceeds certain value, called yield stress(τY ). Whenthe stress exceeds the yield stress, then the fluid deforms as viscous fluid, withcertain rheological models, as generated from varying materials, among whichtwo common ones are what we mainly discussed in this research.• Bingham modelThe Bingham model is the simplest viscoplastic model. Above the yieldstress, the constitutive law followsγ˙ jk = 0, τ < τY ,τ jk =(µ +τYγ˙)γ˙ jk, τ > τY ,(1.1)• Herschel-Bulkley modelThe Herschel-Bulkley model extends the simple power-law model to in-clude a yield stress as belowγ˙ jk = 0, τ < τY ,τ jk =(µγ˙n−1+τYγ˙)γ˙ jk, τ > τY ,(1.2)where τ =√12 ∑ j,k τ2jk and γ˙ =√12 ∑ j,k γ˙2jk. There are other models but we will3focus on the study of these two. The reason for studying the Bingham model isbecause of its simplicity, which reveals the yield stress behavior with minimumnumber of parameters considered. And the reason we use the Herschel-Bulkleymodel is because the most commonly used experiment material, Carbopol, is wellfitted by this model so that we can perform numerical simulations of existingexperiments with Carbopol. We believe that the extension to other yield stressmodels is straightforward, both theoretically and numerically.It is seen from the constitutive laws above that, below the yield stress, the rela-tion between the stress and strain rate remains unknown, which greatly increasesthe difficulty of the analysis and the computation. In the thesis, we will discuss indetail how we overcome this difficulty.1.3 Mathematical ToolsIn this section, we will introduce the main mathematical tools we have appliedin this thesis. For the theoretical part, we shall review the lubrication theory andthe slipline model in plasticity theory. For numerical part, we will describe thebasic idea of the regularization method, augmented-Lagrangian method and thePiecewise-Linear-Interface-Calculation (PLIC) method.1.3.1 Lubrication theoryIn fluid dynamics, lubrication theory describes the flow of fluids (liquids or gases)in a geometry in which one dimension is significantly smaller than the others. Thishelps simplify the momentum equation and allows for consistent asymptotic ex-pansions. Moreover, lubrication theory applies when the flow is slow so that onemay add the Stokes approximation. For free-surface viscoplastic flows, the appli-cation of lubrication theory dates back to Liu and Mei [51]. Further theoreticalinsights and extensions were given by Balmforth et al. [6, 7, 10, 11, 15], Matsonand Hogg [43, 58], and Ancey and Cochard [2], based on which the leading order4solution for the final shape of a viscoplastic slump has been well known ash(x) =√2τYρg(x∞− x)in two-dimensional geometry, andh(r) =√2τYρg(r∞− r)in axisymmetric geometry, where h(x),h(r) refer to the flow height as a function ofthe horizontal coordinate, ρ is the fluid density, g is the gravitational acceleration,and x∞,r∞ are the maximum length of the slump.In this thesis, we will apply lubrication theory to all the three problems con-sidered, and improve those solutions to higher order.1.3.2 Lubrication paradoxThe so-called “lubrication paradox” is from Lipscomb & Denn [50], who appliedlubrication theory and predicted plug regions that were apparently below the yieldstress, while deformation rate is not zero. This inconsistency was later understood,as the paradox can be resolved by recognizing that the lubrication theory is theleading order term of an asymptotic expansion; an exploration of higher-orderterms demonstrates that the problematic plugs are actually slightly above the yieldstress, allowing for their deformation (Putz et al. [65], Walton & Bittleston [81]).The region has been termed a pseudo-plug to emphasize its nature, with the borderreferred to as a fake yield surface. Therefore there is no lubrication paradox.In the pseudo-plug, the stress equals the yield stress to the leading order, sothat it can be considered as a perfect plastic material, where the stress componentmust satisfyτ2xx + τ2zz +2τ2xz = 2τ2Y5in two-dimension, orτ2rr + τ2θθ + τ2zz +2τ2rz = 2τ2Yin axisymmetric geometry. This is the key assumption based on which the leadingorder asymptotic solution can be expanded to higher order. (And this is actuallytrue for slump problem in the steady state, even when it is not shallow. We willsee it in the following chapters.)1.3.3 Slipline modelSlipline field theory is used to model plastic deformation in plane strain only for asolid that can be represented as a rigid-plastic body. Elasticity is not included andthe loading has to be quasi-static. In terms of applications, the approach has beenapplied by Nye [59] to study the “dynamics” of glaciers, and Dubash et al. [33]to study the final shape of viscoplastic slump. In 2D geometry, the equilibriumequations for plane strain are:∂σxx∂x+∂τxz∂ z= 0∂τxz∂x+∂σzz∂ z= 0where σxx,σzz are the total stress components. There will always be two per-pendicular directions of maximum shear stress in a plane. These generate twoorthogonal families of slip lines called α lines and β lines, as shown in 1.1. Thevelocity distribution and stress state in the solid can always be determined fromthe geometry of these lines, and the Hencky relations givep+2kφ = constant along an α linep−2kφ = constant along a β line6Figure 1.1: A sketch of the slipline fieldwhere p is the hydrostatic pressure, k is the shear yield stress(τY in our case), andφ is the angle formed by the x-axis and the α line. Andσxx =−p− k sin2φσzz =−p+ k sin2φτxz = kcos2φBased on the relations above, we are able to construct a slipline field, given ap-propriate initial conditions. This technique has been applied in chapter 2 to solvethe final shape and stress state of two dimensional viscoplastic slump.1.3.4 Numerical toolsIn this thesis, I use the volume-of-fluid method (VOF) to solve the problem. Theadvantage of the VOF method is mass is well conserved. In our study, we usestructured rectangular grids for the simulation. The method works well with our7problem.In every time step, the code does two things:• Solving the Navier-Stokes equations of two fluids, to get the velocity field.• Solving the advection equation by given velocity field, to move the inter-face.There are two options for the Navier-Stokes solver: regularization method oraugmented-Lagrangian method, which are the most commonly used schemes forviscoplastic fluid. The advantage of the regularization method is that the compu-tation converges faster, while the disadvantage is that it is not quite accurate whenτ ∼ τY , so that it is not recommended in studying the incipient failure problem.The advantage of augmented-Lagrangian method is that it solves the exact consti-tutive equations. The disadvantage is that it converges much more slowly than theregularization method.All of the algorithms are implemented in C++ as an application of PELICANS,an object oriented platform developed at IRSN, France, to provide general frame-works and software components for the implementation of PDE solvers. PELI-CANS is distributed under the CeCILL license agreement.1.3.5 Volume of fluidWe use the concept of the volume-of-fluid method to simulate a two-phase fluidproblem. The volume fraction variable c = 1 for the lower layer, and c = 0 forthe upper layer. The problem is discretized into a Navier-Stokes equation and anadvection-diffusion equation. See Equation 1.3Re(c)[∂u∂ t+(u ·∇)u]=−∇p+∇ · τ(u)+Ri(c)g∇ ·u = 0∂c∂ t+∇ · (uc) = 0(1.3)8whereRe(c) :=(c+(1− c)ρ2ρ1)ρ1URµ1Ri(c) :=(c+(1− c)ρ2ρ1)ρ1gR2µ1U0τ(c) := cτ(u)+(1− c)µ2µ1γ(u)where τ(u) is defined in the following section.1.3.6 Regularization methodIn this section, we summarize the regularization method for single-phase fluid,but the generalization to VOF is straightforward. For a general Navier-Stokesequation in a computational domain Ω, with no-slip condition on the boundaryΓD and no-flux condition on the boundary ΓN , as belowρ[∂u∂ t+(u ·∇)u]=−∇p+∇ · τ(u)+ρg in Ω∇ ·u = 0 in Ωu = 0 in ΓD∇u ·n = 0 in ΓN(1.4)At a certain time tn, given the current velocity field un, time step ∆t, the discretescheme for solving un+1 isρ∆tun+1+(ρun ·∇)un+1−∇ · τ(un+1)+∇pn+1 = ρg+ ρ∆tun in Ω∇ ·un+1 = 0 in Ωun+1 = 0 in ΓD∇un+1 ·n = 0 in ΓN(1.5)9The variational form is: Find (u, p)∈ Su×V p so that for all trial functions (v,q)∈V u×V p:1∆t∫Ωρu · v+∫Ω(ρun ·∇)u · v+∫Ωτ(u) : ∇v−∫Ωp∇ · v =∫Ωρg · v+ 1∆t∫Ωρun · v∫Ωq∇ ·u = 0(1.6)where the function spaces Su,V u,V p are defined asSu = {v ∈ H1(Ω)|v = 0 in ΓD}V u = {v ∈ H1(Ω)|v = 0 in ΓD}V p = L2(Ω)(1.7)We can rewrite this in a simple way1∆tm(u,v)+ c(un;u,v)+ k(u,v)+b(v, p) = f (v)b(u,q) = 0(1.8)wherem(u,v) =∫Ωρu · vc(un;u,v) =∫Ωρ(un ·∇)u · vk(u,v) =∫Ωτ(u) : ∇vb(v, p) =−∫Ωp∇ · vf (v) =∫Ωρg · v+ 1∆tm(un,v)(1.9)10This variational equation is then solved using a Galerkin method. On a finitedimensional subspace V uh ⊂V u,V ph ⊂V p1∆tm(uh,vh)+ c(unh;uh,vh)+ k(uh,vh)+b(vh, ph) = f (vh)b(uh,qh) = 0(1.10)After a finite element basis is chosenXh = span{ψuk | 0≤ k < Nudo f }Mh = span{ψ pk | 0≤ k < N pdo f }The variational equation is transformed to a linear system{AU +BT P = FBU = 0(1.11)whereAIJ =1∆tm(ψuJ ,ψuI )+ c(unh;ψuJ ,ψuI )+ k(ψuJ ,ψuI ) 0≤ I,J < Nudo fBIJ = b(ψuJ ,ψpI ) 0≤ I < N pdo f ,0≤ J < Nudo fFI = f (ψuI ) 0≤ I,J < Nudo f(1.12)Notice that k(ψuJ ,ψuI ) is a nonlinear term becauseτ(u) =(µ +τy|∇u|)∇u |τ(u)|> τy∇u = 0 |τ(u)|< τy(1.13)Therefore we use a regularized model and an iteration scheme to overcome thisdifficulty.τ(u) =(µ +τy|∇u|+ ε)∇u (1.14)11is a simple regularized model of Eq 1.13. If we assume |∇u| on the denominatoris known, i.e.τ(u) =(µ +τy|∇uˆ|+ ε)∇u (1.15)where uˆ is from the previous iteration step, then the operator is linear. Thus wecan do the iteration on each time step t = tn. The implementation of this is inAlgorithm 1.Algorithm 1 Regularization methodRequire: GivenUn,Pn, Uˆ =Un, Imax, i = 0Generate matrix A,B and vector F from Un,Pn,Uˆwhile i < Imax doSolve the linear system(A BTB 0)(UP)=(F0)if |Uˆ−U |< tolerance thenreturn Un+1 =UelseUˆ =U Update A from Uˆend ifi = i+1end whilereturn Un+1 =U1.3.7 Augmented-Lagrangian methodA detailed description of the augmented-Lagrangian method can be found in [35].We just list the main ideas here.12Find (u, p,d ,λ ) so that for ∀(v, l,q)ρ∆t∫Ωu · v+∫Ω(ρun ·∇)u · v+∫Ωµ∇u ·∇v−∫Ωp∇ · v+∫Ωλ ·∇v+ r∫Ω(∇u−d) ·∇v =∫Ωf · v∫Ωl∇ ·u = 0d =1r(1− τy|λ + r∇u|)(λ + r∇u) if |λ + r∇u| ≥ τy0 if |λ + r∇u|< τy∫Ω(∇u−d) ·q = 0(1.16)where f = ρg+ ρ∆t un. The implementation of the scheme is in Algorithm 2.1.3.8 PLIC advection schemeWe used the PLIC (short for Piecewise Linear Interface Calculation) to solve theadvection of the fluid volumeC,∂C∂ t+u ·∇C = 0The scheme is designed to avoid numerical diffusion of the interface, and to in-crease accuracy in tracking the interface. It is 2nd order accurate in space. Theequation is discretized by a forward Euler schemeCn+1 =Cn−un∆t ·∇CnThe scheme is decomposed into x and y direction, so we have two sub-steps inevery time step. The key part of this scheme is to compute the line approximationgiven different circumstances. For the details of the this method, we refer toRenardy[57]. An original implementation of PLIC scheme is in Algorithm 3.13Algorithm 2 Augmented-Lagrangian schemeRequire: Given un,r, Imaxi = 0 λ i = 0 d i = 0while i < Imax doi = i+1Solve the problem, find (ui, pi) so that ∀(v, l)ρ∆t∫Ωui · v+∫Ω(ρun ·∇)ui · v+∫Ω(µ + r)∇ui ·∇v−∫Ωpi∇ · v=−∫Ω(λ i−1− rd i−1) ·∇v+∫Ωf · v∫Ωl∇ ·ui = 0Update the velocity gradient dd i =1r(1− τy|λ i−1+ r∇ui|)(λ i−1+ r∇ui) if |λ i−1+ r∇ui| ≥ τy0 if |λ i−1+ r∇ui|< τyUpdate the Lagrangian multiplierλ i = λ i−1+ r(∇ui−d i)if max{|ui−ui−1|, |∇ui−d i|}< tolerance thenreturn un+1 = uiend ifend whilereturn un+1 = u14Algorithm 3 PLIC scheme 1.0Require: Given un,C,∆tEnsure: The CFL condition (udx+vdy)∆t <12Construct the linear approximation y = mx+b in each grid cell, where m = ∇Cand the volume fraction in the grid equals CCompute the new line y = mˆx+ bˆ after advection of un∆tCompute the flux from y = mx+b and y = mˆx+ bˆ and updateCRegularize C so that 0≤C ≤ 11.3.9 Treatment of moving contact linesA common problem in the computation of a two-phase fluid problem with a no-slip boundary condition is the movement of the contact line (the point where theinterface touches the static wall). Failure to address this problem will lead to anunresolved boundary layer (finger) of ambient fluid on the bottom. Commonly aslip boundary condition is used around the contact line to move it. However, theschemes always depend on the contact angle, or what parameterization is used toallow the contact line to move, while surface tension is neglected in our study.Currently we have implemented a simple correction in the PLIC scheme to avoidthe resolution issue associated with the finger. The algorithm is implemented inAlgorithm 4.We find that this scheme successfully removes the finger, leading to an appar-ent motion of the contact line, though none is in fact happening. Currently wehave tested that the solution is not affected by a slight change of the number from0.99 to 0.9. Moreover, the numerical results achieve much better spacial conver-gence using this correction than the original PLIC scheme. The main problem ofthis scheme is the violation of mass conservation, asC is always increasing duringadvection, which can affect the simulation at later times, particularly in axisym-metric geometry. To make the mass conserved, we make further modification of15Algorithm 4 PLIC scheme 1.1Require: Given un,C,∆tEnsure: The CFL condition (udx+vdy)∆t <12Construct the linear approximation y = mx+b in each grid, where m = ∇C andthe volume fraction in the grid equals CCompute the new line y = mˆx+ bˆ after advection of un∆tCompute the flux from y = mx+b and y = mˆx+ bˆ and updateCRegularize C so that 0≤C ≤ 1Check the value ofC on the bottom layer of the grid cell and the one just aboveit, we note them as C0 and C1, respectivelyIf C1 > 0.99, imposeC0 =C1our scheme, see Algorithm 5. In this way, the mass conservation is guaranteedand the solution converges even better.1.4 Literature review1.4.1 Two-dimensional viscoplastic dambreaksDespite wide-ranging practical applications, the theoretical modelling of viscoplas-tic dambreaks remains relatively unexplored. Asymptotic theories for shallow,slow flow have received previous attention and permit a degree of analytical in-sight into the problem (see [11, 52] and references therein). Numerical compu-tations of two-dimensional dambreaks have also been conducted to model flowsthat are not necessarily shallow [77]. However, these simulations do not provide adetailed survey of the flow dynamics over a wide range of physical conditions andhave focused mainly on determining some of the more qualitative aspects of theend state of a slump, such as its final runout and maximum depth. Complement-ing both asymptotics and numerical simulation are cruder predictions of the final16Algorithm 5 PLIC scheme 1.2Require: Given un,C,∆tEnsure: The CFL condition (udx+vdy)∆t <12Construct the linear approximation y = mx+b in each grid, where m = ∇C andthe volume fraction in the grid equals CCompute the new line y = mˆx+ bˆ after advection of un∆tCompute the flux from y = mx+b and y = mˆx+ bˆ and updateCRegularize C so that 0≤C ≤ 1Check the value of C on the bottom layer of grid and its upper layer, we notethem as C0 and C1, respectivelyIf C1 > 0.99, imposeC0 =C1Compute the mass increment ∆C, and the total volume of the interface CIRescale the volume fraction in the interface by a factor of 1− ∆CCIshape based on solid mechanics and initial failure criteria derived from plasticitytheory [20, 21, 61].The key feature of a viscoplastic fluid that sets the problem apart from a clas-sical viscous dambreak is the yield stress. When sufficient, this stress can hold thefluid up against gravity, preventing any flow whatsoever. If collapse does occur,the yield stress brings the fluid to a final rest and can maintain localized rigid re-gions, or “plugs”, during the slump. The evolving plugs and their bordering yieldsurfaces present the main difficulty in theoretical models, particularly in numeri-cal approaches. Augmented-Lagrangian schemes that deal with the complicationsof the yield stress directly are often time-consuming to run, whereas regulariza-tions of the constitutive law that avoid true yield surfaces introduce their ownissues [38]. For the dambreak problem, difficulties are compounded by the needto evolve the fluid surface and impose boundary conditions such as no-slip on thesubstrate underneath the fluid.171.4.2 Axisymmetric viscoplastic dambreaks and the slump testAxisymmetric collapses are exploited widely to gauge fluid rheology in the con-crete, mineral and food industries. The slump test, for example, is commonly usedto measure the yield stress of fresh concrete. In this test, a container filled withconcrete is lifted to release the material and allow it to spread under gravity; atstoppage, the vertical distance over which the concrete falls, the “slump height”,is measured as an indicator of yield stress. A number of experimental studies havebeen directed at establishing the precise relation of the slump height to the yieldstress for a variety of different types of viscoplastic fluids [25, 39, 61, 70].Theoretical studies of the slump test have been performed using either numer-ical computations, asymptotic analysis suitable for the limit of shallow flow, orestimates and bounds based on plasticity theory [11, 17, 21, 25, 31, 34, 61, 63,67, 67, 69, 71]. The numerical simulations have been conducted using a varietyof numerical techniques, although most of the algorithms employed were not spe-cially designed to capture yield-stress rheology or carefully track fluid interfaces.The current status of the modelling of the slump test is reviewed by [69].1.4.3 Viscoplastic surges down an inclineThe spreading of viscoplastic fluid over an inclined surface has been studied ex-perimentally in a number of previous studies, often with the goal of inferringthe yield stress from steady flows [28, 29, 32] or the shape of a final deposit[30, 33, 46, 60]. The most thorough and recent laboratory studies include thetransient dam-break-type experiments of Ancey and co-workers [3, 16, 27] and aseries of investigations on steady viscoplastic surges on inclined conveyor belts[22, 23, 36].Theoretically, a model for shallow viscoplastic flow based on Reynolds lu-brication theory has been widely used to complement such laboratory studies[10, 11, 51]. In this model, the long, thin flow is composed of a fully shearedregion adjacent to the underlying surface buffered from the free surface by a plug-like zone. Importantly, that zone is not truly rigid, but deforms weakly in the18direction of flow and is plug-like in that the transverse velocity profile is largely in-dependent of depth. This structure is common in many shallow viscoplastic flows[4, 81] and results from the separation of length scales in the directions aligned orperpendicular to flow. The border between the fully yield region and the plug-likezone is therefore not a true yield surface; instead, it is often referred to as a fakeyield surface, and the overlying zone as a pseudo-plug. Despite the fact that lubri-cation theory predicts the appearance of pseudo-plugs in shallow flows, it is alsoknown that genuine rigid plugs can appear within these zones surrounding pointsof symmetry [65, 81] or replace them in flow down almost uniform channels [37].This raises the question of whether the superficial regions of a free-surface flowcan also plug up in the far upstream extent of a steady surge flow, where the flowbecomes almost uniform, as discussed in a qualitative way by Piau [63]. Indeed,for a truly steady surge that extends infinitely far upstream, one expects that theflow converges to a uniform sheet flow which, for a yield-stress fluid, is sheathedby a true plug.1.5 Outline of the thesisThe two dimensional viscoplastic dambreak problem is studied in chapter 2. Nu-merical schemes are validated with a Newtonian fluid problem first, after whicha thorough investigation of Bingham slump is illustrated, including the plug phe-nomenology, the final shape in shallow and slender geometry, and the criticalvalue for the slump to happen, etc. In comparison with the numerical result, animproved asymptotic analysis has been derived and presented. The theoreticaland numerical results are then compared with available experimental data fromliterature.In chapter 3, the axisymmetric slump of viscoplastic fluid is studied, withsimilar structure as in chapter 2. In addition, we conduct experiments of slumpand extrusion tests with Carbopol, and present the results, in comparison withnumerical simulations.In chapter 4, the surge problem of viscoplastic fluids on slopes is studied,19theoretically and numerically. The theory starts with a simple sheet flow solution,where a plug exists everywhere, then perturbation study has implied the existenceof a true plug and its breaking condition. The transition zone from a true plugto sheared flow, or the pseudo-plug is shown then, with an improved lubricationtheory. Finally, numerical simulations are conducted and compared with bothasymptotic results, and experimental results from [21].Chapter 5 of the thesis contains a summary of conclusions of the thesis andrecommendations for future work.Appendix A corresponds to Chapter 2. It provides more details of our parame-ter setting in the two dimensional viscoplastic dambreak problem. It also includesthe mathematical derivation of the improved asymptotic analysis, the main idea ofconstructing the slipline field, and ways to compute the critical value with simplefailure modes.Appendix B corresponds to Chapter 3. It includes a detailed resolution studyof the numerical simulation of axisymmetric dambreak problem. Moreover, asymp-totic analysis in axisymmetric geometry is also illustrated here.Appendix C corresponds to Chapter 4. It consists of a study of the effect ofinertia, or the back wall of the surge. Also a detailed asymptotic analysis is given.Most of the material in the appendices is original work; results arising fromprevious work are identified by citations in each section.20Chapter 2Two dimensional viscoplasticdambreak1In the chapter, we present numerical computations of viscoplastic dambreaksspanning a wide range of physical parameters. Our aim is to describe morefully the phenomenology of the collapse and its plugs, the form of the motionat initiation, and the detailed final shape. Our main interest is in the effect of theyield stress, so we consider a Bingham fluid, ignoring any rate-dependence of theplastic viscosity. We mathematically formulate the dambreak problem in section2.1 and outline the numerical strategies we use for its solution. We use both anaugmented-Lagrangian scheme and regularization of the constitutive law to ac-count for viscoplasticity; to deal with the free surface, we use the volume-of-fluidmethod. The latter method emplaces the viscoplastic fluid beneath a less denseand viscous fluid, then tracks the interface between the two using a concentrationfield. This effectively replaces the single-phase dambreak problem with that of atwo-phase miscible fluid displacement (we ignore surface tension), but introducesa significant complication when imposing a no-slip boundary condition: because1A version of chapter 2 has been published. [Liu, Y.], Balmforth, N.J. & Hormozi, S. & HewittD.R. (2016) Two-dimensional viscoplastic dambreaks. J. Non-Newtonian Fluid Mech. 238, 65-79. [53].21the lighter fluid cannot be displaced from the lower surface, the slumping heav-ier fluid over-rides a shallow finger of lighter fluid which lubricates the overlyingflow and thins continually, leading to difficulties with resolution. We expose thiscomplication for a viscous test case in section 2.2, and identify means to avoidit. We then move on to a discussion of Bingham dambreaks in section 2.3, beforeconcluding in section 2.4. The appendices contain additional technical details ofthe numerical schemes, asymptotic theories for shallow or slender flow, and somerelated plasticity solutions.2.1 Formulation2.1.1 Dambreak arrangement and solution strategyTo simulate the collapse of a Bingham fluid, we consider a two-fluid arrangement,with the yield-stress fluid emplaced underneath a lighter viscous fluid. We ig-nore any interfacial tension. The volume-of-fluid method is used to deal with theboundary between the two fluids: a concentration field c(x,y, t) smooths out andtracks the fluid-fluid interface; c = 1 represents the viscoplastic fluid and the over-lying Newtonian fluid has c = 0. The concentration field satisfies the advectionequation for a passive scalar; no explicit diffusion is included although some isunavoidable as a result of numerical imprecision. Figure 2.1 shows a sketch ofthe geometry; the initial block of viscoplastic fluid has a characteristic height Hand basal width 2L , but we assume that the flow remains symmetrical about theblock’s midline and consider only half of the spatial domain.To deal with the yield stress of the viscoplastic fluid, we use both an augmented-Lagrangian scheme [35] and a regularization of the Bingham model. The numer-ical algorithm is implemented in C++ as an application of PELICANS2. We referthe reader to [44, 82] for a more detailed description of the numerical method and2 https://gforge.irsn.fr/gf/project/pelicans/; PELICANS is an object-oriented platform devel-oped at the French Institute for Radiological Protection and Nuclear Safety and is distributedunder the CeCILL license agreement (http://www.cecill.info/).22gxzLxLzLHρ1, µ1τYFluid 1c=1ρ2, µ2Fluid 2c=0u=w=0u=w=0u=0wx=0u=0wx=0Figure 2.1: A sketch of the geometry for the case of a rectangular initialblock.its implementation. We use the regularized scheme as the main computationaltool; the augmented-Lagragian algorithm is slower and was used more sparingly,specifically when looking at flow close to failure or during the final approach torest. In most situations, the agreement between the two computations is satisfac-tory (examples are given below in figure 2.4); only at the initiation or cessationof motion is there a noticeable difference, primarily in the stress field (discount-ing the solution for the plug, which is an artifact of the iteration algorithm in theaugmented-Lagrangian scheme).2.1.2 Model equationsWe quote conservation of mass, concentration and momentum for a two-dimensionalincompressible fluid in dimensionless form:∇ · u = 0, ∂c∂ t+(u ·∇)c = 0, (2.1)23ρRe[∂u∂ t+(u ·∇)u]=−∇p+∇ · τ −ρ zˆ, (2.2)In these equations, lengths x = (x,z) are scaled by the characteristic initial heightof the Bingham fluid,H , velocities u =(u,w) by the speed scaleU = ρ1gH2/µ1,and time t by H /U , where g is the gravitational acceleration; the stresses τand pressure p are scaled by ρ1gH . The Reynolds number is defined as Re =ρ1U H /µ1. Here, the subscript 1 or 2 on the (plastic) viscosity µ and density ρdistinguishes the two fluids, and linear interpolation with the concentration fieldc is used to reconstruct those quantities for the mixture; i.e. after scaling with thedenser fluid properties,ρ = c+(1− c)ρ2ρ1and µ = c+(1− c)µ2µ1. (2.3)In dimensionless form, the unregularized Bingham constitutive law isγ˙ jk = 0, τ < cB,τ jk =(µ +cBγ˙)γ˙ jk, τ > cB,τ =√12∑j,kτ2jk (2.4)whereB =τY Hµ1U≡ τYρ1gH(2.5)is a dimensionless parameter related to the yield stress τY , and the deformationrates are given byγ˙ jk =∂u j∂xk+∂uk∂x j, γ˙ =√12∑j,kγ˙2jk. (2.6)The regularized version that we employ isτ jk =(µ +cBγ˙ + ε)γ˙ jk, (2.7)24where ε is a small regularization parameter. We verified that the size of thisparameter had no discernible effect on the results presented below, as long as it issmall enough; we therefore consider irrelevant the precise form of the regulariza-tion (which is simple, but not necessarily optimal).We solve these equations over the domain 0 ≤ x ≤ ℓx = Lx/H and 0 ≤ z ≤ℓz = Lz/H , and subject to no-slip conditions, u = w = 0 on the top and bottomsurfaces (but see §2.2), and symmetry conditions on the left and right edges, u = 0and wx = 0. The computational domain is chosen sufficiently larger than the initialshape of Bingham fluid that the precise locations of the upper and right-handboundaries (i.e. ℓx and ℓz) exert little effect on the flow dynamics.Initially, both fluids are motionless, u(x,z,0) = w(x,z,0) = 0. We take the ini-tial shape of the viscoplastic fluid to be either a block or triangle; the concentrationfield then begins withc(x,z,0) = 1 for{0≤ x≤ X0 or X0(1− 12z),0≤ z≤ 1 or 2,and c(x,z,0) = 0 elsewhere, where X0 = L /H is the initial aspect ratio. Thedifferent maximum depths ensure that the initial conditions have the same areafor equal basal width X0.The main dimensionless parameters that we vary are the yield-stress parameterB (which we loosely refer to as a Bingham number) and initial aspect ratio X0.Unless otherwise stated, we set the other parameters to beℓx = 5, ℓz = 1.25,ρ2ρ1=µ2µ1= Re = 10−3.By fixing the density and viscosity ratios to be small, we attempted to minimizethe effect of the overlying viscous fluid (but see the discussion of the finger ofover-ridden fluid below). The relatively low Reynolds number reflects our inter-est in the limit of small inertia, although the PELICANS implementation requiresRe 6= 0, evolving the fluid from the motionless state; we established that adopt-25ing Re = 10−3 minimizes the effect of inertia beyond the initial transient. Someadditional technical details of the computations are summarized in A.1. In thisappendix, we also describe a second scheme that we used to study how the initialstate fails at t = 0; this scheme does not solve the initial-value problem, but calcu-lates the instantaneous velocity field at t = 0, assuming that Re = 0 and the initialstresses are in balance (the steady Stokes problem).2.2 Newtonian benchmarkThe collapse of the initial block of the heavier viscous fluid creates a slumpingcurrent that flows out above the bottom surface. However, because of the no-slip condition imposed there, the upper-layer fluid cannot be displaced from athin finger coating the base that is over-ridden by the advancing gravity current.Problematically, the finger becomes excessively thin and difficult to resolve withthe relatively small viscosity and density ratios that we used to minimize the effectof the lighter fluid. We illustrate the formation of the finger and its subsequentdevelopment in figure 2.2. A.1 features further discussion of the finger and itsevolution.The challenges associated with resolving the finger are illustrated in figure2.3, which plots the evolution in time of the flow front, X(t) (defined as the right-most position where c(X ,z, t) = 12), for computations with different grid spacing.The first panel in this figure shows the results using a relatively simple MUSCLscheme for tracking the interface [79], which was previously coded into PELI-CANS [44, 82]. This algorithm fails to track the interface well with the grid reso-lutions used: as the finger develops, it remains erroneously thick and the enhancedlubrication by the lighter viscous fluid causes the heavier current to advance tooquickly (cf. A.1.3).An interface-tracking scheme based on the PLIC algorithm proposed in [66]performs better; see figure 2.3(b). The flow front now advances less quickly.However, the fine scale of the finger still leads to a relatively slow convergenceof the computations with grid spacing ∆x = ∆z (corresponding to finite elements26xz0 0.5 1 1.500.20.40.60.811 1.2 1.400.010.020.030.04Figure 2.2: Snapshots of the evolving interface for a Newtonian heavier fluidat the times t = 1,2,3,4,5; the inset shows a magnification of the fingerof over-ridden lighter fluid.in the PELICANS code with equal aspect ratio). Moreover, the resolution failureis again manifest as an enhancement in the runout of the current that results froma finger that does not thin sufficiently quickly. In A.1.3 we argue that this is anintrinsic feature of the volume-of-fluid algorithm when the interface is containedwithin the lowest grid cell of the numerical scheme.The resolution problems with the finger are compounded in computations withBingham fluid, for which the yield stress further decreases the effective viscosityratio. Although some sort of local mesh refinement and adaptation would be anatural way to help counter these problems, we elected to avoid them in a differ-ent fashion which was more easily incorporated into PELICANS. In particular, bymonitoring c(x,z, t) at z = ∆z for each time step, we determined when the fingerwas expected to be contained within the lowest grid cell. At this stage, to pre-270 50 100 150 200 25011.522.533.54tX(a) MUSCL 0 50 100 150 200 25011.522.533.54tX(b) PLIC 0 50 100 150 200 25011.522.533.54tX(c) PLIC with correction dz=1/40dz=1/80dz=1/160dz=1/320slipasympt.dz=1/40dz=1/80dz=1/160dz=1/320asympt.xz (d)0 0.5 1 1.5 2 2.5 3 3.500.2 PLIC PLIC with correctionPLIC with slipFigure 2.3: Flow front X(t) plotted against time for computations with New-tonian fluids using (a) the simple MUSCL scheme, (b) the PLIC im-provement, and (c) the PLIC scheme with the lower boundary con-dition on c(x,z, t) adjusted according to the algorithm outlined in themain text. In each case, runs with different resolution are shown. For(c), the (red) dashed-dotted line labelled slip shows a solution com-puted with the PLIC scheme, but with no slip imposed on the heavierfluid and free slip imposed on the lighter fluid at z = 0. The circlesshow the prediction of the leading-order, shallow-layer asymptoticsin A.2. In (d), we plot the interface shape at t = 250 for the highestresolution solutions computed with the PLIC scheme with the threedifferent lower boundary conditions.vent the resolution failure from artificially restricting the thinning of the finger inthe volume-of-fluid scheme (see A.1.3), we reset the concentration field to c = 1at z = 0. The finger was thereby truncated and the effective contact line moved.Practically, we reset c(x,0, t) when c(x,∆z, t) exceeded 0.99 (the results were in-sensitive to the exact choice for this value). As shown by figure 2.3(c), this led tocomputations that converged much more quickly with grid spacing and fell closeto both the most highly resolved computations with the original PLIC scheme and28the predictions of shallow-layer theory. Nevertheless, the adjustment destroys theability of the code to preserve the volumes of the two fluids. For the computa-tions we report here, less than about one percent of the volume of the lighter fluidwas lost as a result of the adjustment. But, as the velocity profile was then fullyresolved near the boundary and no other unexpected problems were found, weconsidered this flaw to be acceptable. Hereon, all reported computations use thisadjusted boundary condition.To provide a physical basis for the adjustment scheme, we would need todemonstrate that it corresponds to the addition of another physical effect, suchas van der Waals interactions. We did not do this here, but simply note that theadjustment acts like the numerical devices implemented in contact line problemswith surface tension to alleviate the stress singularity and allow the contact line tomove [78]. Indeed, the scheme is much like limiting the dynamic contact angle tobe about 3pi/4 or less, by adjusting the interface over the scale of the bottom gridcell but without introducing explicitly any interfacial tension.An alternative strategy is to change the lower boundary condition so that thelighter fluid freely slips over the lower surface whilst the heavier fluid still satisfiesno slip. This strategy, which can be incorporated using a Navier-type slip law inwhich the slip length depends on c, leads to results that compare well with thescheme including the concentration correction (see figure 2.3(c)). However, forBingham fluid, the diffuse nature of the interface-tracking scheme and the PLICalgorithm eventually result in the recurrence of resolution problems over longertimes. By contrast, the adjustment scheme successfully survives the long timediffusion process.2.3 Bingham slumpsFor the collapse of a Bingham fluid, we first describe the general phenomenology,then explore the details of failure, and finally categorize the slumped end-states.Along the way, we indicate how the computations approach the asymptotic limitsof relatively shallow (low, squat) or slender (tall, thin) slumps.292.3.1 Slump and plug phenomenologyWhen the heavier fluid is viscoplastic, collapse only occurs provided the yieldstress does not exceed a critical value Bc that depends on initial geometry. ForB < Bc, the viscoplastic fluid collapses, but the yield stress eventually brings theflow to an almost complete halt (slumps with the regularized constitutive law nevertruly come to rest, and iteration errors in the augmented-Lagrangian scheme pre-vent the velocity field from vanishing identically). Figure 2.4 plots the position ofthe flow front X(t) and central depth H(t) = h(0, t) for computations with varyingB, beginning from a square (X0 = 1) initial block. For this shape, the critical valuebelow which collapse occurs is Bc ≈ 0.265, and unlike the inexorable advance ofthe Newtonian current (shown by a dashed line), X(t) and H(t) eventually con-verge to B−dependent constants (the case with the smallest B = 0.01 requires alonger computational time than is shown).Sample collapses from square blocks with B = 0.05 and 0.14 are illustratedin more detail in figure 2.5. The first example shows a slump with relatively lowyield stress, for which the fluid collapses into a shallow current. The case withhigher B collapses less far, with an obvious imprint left by the initial shape. Notethe stress concentration that arises for earlier times in the vicinity of the contactline (a feature of all the slumps, irrespective of initial condition and rheology).In general, for rectangular initial blocks with order one initial aspect ratio,we find that the flow features two different plug regions during the initial stagesof collapse (cf. [77]). First, at the centre of the fluid the stresses never becomesufficient to yield the fluid, and a rigid core persists throughout the evolution.Second, the top outer corner is not sufficiently stressed to move at the initiationof motion. This feature falls and rotates rigidly as fluid collapses underneath, buteventually liquifies and disappears for small B; at higher yield stress, the rigidcorner survives the fall and decorates the final deposit. As the fluid approachesits final shape and flow subsides, further plugs appear, particularly near the flowfront; for the deeper final deposits, these plugs appear to thicken and merge toleave relatively thin yielded zones.300 200 400 600 800 100011.522.533.544.5 B=0.1B=0.06B=0.03B=0.01X(a) X(t)B=0.20 200 400 600 800 10000.30.40.50.60.70.80.91tHB=0.1B=0.01(b) H(t)B=0.2Figure 2.4: (a) Front position X(t) and (b) central depth H(t) = h(0, t) forBingham dambreaks with a square initial block (X0= 1) and the valuesof B indicated. The dashed lines show the Newtonian results. For theviscoplastic cases, the circles show the result using the augmented-Lagrangian scheme and the lines indicate the result with a regularizedconstitutive law. The dotted lines show the prediction of the leading-order shallow-layer asymptotics for B = 0.01.When the initial shape is a triangle with X0 = O(1), only a few of the phe-nomenological details change (figure 2.6): the apex of the triangle now fallsrigidly as material spreads out underneath; this pinnacle decorates the final de-posit unless the yield stress is sufficiently small. Again, the slump features a rigid31Figure 2.5: Snapshots of a collapsing square with (a) B = 0.05 at t = 0, 2.5,5, 10, 20, 40, 70, 150, 450, 950, 4000, and (b) B = 0.14 at t = 0,10, 50, 100, 200, ..., 1000. The insets show density plots of the stressinvariant τ at the times indicated, with the yield surfaces drawn as solidlines (and common shading maps for the final two snapshots in eachcase).core and further plugs form near the nose over late times.With a relatively wide initial condition, the pattern of evolution is slightly32Figure 2.6: Snapshots of the evolving interface for a triangular initial condi-tion with X0 = 1 and (a) B = 0.05 and (b) B = 0.14, at t = 0, 2.5, 5, 10,20, 40, 100, 150, 450, 950, 4000. The insets show density plots of thestress invariant τ at the times indicated, with the yield surfaces drawnas solid lines (the two later time density plots have a common shadingscheme).33different: for the rectangle, the central rigid plug extends over the entire depthof the fluid layer and only the side of the block collapses. The final deposit thenfeatures a flat top at its centre, as illustrated in figure 2.7(a). Such incompleteslumps are predicted by shallow-layer theory to occur for 3BX0 > 1 [12]. Thisestimate can be improved to 3BX0 > 1− 3piB/4 using the higher-order theoryoutlined in A.2 (and specifically the final-shape formula (2.9)), which adequatelyreproduces numerical results for B < 0.15; at higher B, the computations indicatethat BX0 must exceed 0.25±0.015 for the slump to preserve a rigid central blockIncomplete slumps of a different kind arise for an initial triangle. In this case,collapses begin over the central regions where the initial stresses are largest, andmay not reach the edge, where fluid is initially unyielded, if the triangle is toowide. The (leading-order) shallow-layer model predicts that collapse occurs butdoes not reach the fluid edge at x = X0 if 4> BX0 > 9/8. As illustrated in figure2.7(b), such incomplete slumps are also observed numerically, though again for aslightly different range of initial widths (the plugged toe of the triangle is relativelysmall in the example shown).Slender (i.e. tall, thin) initial blocks, with X0 ≪ 1, also collapse somewhatdifferently, with fluid yielding only at the base of the fluid and remaining rigid inan overlying solid cap; see figure 2.8, which shows computations for rectangles(slender triangular slumps are much the same). The lower section of the columnsubsequently spreads outwards with the rigid cap descending above it in a mannerreminiscent of the (axisymmetrical) slump test [61, 68]. Interestingly, the slenderslump also generates undulations in the thickness of the column. As illustrated byfigure 2.8(a), these undulations (which do not occur in the Newtonian problem)appear towards the end of the collapse and are linked to zigzag patterns in thestress invariant and yield surfaces. The features follow characteristic curves ofthe stress field (the “sliplines”) that intersect the side free surface, and which haveslopes close to forty-five degrees (see A.4.1); the wavelength of the pattern istherefore approximately the width of the column.34Figure 2.7: Snapshots of the evolving interface for an initial (a) rectangularwith (X0,B) = (3,0.11), and (b) triangle with (X0,B) = (8,0.14), atthe times t = 0, 10, 20, 30, 50, 100, 400, 700, 1000. The dotted line in(a) shows a modification of the prediction in (2.9) that incorporates thecentral plug (and which terminates at finite height). The insets showdensity plots of the stress invariant τ at the times indicated, with theyield surfaces drawn as solid lines.2.3.2 Shallow flowAs illustrated above, when B ≪ 1 the fluid collapses into a shallow current with|∂h/∂x| ≪ 1 whatever the initial condition. Asymptotic theory [11, 52] then pro-vides a complementary approach to the problem. As illustrated in figure 2.4, the35Figure 2.8: Slumps of slender columns: the four images on the left showcollapsing columns for B = 0.3 and X0 = 0.025 at the times indicatedin the top right corner. Shown are the interface, yield surfaces andcolormap of the stress invariant τ . The evolving side profile for thiscollapse is shown in panel (a) at the times t = 0, 4, 9, 16, 25, 36, 49,100, 400, 1000. Panels (b)–(e) show columns at t = 50 for the valuesof B indicated, all with X0 = 0.025. Panel (f) compares the final sideprofiles (blue) with the predictions of slender asymptotics (red). Theshading scheme for τ for all the colormaps is shown in (e).sample collapse with B = 0.01 is relatively shallow and the numerical solutionsfor the runout and central depth match the predictions of the shallow-layer asymp-totics.The asymptotics also predict that flow becomes plug-like over a region under-neath the interface (see [6] and A.2). This “pseudo-plug” is not truly rigid but isweakly yielded and rides above a more strongly sheared lower layer. Horizontalvelocity profiles for a collapsing square block with B = 0.01 are shown in figure362.9(a) and compared with the predictions of the shallow-layer theory. Althoughthe pseudo-plugs are less obvious in the numerical profiles, the horizontal veloc-ity does become relatively uniform over the predicted pseudo-plug. Figure 2.9(b)illustrates how the superficial weak deformation rates associated with the pseudo-plug feature in a snapshot of log10 γ˙ , and how the transition to a more obviouslysheared layer underneath approximately follows asymptotic predictions.371 2 3 4 5 6x 10−300.050.10.150.20.25uzt=100200600 300(a) Profiles of uxh(c) Final profileB=0.01B=0.050 1 2 3 4 500.10.20.30.40.50.6Figure 2.9: Comparison of shallow-layer theory with numerical results for acollapsing square block (X0 = 1) with B = 0.01: (a) Horizontal veloc-ity profiles at x = 2.5 and the times indicated. (b) Logarithmic strain-rate invariant, A density map of log10 γ˙ on the (x,z)−plane, at t = 600.(c) Final shape. In (a) the crosses plot the numerical solution, and thepoints indicate the asymptotic profile (A.8); the star locates the bot-tom of the pseudo-plug. In (b), the solid (green) and dashed (blue)lines show the surface and fake yield surface (z = Y = h+B/hx) pre-dicted by asymptotics. In (c), the final profile for B = 0.05 is included.The solid lines show computed final states, the dotted lines denote theleading-order result (2.8) and the dashed line shows the higher-orderprediction (2.9). The dashed-dotted line is the asymptotic prediction(2.13) from [33].38Nevertheless, the numerical computations show notable disagreements withthe shallow-layer asymptotics. None of the true plugs appear in the asymptoticsolution, reflecting how they are associated with non-shallow flow dynamics: atthe midline, the asymptotics fail to incorporate properly the symmetry conditions,and at the flow front and the relic of the upper right corner, the surface alwaysremains too steep for a shallow approximation. The depth profiles predicted bythe asymptotics consequently do not show any of the finer secondary featuresimprinted by the true plugs, as illustrated by the profile for B = 0.05 also plottedin figure 2.9(c). Only when B is smaller are such features largely eliminated bythe fluid flow and the final shape well predicted by the asymptotics.Despite these details, the broad features of the numerical solutions are repro-duced by the shallow-layer asymptotics, particularly when first-order-correctionterms are included: for the final shape and if the initial block collapses completely,the leading-order solution ish(x) =√2B(X∞− x); (2.8)with the next-order corrections, we findh =√2B(X∞− x)+ pi2B (2.9)(A.2). The final runout X∞, or slump length, is dictated by matching the profile’sarea with that of the initial condition. This implies X∞ = (9X20/8B)1/3 for (2.8)and furnishes an algebraic problem to solve in the case of (2.9), with approximatesolution X∞ ≈ (9X20/8B)1/3[1− pi(B2/81X0)1/3]. As shown in figure 2.9(c), theprediction (2.9) agrees well with the shapes reached at the end of the numericalcomputations, even though the profile ends in a vertical cliff which violates theshallow-layer asymptotics.392.3.3 Slender slumpsFor a slender column with X0≪ 1, we can again make use of the small aspect ratioto construct an asymptotic solution. As summarized in A.3, this limit correspondsto theory of slender viscoplastic filaments [13] and indicates that the final state isgiven byx = ξ (z) =X02Bexp(− z2B)for 0≤ z≤ Z, (2.10)whereZ =−2Bh(0,0) log(2B) (2.11)is the height dividing yielded fluid from the overlying rigid plug. The fluid adoptsits original shape over Z < z < Z + 2Bh(0,0), having fallen a vertical distance(1−2B)h(0,0)−Z. It follows that the column will not slump if B > 12,X∞ =X02Band H∞ = 2Bh(0,0)[1− log(2B)]. (2.12)The latter is equivalent to the “dimensionless slump” reported previously [61, 68],although it was not declared as an asymptotic result relying on the column beingslender. The profile (2.10) is compared with the final profiles of numerical com-putations in figure 2.8. Aside from the relatively short-wavelength undulationsin column thickness over the yielded base of the fluid (whose lengthscale vio-lates the slender approximation), the asymptotics are in broad agreement with thenumerical results.Note that overly slender configurations are likely susceptible to a symmetry-breaking instability in which the column topples over to one side [8]. This is ruledout here in view of our imposition of symmetry conditions along the centrelinex = 0.402.3.4 FailureCritical yield stressThe critical yield stress, Bc, above which the fluid does not collapse is plottedagainst initial width, X0, in figure 2.10(a) for rectangular initial blocks. We cal-culate Bc in two ways: first, the final runout X∞ recorded in the slumps computedwith the PELICANS software (defined as in §2.3.5) converges linearly to the ini-tial width X0 as B → Bc. Second, in the inertialess problem, the initial stressesdictate the initial velocity field, and the maximum speed also falls linearly to zeroas B → Bc. Hence, we can determine Bc without performing any time steppingusing the scheme for Re = 0 described in A.1.2.410 0.2 0.4 0.6 0.8 1 1.20.250.30.350.40.450.5B c (a) Rectangles0 1 2 3 40.250.30.350.40.450.5X0B c(b) Triangles Figure 2.10: Critical yield stress, Bc, plotted against initial width X0 for (a)rectangular and (b) triangular blocks, found by from monitoring ei-ther the final runout X∞ (stars) or the initial maximum speed forRe = 0 (solid line). In (a), also shown are the bounds Bc = 0.2642and 0.2651 obtained from plastic limit-point analysis [55, 56, 62](dashed-dotted), results from slip-line theory [20] (dashed line withpoints) and the simple lower bound Bc =12(√X20 +1−X0) [21] (dot-ted line). The analogues of the latter two for the triangular blocks (seeA.4) are shown in (b). The dotted lines with open circles show im-provements of the simple lower bounds allowing for rotational failure(see A.4.2).42As illustrated in figure 2.10(a), Bc ≈ 0.2646 independently of X0 when X0 > 1.For such initial widths, collapses are incomplete and a solid core spans the fulldepth of the fluid, rendering the failure criterion independent of X0. The initialwidth matters for X0 < 1, leading to an increase of Bc. Eventually, Bc → 12 forX0→ 0, as expected from the slender column asymptotics in §2.3.3.Just below the critical yield stress, velocities are small and the Bingham prob-lem reduces to an analogous one in plasticity theory except over thin viscoplasticboundary layers. The incomplete slump is analogous to the classical geotechnicalproblem of the stability of a vertical embankment (e.g. [41]), provided no defor-mation occurs in z < 0. Classical arguments dating back to Coulomb, describe theform of failure in terms of the appearance of a slip surface dividing rigid blocks,allowing analytical estimations of Bc from balancing the plastic dissipation acrossthe slip surface with the potential energy release. In particular, assuming that thefailure occurs by the rotation of the top right corner above a circular arc, one ar-rives at an estimate Bc ≈ 0.261, after maximizing over all possible positions ofthe centre of rotation (cf. [24] and figure 2.11(a)). However, this type of solu-tion strictly provides only a lower bound on Bc and may not be the actual modeof failure. Indeed, the bound has been optimized and improved to 0.2642, and acomplementary upper bound computed to be 0.2651 [55, 56, 62]; the optimizationsuggests that failure occurs over a relatively wide region of plastic deformation[56]. The upper and lower bounds are included in figure 2.10(a), and are indistin-guishable on the scale of this picture, but bracket the value of Bc ≈ 0.2646 foundfor our Bingham slumps with X0 > 1.43−2 −1 0 100.511.52xz(a)−0.5 0 0.5 100.511.52zx(b)sα−0.5 0 0.5 1x(c)ssβFigure 2.11: Trial velocity fields to compute lower bounds on Bc for (a) thevertical embankment with a circular slip curve, and relatively slender(b) rectangular and (c) triangular initial blocks. In (b) and (c), thestraight (dashed) and circular (solid) failure surfaces for linear or ro-tational sliding motion are plotted; these surface can be parametrizedby the local slopes at the bottom corner, sα , and midline, sβ and s (forlinear sliding sα = sβ ).For a slender column, Chamberlain et al. [21] provide an estimate of thecritical yield stress by assuming that two lines of failure occur: the lowest cutsoff a triangular basal section, whereas the upper cleaves off a second triangle thatslides away sideways, leaving the remaining overlying trapezoid to fall vertically;44see figure 2.11(b). Optimizing the slopes of the two cuts furnishes the lowerbound, Bc =12(√1+X20 −X0), which is included in figure 2.10. Chamberlain etal. [20] also provide a numerical solution of the slipline field for a failure withthe form of an unconfined plastic deformation. This estimate converges towardstheir lower bound as X0 → 0, as indicated in figure 2.10(a). Our numericallydetermined values of Bc match Chamberlain et al.’s slipline solutions for X0 <0.5. For wider initial blocks, the slipline solutions deviate from the numericalresults and become inconsistent with the bounds for the vertical embankment forX0 > 0.8, highlighting how a different failure mechanism must operate.For triangles, no corresponding plasticity solutions exist in the literature. How-ever, Chamberlain et al.’s [20, 21] slipline solution and simple lower bound caneasily be generalized, as outlined in A.4 and illustrated in figure 2.11(c). Theslipline solution and bound are compared with numerical data in figure 2.10(b).Again, Bc → 12 for X0 → 0. Now, however, the slope of the initial free surfacecontinues to decline as X0 is increased, and so there is no convergence to a limitthat is independent of width. Instead, the shallow-layer theory of A.2 becomesrelevant and predicts that Bc → 4/X0 for X0 ≫ 1 (a limit lying well beyond thenumerical data in figure 2.10(b)).Note that Chamberlain et al.’s lower bounds can be improved by allowing thetriangle at the side to rotate out of position rather than slide linearly; see A.4.2.The failure surfaces then become circular arcs, as illustrated in figure 2.11(b,c).For the rectangle, the resulting improvement in the bound on Bc amounts to afew percent and is barely noticeable in figure 2.10(a). More significant is theimprovement of the bound for the triangle, which is brought much closer to thenumerical and plasticity results; see figure 2.10(b).Flow at failureThe failure modes of our rectangular viscoplastic solutions (for t = Re = 0) areillustrated in figure 2.12. The thinner initial columns yield only over a triangularshaped region at the base which closely matches that predicted by slipline theory45Figure 2.12: Strain-rate invariant plotted logarithmically as a density on the(x,z)−plane for solutions with B = 0.99BC and t = Re = 0 (A.1.2),at the values of X0 = 0.2,0.5,0.6,0.9,1,1.4 indicated by the x−axis.Also shown are a selection of streamlines. In (a)–(c), the solid bluelines indicate the border of the plastic region and expansion fan of thecorresponding slipline solutions (A.4). In (e) and (f), the blue linesindicate the circular failure surface of the lower bound solution.[20] (see also A.4.1). The failure mode changes abruptly when X0 slightly exceedsabout 0.5. The failure zone then takes the form of a widening wedge extendingfrom the lower right corner of the initial block up to the centre of the top, with theentire side face remaining rigid. Evidently, this mode of failure is preferred overthe slipline solution at these values of X0, leading to the departure of the observedvalues of Bc from Chamberlain et al.’s curve in figure 2.10(a). The failure modechanges a second time for X0≈ 1: wider initial blocks fail chiefly over a relatively46Figure 2.13: A series of pictures similar to figure 2.12, but for initial trian-gles (with solid blue lines showing the yield surfaces of the sliplinesolution in (a), and the circular arcs of the bound of A.4.2 for rota-tional failure in (b) and (c)).narrow layer connecting the lower left corner to an off-centre location on the topsurface, which lies close to the circular failure surface of the simple lower boundin figure 2.11(a).For both the narrower and wider examples in figure 2.12, the failing deforma-tions are dominated by sharp viscoplastic boundary layers that likely correspondto yield surfaces. Spatially extensive regions (in comparison to the fluid depth orwidth) of plastic deformation do, nevertheless, occur, and the failure modes nevertake the form of a patchwork of sliding rigid blocks. Note that it is difficult tocleanly extract the yield surfaces as B→ Bc, which complicates the identificationof the failure mode. In figure 2.12, we have avoided showing these surfaces andinstead displayed the deformation rate and a selection of streamlines. Curiously,for X0 > 0.5, it is difficult to envision how one might construct corresponding sli-pline fields (there are no surfaces bounding the plastic zone with known stressesthat can be used to begin the slipline construction).47For a triangular initial condition, failure for smaller widths again occurs throughthe appearance of a lower plastic zone that compares well with slipline theory;see figure 2.13. This agreement is confirmed by the match of the observed criti-cal yield stress, Bc(X0), with the slipline predictions in figure 2.10(b). Unlike therectangle, however, there is no abrupt change in failure mode as X0 is increased,at least until the surface slope of the triangle becomes too shallow to apply Cham-berlain et al.’s construction [20] (see A.4.1). At the largest widths, the trianglefails at the centre but not the edge, as already noted in section 2.3.1.2.3.5 The final shape and slump statisticsTo extract statistics of the final shape, we define a stopping criterion accordingto when the stress invariant first becomes equal or less than B at each point inthe domain. The resulting “final state” appears to be reached in a finite time(for both augmented-Lagrangian and regularized computations), in disagreementwith asymptotic theory [58], which predicts that flow halts only after an infinitetime (see also A.3). A selection of final profiles for varying Bingham number isdisplayed in figure 2.14 for both square and triangular initial conditions.Plasticity theory is again relevant in the limit that the slump approaches itsfinal state. This fact was used previously [33] to construct the final profiles withslipline theory, following earlier work by Nye [59]. A key assumption of thisconstruction is that the flow is under horizontal compression. The assumption canalso be used to continue the shallow-layer asymptotics to higher order to predictthe final profile,h =√2B(X∞− x)+ pi24B2− pi2B, (2.13)which agrees well with the slipline theory [33]. Unfortunately, neither the sliplineprofiles or (2.13) compare well with laboratory experiments.48xzB ↑(a) Square0 0.5 1 1.5 2 2.5 3 3.5 4 4.500.20.40.60.81xz(b) TriangleB ↑0 0.5 1 1.5 2 2.5 3 3.5 4 4.500.20.40.60.811.21.41.61.82Figure 2.14: Profiles of the final deposit, starting from (a) a square blockand (b) a triangle, with X0 = 1, for B = 0.01, 0.02, 0.04, ..., 0.22 and0.24, together with the initial states.Figure 2.15: Final numerical solution for the slump of an initial square withB = 0.02, showing density maps of (a) pressure, (b) τxx and (c)τxz, and (d) the slipline field diagnosed from the numerical solution(solid) and built explicitly by integrating the slipline equations start-ing from the curve (2.9) (dotted). In (d), the plugs are shaded black,and no attempt has been made to match up the two sets of sliplines.49A sample final state from the numerical computations with a diagnosis of theassociated slipline field is displayed in figure 2.15. For the latter, we map outcurves of constant p− z±2Bθ , which are the invariants that are conserved alongthe two families of sliplines, where θ = −12tan−1(τxx/τxz) [64] (see also A.4.1).As also shown by figure 2.15(d), the resulting curves compare well with an explicitcomputation of sliplines launched from the surface position predicted by (2.9),where p = 0 and the sliplines make an angle pi/4 with the local surface tangent.The sliplines in figure 2.15 follow a different pattern to Nye’s construction(see figures 5 and 6 in [33]). The reason for this discrepancy is that almost all ofour slumps comes to rest in a state of horizontal expansion, rather than compres-sion (we have observed regions under horizontal compression only in the slumpsof relatively wide triangles, as in figure 2.7(b)). Correcting this feature of the dy-namics leads to the revised higher-order asymptotics summarized in A.2, whichculminates in the prediction for the final profile in (2.9). As is evident from figure2.9(c), the asymptotic predictions for horizontal expansion compare much morefavourably with the numerical results than the slipline theory and asymptotics forhorizontal compression.The comparison is illustrated further in figure 2.16, which shows scaled finalrunouts and depths, X∞/√X0 and H∞/√X0, as functions of B/√X0 for the numer-ical computations, slipline theory and the various versions of the shallow-layerasymptotics. Scaling the runout and depth in this fashion corresponds to choosingthe square root of the initial area as the lengthscale in the non-dimensionalizationof the problem. The slipline theory and various versions of the shallow-layerasymptotics furnish curves of X∞/√X0 and H∞/√X0 against B/√X0 that are in-dependent of X0, implying an insensitivity to the initial condition. By contrast,the deeper final profiles of the numerical solutions with higher B, and their scaledfinal runout and depth, do depend on X0 and the initial shape. This dependence ishighlighted in figure 2.17, which compares data for initial triangles and rectangles.500.05 0.1 0.15 0.2 0.2511.522.533.544.55B /√X0X∞/√X0 (a)kaolin Carbopol 0.30 wt.%Carbopol 0.15 wt.%Joint Compound Cochard & Ancey0.05 0.1 0.15 0.2 0.250.20.40.60.811.2B /√X0H∞/√X0 (b)NumericsStaron et alSliplineH∞<<1H∞<<1, expansionH∞<<1, compressionX∞<<1 (Pashias et al)Figure 2.16: Scaled final (a) runout X∞/√X0 and (b) central depth H∞/√X0as a function of B/√X0 for Bingham slumps from square initial con-ditions (solid lines with dots). Also plotted using the symbols indi-cated are experimental data from [33] and [26] for slumps of aqueoussolutions of Carbopol, kaolin and “Joint Compound”. The leadingorder asymptotic prediction (2.8) is shown by the dotted line; thesolid lines with circles or squares plot the predictions in (2.9) and(2.13); the dashed lines shows the results of slipline theory [33]. Theslender-column asymptotic prediction in (2.12) with X0 = 1 is shownby the dotted line and pentagrams. The dotted line and hexagramsshow the fit proposed by Staron et al. [77].510 0.1 0.2 0.3 0.411.522.533.544.5B/√X0X∞/√X0(a) 0 0.1 0.2 0.3 0.40.511.522.5B/√X0H∞/√X0(b)X0=1.5, blockX0=0.5, block X0=1, blockX0=1.5, triangleX0=0.5, triangleX0=1, triangleFigure 2.17: Scaled final (a) runout X∞/√X0 and (b) central depth H∞/√X0as a function of B/√X0 for Bingham slumps from rectangular andtriangular initial conditions with X0 = 0.5, 1 and 1.5. The solid linesshow the predictions of the higher-order asymptotics from (2.9) for acomplete slump.Figure 2.16 also includes data from laboratory experiments with Carbopol[26, 33] and some other fluids, which were conducted by Dubash et al. though notreported in their paper. None of these fluids are well fitted by the Bingham model,with a Herschel-Bulkley fit being superior. However, the final state is controlledby the yield stress and likely independent of the nonlinear viscosity of the mate-rial (at least provided inertia is not important), permitting a comparison betweenthe experiments and our Bingham computations. Although the numerical resultscompare more favourably with the experiments than the slipline theory, the com-parison with the leading-order asymptotic prediction is just as good. Thus, thediscrepancy between theory and experiment noted by [33] is only partly due tothe assumption that the slump came to rest in a state of horizontal compression,but other factors must also be at work, such as stresses in the cross-stream di-rection and non-ideal material behaviour. Note that the range of dimensionless52yield stresses spanned by the experimental data lie well into the regime wherethere should be no significant dependence on the initial shape. This is comfort-ing in view of the fact that the experiments were conducted using different initialconditions, either by raising a vertical gate or tilting an inclined tank back to thehorizontal (which correspond roughly to our rectangular or triangular initial con-ditions).Finally, figure 2.16 includes the predictions of the slender column asymptoticsin (2.12) (see A.3) for X0 = 1, and a formula presented by Staron et al. [77] basedon their volume-of-fluid computations with GERRIS and a regularized constitu-tive law. Given that the slumps from which the data in figure 2.16 are taken are notslender, it is not surprising that (2.12) compares poorly with the numerical results.We suspect that the disagreement between our results and the fit of Staron et al.[77] originates either from an inadequate resolution of the over-ridden finger ofless dense fluid or the significance of inertial effects in their computations. Indeed,these authors quote a final shape that depends explicitly on the plastic viscosityof the heavier fluid, whereas this physical quantity is completely scaled out in ourcomputations when Re→ 0.2.4 Concluding remarksA yield stress introduces two key features in the collapse and spreading of a vis-coplastic fluid: failure occurs only provided the yield stress can be exceeded, and,when flow is initiated, the yield stress eventually brings motion to a halt. Here wehave provided numerical computations of the two-dimensional collapse of Bing-ham fluid, exploring the phenomenology of the flow for different initial shapes.We compared the results with asymptotic theory valid for relatively shallow (low,squat) or slender (tall, thin) slumps, and with solutions from plasticity theory ap-plying near the initiation and termination of flow. We verified that the compu-tations converge to the asymptotic solutions in the relevant limits and identifiedwhere the plasticity solutions apply. We studied both the initial form of failure,extracting criteria for a collapse to occur, and the shape of the final deposit, com-53paring its runout and depth with previous experiments and predictions.There are three key limitations of our study with regard to the collapses of vis-coplastic fluids in engineering or geophysical settings. First, our two-dimensionalgeometry is restrictive and an axisymmetric assumption prefereable for a range ofapplications such as the slump test. Such a generalization raises the interestingquestion of how incompressible viscoplastic flow avoids the inconsistency of thevon-Karman-Haar hypothesis [64] (This will be further explained in Chapter 3).Second, the issues associated with the no-slip condition on the underlying sur-face are not merely numerical: viscoplastic fluids can suffer apparent slip [48],demanding the inclusion of a slip law. Finally, inertia is important in many ap-plications, an effect that allows slumps to flow beyond the rest states we havecomputed. Other interesting generalizations include the incorporation of differentrheologies, such as thixotropy, and surface tension at small spatial scales.54Chapter 3Axisymmetric viscoplasticdambreaks and the slump test1One goal of this chapter is to provide a reliable solution of the benchmark problemproposed in [69], and to offer a more complete description of the slump behaviourover a wider range of physical conditions. For the task, we perform computationsbased on the VOF method to deal with the fluid interface and exploiting speciallydesigned codes to capture the yield-stress rheology. Our study follows on from thechapter 2[53] in which we considered two-dimensional dambreaks of viscoplas-tic fluid. We complement the computations with asymptotic theory for shallowgravity currents and slender vertical columns, and bounds from plasticity theoryto constrain the mode of failure for slumps near the critical yield stress whereatno collapse actually occurs.A second goal is to compare our theoretical modelling with experiments, col-lating some of the existing measurements from the literature [61, 70]. These ex-periments have not previously been performed sufficiently thoroughly to disentan-gle the effects of material rheology, the mechanism of release, and any interaction1A version of chapter 3 has been published. [Liu, Y.], Balmforth, N.J. & Hormozi, S. (2018)Axisymmetric viscoplastic dambreaks and the slump test. J. Non-Newtonian Fluid Mech. 258,45-57. [54].55with the underlying surface (effective slip). Therefore, we also perform our ownsuite of experiments using aqueous suspensions of Carbopol. This suspension isoften suggested to be well characterized by a Herschel-Bulkley rheology and po-tentially eliminates some of the confounding effects brought into experiments bynon-ideal material behaviour [14]. We thereby provide a demanding test of thetheory whilst gauging the effects of the release mechanism and any effective slip.3.1 Formulation3.1.1 Problem set-up and solution strategyFigure 3.1: Sketch of the geometry.The geometry of the problem is shown in figure 3.1: we use an axisymmetriccylindrical polar coordinate system (r,z) to describe the sudden release of a cylin-der of incompressible Bingham fluid with radius Rˆ and height Hˆ. The fluid hasdensity ρ1, yield stress τY and plastic viscosity µ1 and is immersed in an ambientNewtonian fluid with density ρ2 and viscosity µ2. The density and viscosity ratiosare set at the small valuesρ2ρ1= µ2µ1 = 0.002 in order to minimize the effects ofthe ambient fluid (we have verified that the precise values of these ratios have nosignificant effect on the computations once one deals with the resolution issuesdescribed in §2.3). We use the VOF method to track the fluid interface, which56introduces an advected concentration field c(r,z, t) to distinguish the fluid phase:c = 1 represents the viscoplastic fluid, and c = 0 denotes the Newtonian ambient.Given c, bulk material parameters are computed using linear interpolation. Theinitial configuration isc(r,z,0) ={1 for 0≤ r ≤ Rˆ,0≤ z≤ Hˆ,0 elsewhere.The two fluids are miscible, eliminating interfacial tension.The computation domain extends to a height Lz and radius Lr, which are cho-sen to ensure that these boundaries remain remote and do not affect the evolutionof the slump. We impose no slip (u = 0,w = 0) on the bottom z = 0, regularityconditions on the axis r = 0 (i.e. u = ∂w/∂ r = 0), and no normal flow and freeslip conditions along r = Lr and z = Lz.We use two methods to deal with the yield-stress constitutive law: a regu-larization scheme that treats the unyielded region as a highly viscous fluid, andan augmented-Lagrangian scheme that explicitly treats the yield stress within aweak formulation of the problem [35, 80]. Both are implemented in C++ as anapplication of the PELICANS platform (e.g. [44]).3.1.2 Dimensionless model equationsWe scale lengths by the initial height Hˆ, velocities by the speed scaleU = ρ1gHˆ2/µ1,and time by Hˆ/U , where g is the gravitational acceleration; the stresses and pres-sure are scaled by ρ1gHˆ. The governing equations for the concentration field c,velocity u = (u,w), deviatoric stress tensor τ , and pressure p are then∇ · u = 0, ∂c∂ t+(u ·∇)c = 0,ρRe[∂u∂ t+(u ·∇)u]=−∇p+∇ · τ −ρ(01),(3.1)57whereρ = c+(1− c)ρ2ρ1and µ = c+(1− c)µ2µ1. (3.2)The unregularized Bingham constitutive law, used in the augmented-Lagrangianmethod, is γ˙ jk = 0, τ < cB,τ jk =(µ +cBγ˙)γ˙ jk, τ > cB,(3.3)whereas the regularization method uses the variant,τ jk =(µ +cBγ˙ + ε)γ˙ jk, (3.4)where τ =√12 ∑ j,k τ2jk and γ˙ =√12 ∑ j,k γ˙2jk denote second tensorial invariants, andthe deformation rates are given byγ˙rr γ˙rθ γ˙rzγ˙θ r γ˙θθ γ˙θ zγ˙zr γ˙zθ γ˙zz=2ur 0 uz +wr0 2u/r 0uz +wr 0 2wz , (3.5)with subscripts represent partial derivatives, except in the case of tensor compo-nents. The scalings introduce the dimensionless initial radius (or aspect ratio), andthe Reynolds and Bingham numbers,R =RˆHˆ, Re=ρ1UHˆµ1and B =τY Hˆµ1U. (3.6)The regularization parameter ε in (3.4) is taken to be 10−8, which was verified tobe sufficiently small that the modification of the constitutive law had an insignifi-cant effect on the results reported below (but see the comment at the end of §2.3).Given the concentration field, we define the instantaneous position of the surfaceof the slump to be given by c(r,z = h) = 12. For most of the computations, weselect parameters so that inertial effect is small, Re= 10−3, and vary B and R.583.1.3 PLIC scheme with interface correctionThe piecewise-linear-interface-construction (PLIC) [40, 66] is a contemporarystandard in the VOF method. The interface is represented by a line segment ineach grid cell, which is computed using the volume fraction c, as in [40]. Then,given the velocity field, the line segment is advected to a new position, and c up-dated accordingly. In view of the boundary conditions, there is no flux of c into orout of the domain, which is incorporated into the scheme in the manner in whichthe solver advects c along the boundary. Importantly, the bottom boundary is noslip, which does not permit the contact line to move and introduces an awkwardresolution issue, as in the 2D problem in chapter 2[53]. More specifically, as theviscoplastic fluid collapses and spreads out, a finger of ambient fluid adhering tothe bottom surface is over-ridden. For the relatively low density and viscosity ra-tios that we employ, this finger lubricates the slumping viscoplastic current andthins dramatically to introduce the resolution issue [53]. A key problem is thatthe scheme fails to accurately evolve c when the interface is inside the lowest gridcell, leaving the finger artificially thick and lubricating.A common way of moving the contact line in problems with surface tensionand Newtonian fluid is to replace the boundary condition with another that permitsslip. However, numerical solutions may not converge with mesh refinement [1].Instead, in chapter 2[53] we suggested a correction scheme that eliminates thefinger in a different way, allowing the computations to remain well resolved overlong times. The main point is that, with no slip, a finger of ambient fluid must stillcoat the underlying surface. However, counter to the un-corrected VOF scheme,the finger actually becomes too thin to lubricate the slump and should instead beignored. Practically, the scheme implements this idea by removing all the ambientfluid from a grid cell adjacent to the base when c exceeds a threshold near unity(chosen to be 0.99). This procedure clips the interface when it invades the lowestgrid cells, thereby truncating the finger and rendering the computation convergentin grid spacing [53].Despite the success of the scheme for 2D dambreaks, the algorithm is not con-59servative, with the mass of viscoplastic fluid growing with time. This awkwardfeature does not impair computations in 2D, but it does become more problematicin axisymmetric geometry, for which the convergence of the solutions with meshrefinement is weakened. For the current computations we therefore modified thecorrection scheme so that it conserved mass. In particular, whenever a correctionto c was implemented, and some of the ambient fluid removed from the one of thelowest grid cells, the lost material was added back by uniformly redistributing itinto the grid cells containing the interface (where c = 12). Essentially, this redistri-bution incurs an error in the position of the interface that is of the order of a smallfraction of the grid spacing, but in such a way that mass is conserved. Thoughhard to justify from a physical perspective, the resulting conservative correctionscheme converges more satisfyingly with mesh refinement for slumps with New-tonian fluid than the original non-conservative scheme (see B.1), chiefly becausethe latter suffers a resolution-dependent mass loss.The conservative correction scheme also performs better for Bingham fluidusing either the regularized constitutive law or the augmented-Lagrangian solver.Moreover, both produce comparable results for the global properties of the slumps(see, for example, figure 3.2 below), and their interaction with the PLIC schemedoes not introduce any additional unexpected issues. In more detail, the regular-ization scheme has the drawback that the fluid can never truly come to rest andcan fail to correctly predict the positions of the yield surfaces [38]. However, theregularization scheme is faster than the augmented-Lagrangian code. Therefore,we use regularization scheme for exploring global features of the slumps (definingthe flow to have come to rest when Max(|v|) < 10−5), while for the finer detailssuch as the yield surfaces, we compare both schemes.3.2 Bingham slumpsFigure 3.2 shows the slumps of cylinders of Bingham fluid for R = 1 and varyingyield stress B. Displayed is the central depth H (t) along with the radial positionof the flow front R(t) (defined as the largest radial extent of the interface, c = 12).600 100 200 300 400 500 600 700 800 900 10000.30.40.50.60.70.80.91(a)tH(t)0 200 400 600 800 100011.522.5tR(t)(b)Figure 3.2: Time series of flow height H (t) and front position R(t) forBingham slumps with R = 1 and B = 0.01, 0.03, 0.05, 0.08, 0.125, 0.2and 0.3. Blue curves are from regularization, and red circles from theaugmented-Lagrangian method.The computations suggest that flow is arrested in finite time and illustrate howcollapse only occurs when the yield stress is below a critical value Bc. A selectionof the final profiles is illustrated in figure 3.3(a).61(a)B ↑0 0.5 1 1.5 2 2.500.20.40.60.810 0.5 1 1.5 2 2.5 3 3.5 4 4.500.51 r(b)B=0.19B=0.1Figure 3.3: Profiles of the final deposit of cylindrical slumps with (a) R = 1,B= 0.01, 0.02, 0.04, 0.07, 0.125, 0.2 and 0.275, and (b) R= 4, B= 0.1and 0.19. The solid curves show numerical computations, the dashedcurves show the improved shallow-layer asymptotic result in (3.11)(for the lowest four values of B in (a)), and the dotted lines show theinitial cylinders.Further details of the phenomenology of a slump are shown in figure 3.4. Ini-tially, the stresses exceed the yield value throughout the cylinder except over asmall conical region at its core. Fluid subsequently slumps outwards, reducingthe stresses and allowing that plug to expand with time. Eventually, stresses de-cline towards B all the way to the flow front, bringing the fluid to rest. Notethat, there is only a single central plugged region; for equivalent two-dimensionaldambreaks [53] plugs persist at the periphery of the initial configuration, leadingto sharp corners that decorate the final deposit. Rigid features of this sort cannotoccur in axisymmetric slumps because the fluid edge must expand in order to fall.However, analogous weakly yielded zones persist during the collapse, then plug62up to create sharp rings that disfigure some of the final shapes (see figure 3.3).Figure 3.4: Evolution of the interface for a slumping cylinder with (R,B) =(1,0.1). The main panel shows the interface at the times t = 0, 2, 4,6, 8, 20, 200 and 1000. The three insets show the stress invariant τ asa density on the (r,z)−plane for the times t = 2, 20 and 1000 (with acommon scale for the last two cases). The green curve shows the yieldsurface.As for 2D dambreaks [53], when the initial configuration is sufficiently wide,fluid only collapses near the edge, leaving an unyielded flat-topped central section.Such “incomplete slumps” arise when the central plug spans the entire fluid layerand are illustrated in figure 3.3(b).3.2.1 Shallow flowWhen flow is shallow and inertia is negligible, lubrication theory can be appliedto obtain analytical results [11]. With our current scaling of the problem, thislimit is achieved when B ≪ 1. To account for the low aspect ratio, we rescalethe horizontal coordinate r = ε−1χ , deviatoric stress components τi j = ετˇi j and63Bingham number B = εBˇ using a small parameter ε ≪ 1, and then restate theforce balance equations for the viscoplastic layer:−pχ + εχ(χτˇrr)χ − εχτˇθθ + τˇrz,z = 0,−pz + ε2χ(χτˇrz)χ + ετˇzz,z = 1.(3.7)Neglecting the ambient fluid, the surface of the current can be located by theelevation z = h(χ , t), and is force free, demanding thatτˇrz+hχ(p− ετˇrr) = 0p− ετˇzz + ε2hχ τˇrz = 0}at z = h. (3.8)At leading order, we find that p ∼ h− z and τˇrz ∼ −hχ(h− z). The constitutivelaw and depth-integrated continuity equation can then be used to derive an evo-lution equation for h(χ , t) [11, 51]. However, the final profile arises when radialspreading speeds subside and τˇrz → Bˇ at the base of the fluid layer z = 0. Thenceτˇrz(χ ,0)∼−hhχ ∼ Bˇ, which givesh(r)∼√2B(r∞− r), (3.9)in terms of the original variables, where r∞ is the final radius.In chapter 2[53] a higher-order approximation for the final profile of a 2Dslump was developed by continuing the asymptotic solution to O(ε), assumingthat τˇ → Bˇ throughout the fluid layer. We follow suit here, although a signifi-cant complication emerges owing to the axisymmetric geometry. In particular, inaddition to the two force balance equations, we must also satisfyτˇrr + τˇθθ + τˇzz = 012(τˇ2rr + τˇ2θθ + τˇ2zz)+ τˇ2rz = Bˇ2,(3.10)leaving us one equation short for determining the full stress state (i.e. we have four64equations for the five unknowns p, τˇrr, τˇθθ , τˇrz and τˇzz). The origin of this indeter-minacy is the component τˇθθ , which does not arise in 2D and highlights how thestress field cannot, in general, be constructed independently of the velocity field.The situation is identical to classical plasticity theory where the so-called vonKarman-Haar hypothesis is often invoked to avoid this problem. The hypothesis,which states that τˇθθ must equal one of the principal stresses in the (r,z)−plane,implies that τˇ2θθ ≡ 13 Bˇ2 [42]. However, τˇ2rz → Bˇ2 at z = 0 for the leading orderlubrication solution, indicating that τˇθθ must vanish at the base of the fluid layer.Therefore, the von Karman-Haar hypothesis contradicts the leading-order asymp-totic solution and cannot be invoked here.Instead, we add the approximation τˇrr ∼ τˇθθ , which is suggested by both thevelocity field of the leading-order asymptotic solution and the numerical compu-tations; see B.2. With this alternative hypothesis, the asymptotic analysis can becontinued to O(ε) in order to arrive at the higher-order asymptotic approximation,h(r)∼√2B(r∞− r)+√34piB (3.11)(again in terms of the original variables).The predictions in (3.9) and (3.11) are compared to a numerical simulation fora slump with B = 0.0074 and R = 0.2546 in figure 3.5. Figure 3.3 also comparesthe improved approximation (3.11) with computed final shapes over a wider rangeof B. Note that (3.11) predicts that the final profile ends in a vertical cliff, violatingthe shallow-layer asymptotics. Nevertheless, (3.11) provides a meaningful pre-diction along the flow body that furnishes a better approximation than the leadingorder result (3.9), even when the flow is not particularly shallow (B > 0.04).In the example of figure 3.5, the fluid yields significantly almost everywhere,removing any sign of the initial shape in the final profile. Indeed, computationsthat begin with the same amount of fluid but conical initial shape also lead tosimilar final profiles. Figure 3.5 includes a computation using an initial cone witha top radius of Rtop = 1/6 and bottom radius of Rbase = 1/3, which corresponds65rz−1 −0.500.050.1rˆ (m)zˆ(m)0 0.1 0.2 0.300.010.020.030.04Figure 3.5: Final profile for a cylindrical dambreak with B = 0.0074 andR = 0.2546. The thin black curve shows the numerical simulation.On the left, the data is plotted using dimensionless variables, and theleading-order and improved asymptotic predictions in (3.9) and(3.11)are included as the dotted (green) line and (red) points, respectively.On the right, the data is replotted in dimensional variables assumingHˆ = 0.3m. Also shown are three other final profiles: one from a sim-ulation beginning with a cone with the same volume, with a top radiusof Rtop = 1/6 and bottom radius of Rbase = 1/3 (the ASTM standardgeometry; thick light grey line); a second from a simulation with theregularized Bingham model in which the regularization parameter isincreased to 10−4 (dashed line); and a third from a simulation thatdoes not apply the interface correction scheme to remove the underly-ing finger of ambient fluid (thin red contour).to the ASTM standard geometry; the final state cannot be distinguished from thatof the cylindrical dambreak in the plot. Thus, the example shown in this figurecorresponds to the benchmark problem proposed in [69]. Indeed, on the right-hand side of the figure, the profiles are replotted in dimensional variables usingHˆ = 30cm; the results can then be directly compared with figure 5 in [69].The agreement amongst the computations in [69] is relatively poor, a discrep-ancy that may arise from different treatments of the contact line and the fluidrheology. By contrast, the current computations align satisfyingly with (3.11),have converged with respect to mesh refinement, and are independent of the nu-merical algorithm (augmented-Lagrangian or regularization). Our computationindicates that the final (dimensional) radius is 28cm, in comparison to the averageof 29.5cm quoted in [69]. Figure 3.5 also tells the cautionary tale of the resultswhen we fail to resolve the finger of upper-layer fluid by applying the uncorrectedPLIC scheme, or overly regularize the constitutive model. Both deficiencies lead66to an enhancement in spreading. Note that there is essentially no effect of iner-tia on the profiles shown in figure 3.5: recomputing the results with a Reynoldsnumber (O(1)) based on the benchmark conditions, rather than the artificially lowvalue used for the bulk of our simulations, leads to no significant differences.When a complete slump does not occur but a central section of the initialcylinder survives, the shallow-layer solution is modified accordingly:h =1, 0< r < r∞− 12B√2B(r∞− r)+ pi√34B, r∞− 12B< r < r∞(3.12)This approximation is again compared with numerical final shapes in figure 3.3(b).3.2.2 Slender columnsFor a tall slender column (R ≪ 1), the asymptotic analysis of [53] can be gen-eralized to determine the instantaneous radius r = r(a, t)≪ 1 in terms of a La-grangian coordinate a ∈ [0,1] corresponding to initial height: where the fluid islocally yielded, we findr2(a, t) = E(t)r20(a)+1−E(t)√3B∫ 1ar20(x)dxz(a, t) =∫ a0r20(x)r2(x, t)dx, E(t) = e− Bt√3 ,(3.13)where r0(a) represents the shape of the initial column. As t → ∞, we then obtainr2(a) =1√3B∫ 1ar20(x)dx,z(a) =√3B ln[∫ 10 r20(x)dx∫ 1a r20(x)dx].(3.14)67The solution in (3.13) and (3.14) must be matched to a plugged upper sectionof the column, which always remains unyielded. The yield surface is given byr0(aY ) = r(aY ). The dimensionless slump height s (the difference between theinitial and final heights) is therefores = aY +√3B ln[√3Br20(aY )∫ 10 r20(a)da]. (3.15)For an initial cylinder, r0(a) = R, ands = 1−√3B+√3B ln(√3B), (3.16)which is widely used as a prediction of the slump test [61]. This immediatelyimplies that the column will not yield anywhere if B > Bc = 1/√3≈ 0.577.Figure 3.6 compares numerical computations with the slender-column asymp-totics for initial shape with r0(a) = R = 0.025. First, the dynamical evolution of acolumn with B= 0.15 is shown; second, the final shape is compared for cases withvarying yield stress B. Note that undulations in the surface profile do not appearnear the base of the column in the axisymmetric simulations, unlike in 2D [53].As a result, the computations and slender-column theory agree more satisfyingly,except at the very base of the column where the no slip condition is not correctlycaptured by the asymptotics.68rz(a) B=0.15, t=1, 2, 4, 8, 10000 0.02 0.0400.10.20.30.40.50.60.70.80.910 0.02 0.0400.10.20.30.40.50.60.70.80.91(b) Final shapes, varying BrFigure 3.6: (a) Collapse of a cylindrical column with B = 0.15, showingthe interface at the times indicated. (b) Final shapes of cylinder withR = 0.025 and B = 0.15, 0.20, 0.25, 0.30, 0.35, 0.40 and 0.45. Dashedlines indicate the slender asymptotic predictions (3.13) and (3.14);solid lines are from numerical simulations.3.2.3 Failure modeJust below the critical value Bc(R) for which no slump will occur, the collapse ischaracterized by a particular mode of failure that depends on the initial shape.To compute the critical yield stresses and find the failure modes, we use theaugmented-Lagrangian method and monitor the rate at which the iterations ofthe scheme converge during short-time computations ending at t = 1; iterationsconverge significantly faster when B > Bc. The results are plotted in figure 3.7.As R → 0, Bc converges to 1/√3, the limit for a slender column; as R → ∞, Bc69approaches the value 0.265, corresponding to the failure of a 2D vertical embank-ment [55, 56].0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10.250.30.350.40.450.50.55RB c (a)Numerical velocity fieldChamberlain et al.(2004)Lower boundBc for R→∞Chamberlain et al. (2003)0 0.5 1 1.5 2 2.5 3 3.5 40.20.250.30.350.40.450.50.55RB c (b)Numerical velocity fieldChamberlain et al.(2004)Lower bound IBc for R→∞Lower bound IIFigure 3.7: Critical yield stress for collapse, Bc, as a function of initial radiusR. Blue stars indicate the results from numerical simulations; the solidline shows the lower bound assuming a mechanism involving basalcollapse (labelled I); the dotted line shows the result from [21]. In (a)we show the range R < 1 and include experimental data from [19]. In(b) we show the range R < 4 and an alternative lower bound assuminga peripheral collapse (labelled II).70Figure 3.8 collects together computational results for sample failure modes,extracted from the initial velocity field for simulations with B ≈ Bc. Much asexpected on physical grounds and found for 2D dambreaks, the fluid collapsesby failing first over a basal region for thinner initial cylinders, and only at theperiphery with a wider initial state; cylinders with radii around unity show mixedtypes of failure modes.Figure 3.8: Strain-rate invariant plotted logarithmically as a density on the(r,z) plane for solutions with B ≈ Bc,Re ≈ 0, t ≈ 0, at the values ofR = 0.1,0.25,0.5,1.4,4 indicated by the r-axis limits. Also shown area selection of streamlines. In (a)−(c), the solid blue lines indicate theborder of the plastic region predicted by limit analysis (3.17). In (e),the blue line indicates the circular failure surface of the lower boundsolution in 2D geometry [53].The critical yield stress can be bounded using methods from plasticity theory71[64]. In particular, lower bounds on Bc can be established by maximizingB =∫∫wrdrdz∫∫Bγ˙rdrdz, (3.17)over families of trial velocity fields. Guided by the numerical failure modes infigure 3.8, we find a lower bound for slender initial cylinders by considering trialvelocity fields that are composed of a rigid central section in z < f (r), an overly-ing descending plug for z > g(r), and a plastically deforming region sandwichedinbetween. We take the “yield surfaces” f (r) and g(r) to have parabolic shapeand the (incompressible) plastic deformation to have a cubic horizontal velocityprofile. Thus, our trial isv =(0,−1) g(r)< z < 1(u,w) f (r)< z < g(r)(0,0) 0< z < f (r)(3.18)with u =6rg(r)− f (r)η(1−η)2w =−(6−8η +3η2)η2−6rη(1−η)2ηrη =z− f (r)g(r)− f (r)(3.19)Optimization of B can then be performed over the parameters c j of the parabolasdefining f = c1(R− r)+ c2(R− r)2 and g = c3+ c4(R− r)+ c5(R− r)2.Upper bounds on Bc can also be constructed using trial admissible stress fields[64]. However, in axisymmetric geometry, without a velocity field to determineall the stress components, an additional assumption is needed such as the vonKarman-Haar hypothesis. This artifice permitted Chamberlain et al. [18] to con-struct upper bounds on Bc. Unfortunately, this hypothesis is not appropriate forour axisymmetric slump and so their upper bound does not apply here.Note that the trial velocity field in (3.19) is continuous across the curves72z = f (r) and g(r), in line with the structure of the failure modes of figure 3.8.This contrasts sharply with the 2D problem in which failure can occur over dis-tinct curves that support velocity jumps and which become smoothed into viscousboundary layers in numerical solutions. Indeed, Chamberlain et al. [21] havepreviously computed bounds for the failure of cylinders using trial velocity fieldsmore similar to the 2D failure modes. However, distinct failure lines and viscousboundary layers do not characterize our axisymmetric solutions, leading to thechoice of the cubic velocity field in (3.19). This choice complicates the optimiza-tion computation but significantly improves the results.The lower bound on Bc for basal failure is included in figure 3.7, and compareswell with numerical computations when R is small. Figure 3.8 also displays thepredictions of the optimization calculation for the surfaces z = f (r) and g(r),which have some correspondence with the computed yield surfaces. For largervalues of R, the bound diverges from the computations and the trial velocity fieldis less similar to the actual failure modes. Both occur because of the switch in theform of the failure mode, from a basal collapse to a peripheral one.For a lower bound on Bc for peripheral failure we require a different trialvelocity field. In 2D, a useful bound is found by assuming that a circular arc offailure connects the foot of the vertical face at the edge with some point on the topsurface; above this arc, material rotates rigidly out of initial position. For a simpleaxisymmetric generalization of this trial, we again assume that a circular arc offailure arises, but divide the 2D velocity field by r (with the horizontal coordinate xreplaced by r), which ensures that the trial is incompressible. The resulting boundis shown in figure 3.7, and always lies below the 2D bound (though converges toit for R→ ∞).733.3 Comparison with experiments3.3.1 MethodsTo complement the theory, we performed experiments using an aqueous suspen-sion of Carbopol Ultrez 21 (with a concentration of about 0.5 % by weight, andneutralized by sodium hydroxide). A Herschel-Bulkley fit to the flow curve mea-sured in a rheometer (MCR501, Anton Paar, with roughened parallel plates) gaveτY = 39Pa, n = 0.3 and K = 32Pa sn.Slump tests were conducted by filling a variety of cylinders with different ge-ometry with the suspension, smoothing the fluid’s upper surface with a sharp edge.The cylinders, with radii varying from 10mm to 95mm and heights over the rangeof 11mm to 562mm, were attached to a pivoted arm that raised each container ina relatively controlled and reproducible fashion. For most of the tests, the armwas raised relatively quickly, releasing the fluid over a time of about 1 second,and the surface over which the fluid slumped was a plexiglass plate roughenedwith sandpaper. Measurements of the final deposit were taken after waiting fora few minutes. However, we also conducted tests in which we varied the rate atwhich the cylinder was raised, or replaced the roughened plexiglass with either asmooth sheet or covered it with 60 grit sandpaper. Note that for the smaller aspectratio cylinders that we used (R < 12), there was a tendency for the slumps to toppleover sideways if the cylinder was not raised sufficiently slowly (which typicallylengthened the release time to a couple of seconds), highlighting an instability oftall thin columns [8]. We abandoned any tests that showed substantial sidewaysmotion.To remove the effect of the manner in which the fluid was released, we con-ducted a second series of extrusion experiments. Here, the Carbopol was pumpedup onto the surface through a vent with a diameter of about 0.5mm. For boththe slump tests and extrusions, thickness profiles h(r, t) were extracted by takingphotographs from the side, allowing measurements of the final central height andradius (except for the extrusions over the smooth surface, as discussed below, the74tests were axisymmetric).Note that, for the slump tests, we continue to use the Bingham model to pro-vide theoretical solutions to compare against the experiments, despite the fact thatthe Carbopol is better fitted by the Herschel-Bulkley law. Thus, we implicitlyassume that the power-law viscosity plays a minor role in controlling the formof the final state of a slump. We explicitly confirmed this for isolated exam-ples of dambreak computations in which we implemented the Herschel-Bulkleymodel with a power-law index suggested from the flow curve of the Carbopol (i.e.n = 0.3, see the end of §4.2 and figure 3.12). In order to improve the comparisonwith experiments, we also used this alternative rheological model for our compu-tations of time-dependent extrusions. Specifically, the computations exploited theregularization,τ jk =[cγ˙n−1+µ2µ1(1− c)+ cBγ˙ + ε]γ˙ jk. (3.20)The characteristic velocity scale in this case is given by U = (ρgHˆn+1/K)1/n(ρ = 103 kg m−3, g = 9.81m s−2).3.3.2 Cylindrical dambreaksThe conventional slump test focuses on the distance fallen at the centre of thefluid, the so-called slump height S = 1−H∞, where H∞ is the final central depth.A summary of our results for this diagnostic is presented in figure 3.9. An impor-tant experimental parameter that we varied in the current suite of experiments isthe aspect ratio of the cylinder, R = Rˆ/Hˆ. Most previous experiments have con-ducted tests with aspect ratios close to a half, concluding that this parameter haslittle significant effect. However, the theoretical results summarized in §3 clearlyexpose how R plays an important role if varied over a sufficiently wide range.This is confirmed in our experiments, which clearly demonstrate how S decreaseswith R at fixed B; see the colour coding of the points in panel (a), and the plot forcylinders with similar initial height Hˆ in (c).The decrease of the slump height with increasing initial radius is a natural750 0.05 0.1 0.15 0.2 0.25 0.3 0.3500.10.20.30.40.50.60.70.80.91(a)R = Rˆ/HˆBS=1−H∞ 0.5 1 1.5 2GaoPashias−aPashias−bSaakClaytonCarbopol10−1 1000.250.30.350.40.45 (b) Fall vs rate of release (s−1)510Residual weight (g)0 0.5 1 1.5 2 2.5 300.10.20.30.4R (c) Fall against R1 − H∞1 − H cFigure 3.9: (a) Comparison of experimental slump height with previouslypublished data for R < 2. The points are colour coded by R, as indi-cated by the colour bar. The published data is taken from [39] (Gao;R = 12), [61] (Pashias-a: data at R≈ 12for differing materials; Pashias-b: data for red mud with varying R), [70] (Saak; R = 0.52) and [25](Clayton; R = 12). On the right, slump height is plotted against (b) therate of release (i.e. the inverse of the time taken to lift up the cylinder)for Rˆ = 3cm and Hˆ = 3.8cm, and (c) R for Hˆ ≈ 3.1cm. Most slumpswere conducted by releasing the fluid over about a second, but notrecording the time precisely; several repeated slumps with this proto-col are plotted at release rate 1.25 s−1. In (b) and (c) the fall of thecorner of the initial cylinder 1−Hc is also recorded (open circles).Panel (b) also shows the residual mass attached to the cylinder afterthe release of the fluid, and data for slumps that were conducted over asmooth surface (red stars; in these cases the slumps were left to cometo rest for five minutes).consequence of the growth of the rigid core of the dambreak and the eventualemergence of an incomplete slump. Indeed, for the latter, H∞ = 1 and S = 0,rendering this diagnostic useless for sufficiently large R. Instead, a measure ofthe degree of slump at the edge of the initial cylinder is provided by the distance76Figure 3.10: (a) Experimental and (b) numerical slump heights against Bfor R < 2. Shown in (a) are both S (filled circles) and the distancefallen by the corner of the initial cylinder 1−Hc (open circles), asillustrated by the experimental image for (R,B) ≈ (1.25,0.13) insetabove (a) (and also shown in figure 3.12). Inset above (b) are thefinal profiles (red) of the computations with R = 1.25; that matchingthe value of B of the experiment is highlighted (blue). The colourcoding of the points by R is the same as in figure 3.9, and a collectionof previous data (with R = 12) from numerical computations are alsoplotted in (b) (crosses [39]; stars [31]; pentagrams [63]).77fallen by the upper circular corner, which clearly decorates the final deposit exceptfor very low yield stresses (see figures 3.3 and 3.10). This alternative diagnostic,denoted by 1−Hc where Hc is the final height of the corner, is compared with Sin figures 3.9 and 3.10. Evidently, 1−Hc is less sensitive to variations in R. Notethat, for some of the cylinders with relatively small aspect ratio, our side imagingof the thickness profile obscures the centre of the deposit when the circular cornerfalls less far; in this situation, the measurement of S is actually given by 1−Hc.Despite the reduced sensitivity of 1−Hc to R, this diagnostic does depend onthe rate at which the fluid is released (figure 3.9(b)). Evidently, the speed at whichthe cylinder is raised affects how much fluid yields at the fluid edge, which partlycontrols the fall there. Indeed, the release rate correlates closely with the amountof material left on the cylinder after it is raised (which we measured on a weighscale in this series of experiments; see figure 3.9(b)). The impact of the releasemechanism on the degree of slump has been reported previously [39], and includesthe possibility that both inertia and adhesion to the cylinder play important rolesin the dynamics. Here, our goal is to complement our inertialess axisymmetricdambreak computations with similar, if not identical experiments, and so the effectof the release mechanism is a distraction that precludes a quantitative comparisonof theory and experiment. In particular, the action of lifting the cylinder mustforce the fluid to yield even if B > Bc, which may well explain why the existingexperimental data at higher B does not progress to the unyielded slump value S= 0(0.265< Bc < 1/√3≈ 0.577 over the full range of R).Figure 3.10 compares in more detail the experimental slump heights with theresults of our simulations. The two agree qualitatively, if not quantitatively, withthe wider initial cylinders of the simulations slumping less far. The comparisonis worse for previously reported numerical data [31, 39, 63] for Bingham slumps,which mostly collapse even further. We suspect that this is due to resolution is-sues, as illustrated in figure 3.5.Figure 3.11 shows another comparison of the final central depth and radiusof the deposit against yield stress. Here, the data is scaled by R2/3, which cor-78responds to choosing a length scale based on the initial volume of the cylinderrather than its height, and constitutes a more natural choice in the shallow limitwhere much of the memory of the initial shape is lost. The rescaling collapsesboth the height and radius data close to a common curve matching the asymptoticprediction from (3.11) (at least for the aspect ratios used in the plot, with R < 4),in contrast with the slump height diagnostic of figures 3.9 and 3.10.7900.511.522.5H∞/R2/3(a)0 0.2 0.4 0.6 0.8 10.511.52(b)B /R 2/3R∞/R2/3Figure 3.11: The central depth and radius of the final deposit against Bing-ham number, all scaled by R2/3, for R < 4. Our results for Carbopolslumps are shown by the filled circles; the crosses show data from[39]. The results from the numerical simulations also plotted in fig-ure 3.10 are shown by the solid lines; the improved asymptotic pre-diction from (3.11) is shown by the dotted line. The colour coding byR is the same as in figure 3.9.Previous experiments have also invariably been conducted over smooth sur-80faces. However, it is known [48] that spreading drops of viscoplastic fluid cansuffer effective slip unless the surface is either chemically treated or roughened.Indeed, when we conduct slumps over the smooth plexiglass, the fluid collapsesnoticeably further, implying a degree of slip (the contact angle is also noticeablydifferent, exceeding 90◦ for the roughened surface, but not for the smooth one).By contrast, slumps over a surface covered with sandpaper are similar to thoseabove our roughened plexiglass, providing confidence that our surface roughen-ing significantly reduces slip in the bulk of our tests. Figure 3.12 compares thefinal shapes of a particular slump conducted over the different three surfaces, andfigure 3.9(b) includes data for the same slump over the smooth plexiglass.810 0.01 0.02 0.03 0.04 0.05 0.0600.0050.010.0150.020.025Free slipNo slipFigure 3.12: Side images from the slump experiments and theoretical fi-nal profiles: the photographs compare slumps over the roughened,smooth and (green) sandpaper surfaces (from top to bottom); ver-tical dashed lines indicate the final radii, and the scale of the pho-tographs is indicated in the first panel, for which the rectangle showsthe initial fluid cylinder, of height 3cm and radius 3.8cm. The plotshows computations matching the experiments (R = 1.29, B = 0.13and Re = 0.28) and using the surface boundary conditions indicated(blue line for no slip, red for free slip); the dashed lines indicate theshallow-layer asymptotic predictions in (3.11) and (B.14). For the noslip case, two other computations are shown: a computation with arti-ficially high Reynolds number (Re= 35; thick grey line) and one withHerschel-Bulkley model (n = 0.3; green dots). For the free slip case,the green dots again show a computation with the Herschel-Bulkleymodel (n = 0.3).82For a theoretical examination of the effect of surface slip, we conduct simula-tions in which the no slip condition is replaced by free slip (i.e. we impose τrz = 0at z = 0). Figure 3.12 compares the result with that using no slip for simulationsin which the geometrical parameters and yield stress are matched to the experi-ments. A free slip version of the asymptotic analysis for shallow flows can also beprovided, as summarized in B.3; the final shape is predicted to remain cylindrical,with a depth of 2B√3 and a radius of R/(12B2)1/4. Both this prediction and thatof the improved no slip theory in (3.11) are also plotted in figure 3.12. The fluidspreads noticeably further with free slip, mirroring the experiments. The spreadover the smooth plexiglass is, however, rather less significant than suggested bythe free slip computations, as would be the case if there was a residual surfaceinteraction.Figure 3.12 also includes results from a simulation in which inertial effects arepromoted by taking a higher Reynolds number than that matching the experimen-tal conditions (Re is about 125 times higher in this case). The two solutions canbarely be distinguished in the figure, highlighting how inertia plays little role incontrolling the final shape in the simulations. Indeed, our computations suggestthat inertial effects are relatively minor over most of the range of physical condi-tions of our experiments: figure 3.13, shows the slump height S and sample finalprofiles for computations with R = 1 and varying B and Re. Inertia plays no rolein controlling the final shape for Re < 1; the slump height increases for higherReynolds number, but the effect is modest for our experiments (with Re< 125).830 0.05 0.1 0.15 0.2 0.25 0.300.10.20.30.40.50.6BS Re=1125103B=0.2250.1250.05rz0.5 1 1.50.51Figure 3.13: The dimensionless slump height for cylindrical dambreaks withR = 1, varying B and the Reynolds numbers indicated. The insetshows three sample triplets of the final shape at the Bingham numbersindicated.A final feature of figure 3.12 is the inclusion of computations using the Herschel-Bulkley model in (3.20), rather than the Bingham law (for both no slip and freeslip surfaces). Evidently, as alluded to earlier, the power-law viscosity has littleeffect on the final shape. However, there do appear to be some minor quanti-tive differences, especially in the vicinity of the corner stemming from the upperedge of the initial cylinder. It is conceivable that these result from differences inthe time evolution of the plugs during the collapse, which then impacts the finaldeposit. Otherwise, the viscous stress is not expected to feature directly in theforce balance controlling the final shape. Nevertheless, any differences in the fi-nal radius and height that may be introduced in this manner from the power-lawrheology are unlikely to upset the comparison of experiments and Bingham theory84in figures 9-13.3.3.3 ExtrusionsThe impact of the release mechanism on the slump is removed in experimentsin which fluid is extruded slowly from a vent onto the underlying surface, al-lowing a clearer examination of the effect of surface slip and a more quantita-tive comparison with theory. Results of a sample extrusion are shown in fig-ure 3.14, which plots the instantaneous radius R(t)/B = ρgRˆ/τY against centraldepth H (t)/B = ρgHˆ/τY , both scaled by B. Rescaling the results in this wayremoves the scaling by Hˆ implicit in the non-dimensionalization of the problemin §2, which has no meaning for the extrusions. Instead, distance is scaled by thelength scale τY/ρg.850 5 10 15 20 25 30 35012345678910 R ( t )/BH(t)/B RoughSmooth√2H /B√2H /B+ pi√3/4No-sl i pFree -sl i pFigure 3.14: Maximum radius R(t)/B = ρgRˆ/τY versus depth H (t)/B =ρgHˆ/τY for experimental extrusions with a pump rate of 31.5 ml/minover the rough (stars) and smooth (squares) plexiglass. Also plot-ted are the leading-order and improved asymptotic predictions, andHerschel-Bulkley simulations matching the physical parameters ofthe experiment with either no slip or free slip boundary conditions.Figure 3.14 includes data from an extrusion with similar pump rate over thesmooth plexiglass. The radius-height relation,R/B versusH /B, is very differentand reflective of a shallower flow. This is illustrated further in figure 3.15, whichdisplays sample images and height profiles taken from the experiments. In fact,the extrusions over the smooth plexiglass became noticeably non-axisymmetricover late times, developing interesting lobe patterns towards the rim of the ex-truded dome (see the images shown at the top of figure 3.15). It is not clearwhether these patterns reflects an intrinsic instability of sliding extensional flow86[cf. 75] or a more pronounced sensitivity to surface imperfections at the contactline.When the extrusion flux is relatively small, the developing dome above thevent is expected to be close to a steady equilibrium state, allowing us to recyclethe shallow asymptotics of §3.1. In particular, the radius-height relation from(3.11) isHB∼√2RB+pi√34, (3.21)which is also included in figure 3.14, along with the leading-order result H /B =√2R/B. The improved model in (3.21) is unrealistic for R → 0, where theinvalid treatment of the edge predicts a finite central depth, which is the analogueof the vertical cliff predicted at the fluid edge by (3.11).For equivalent numerical computations, we perform simulations as outlined in§2, but using the Herschel-Bulkley law and replacing the no-penetration condi-tion w(r,0, t) = 0 with w(r,0, t) = 2pi−1(R2v − r2)/R4v for r < Rv, where Rv is thedimensionless radius of the vent. As initial condition, we take c(r,z,0) = 0 every-where except over a shallow prewetted film spanning the vent and of depth 0.05.In these computations, Hˆ no longer has any meaning as the characteristic heightof the initial configuration, but can be defined using the net dimensional flux Qˆ: inview of our prescription for the vertical velocity, the dimensionless flux is unity,so that Qˆ = Hˆ2U = Hˆ2(ρgHˆ1+n/K)1/n. Thus, Hˆ =(KQˆn/ρg)1/(1+3n).A numerical simulation designed to match the experimental conditions is in-cluded in figures 3.14 and 3.15. The numerical simulation over-predicts the cen-tral height in comparison to the experiments, but otherwise tracks the observedradius-height relation. The comparison of surface profiles illustrates how the the-ory reproduces the shape of the observed extruded dome, but again reveals theslight discrepancy between theory and experiment. We suspect that this originatesfrom the incomplete removal of slip over the underlying roughened plexiglasssurface (fluid in the experimental extrusion flows further than in the simulationand the side profiles are less steep). However, errors in the fluid rheology (theHerschel-Bulkley fit, or non-ideal properties of the Carbopol) might also be re-87−20 −15 −10 −5 0 5 10 15 20 25 300510No slip Free slipr/Bh/BFigure 3.15: Images and profiles from the extrusion experiments also shownin figure 3.14. The top images show top and side images of an ex-trusion over the rough (left) and smooth (right) plexiglass (a rulerstraddles the extrusions to add a scale for the side images). The firstrow of plots underneath show a sequence of height profiles at simi-lar times (extruded volumes; the side photographs correspond to thesixth profiles). The insets show magnifications of the contact line.The lowest row of plots shows height profiles from correspondingnumerical simulations applying either a no slip (left) or a free slip(right) boundary condition on z = 0 (less one profile for the free slipcase). The times of the snapshots are not perfectly matched owing toirregularities in pump rate in the experiments.88sponsible.Figures 3.14 and 3.15 also include simulations results for a computation inwhich the no slip condition on z = 0, r > Rv, is replaced by free slip. As forthe axisymmetric dambreaks, the enhanced spread of the fluid and the attendantmodification in the height profile are reminiscent of how the experiments on thesmooth plexiglass differ from those above the roughened surface, reinforcing ourconclusions regarding surface slip. Once more, however, the freely sliding com-putations spread even further than the extrusions on smooth plexiglass, suggestiveof residual surface traction.3.4 Concluding remarksIn this chapter, we have presented a theoretical analysis of the axisymmetricdambreak of a cylinder of viscoplastic fluid, and compared the results with ex-periments using a Carbopol gel. In the theory, we used computations with eithera regularized Bingham model or an augmented-Lagrangian scheme to study thefluid slump, complemented by asymptotic analyses relevant for shallow gravitycurrents or tall thin columns. We also examined the states at the brink of failureto determine the conditions under which the initial cylinder does not collapse.The computations, which are based on the VOFmethod to track the fluid inter-face, are complicated significantly by the need to resolve the flow adjacent to theunderlying no slip surface. Without special attention to this detail, computationsinevitably become unresolved and spread excessively far as a result. Both this andinaccurate treatments of the yield stress likely reduce the reliability of slump testcomputations. Here, we adopted a numerical device which ensures that our com-putations remain resolved, and acts by adjusting the fluid interface and allowingthe contact line to move over the surface.For dambreaks (cylindrical slump tests), theory and experiment agree qualita-tively and are consistent with previous experiments that measure the dimension-less “slump height” (the distance fallen by the centre of the cylinder divided bythe initial height), and use a variety of different kinds of (less ideal) viscoplastic89fluids. However, one must be careful to eliminate slip over the underlying surface,which can significantly enhance the collapse. Moreover, the mechanism by whichthe fluid is released introduces quantitative differences between theory and exper-iment, either through the interaction with the lifted container, or via the amountof inertia imparted at the moment of release. The effect of the release mechanismis avoided when fluid is pumped slowly through a vent onto the surface. In suchextrusions, the agreement between theory and experiment is improved, althoughthere are some remaining differences that are most likely due either to residualslip or an inaccurate treatment of fluid rheology.A main motivation of the current work was to shed further light on the fluiddynamics of the slump test. That practical device exploits a conical initial shaperather than a cylinder, begging the question of how the present results carry over tothis other geometry. Figure 3.16 plots simulation results for the slump height com-puted for Bingham fluid with an initial shape given by an ASTM standard cone(which is 30cm high, with a top radius of 5cm and a basal radius of 10cm). Asfor the cylinder, the theoretical slumps do not collapse as far as experiments (thistime taken from [25]) at higher yield stresses, a discrepancy that could arise fromeither the release mechanism or slip. However, although the threshold for failuredepends on the slip condition on the underlying surface (Bc ≈ 0.27 or 0.32 for thisgeometry with either no or free slip, respectively), the continued significant slumpof the experimental cones for higher yield stress suggests that the release mech-anism, not slip, is more important. Except for low yield stresses, the simulationsresults also disagree significantly with those reported by [67], especially for thefailure criterion.900 0.05 0.1 0.15 0.2 0.25 0.300.10.20.30.40.50.60.70.80.91Dimensionless yie ld stress B = τY/(ρgHˆ)DimensionlessslumpheightS Clayton − conesSlender coneNumerics − no−slipNumerics − free−slipRousselFigure 3.16: Slump heights for Bingham dambreaks with an initially con-ical shape corresponding to the ASTM standard, computed using ano slip (solid line with points) or free slip (dashed-dotted line withstars) condition on the underlying plane. Squares show experimentaldata from [25]. The dashed line shows the prediction for a slenderconical column in B.4. The diamonds show the results of numericalcomputations of [67].91Chapter 4Surges of viscoplastic fluids onslopes1The purpose of this chapter is to explore models for steady, shallow viscoplasticsurges with a length that is sufficiently long that the flow converges to a steadysheet flow well upstream of the flow front; i.e. the theoretical analogue of theexperiments by Chambon et al. [22, 23, 36]. For this task, we reconsider the lu-brication analysis of Liu &Mei [51]. First, we consider the upstream extent of thesurge to examine how the surge converges to uniform sheet flow, and how the cor-responding superficial plug breaks as one progresses downstream. This demandsa variation of the lubrication analysis that is designed for almost uniform flows,and parallels theory for flow down weakly varying channels [37]. Second, oncethe upstream plugged flow gives way to a fully yielded surge with a pseudo-plug,standard lubrication theory applies; for this region, we improve that theory bycontinuing the analysis to higher order in order to account better for non-shalloweffects.We complement the shallow-flow analysis with computations using the Volume-of-Fluid (VOF) method and an augmented-Lagrangian scheme to deal with the1A version of chapter 4 has been submitted for publication (under review). [Liu, Y.], Balmforth,N.J. & Hormozi, S. Viscoplastic surges down an incline.92xzPF FFLZh(x)Y(x)ugUbFigure 4.1: Sketch of the geometry of the surge in the frame of referencein which there is no net flux and the flow is steady. The free surfaceis located at z = h; the level z = Y divides a fully yield region under-neath from either a true plug or a weakly yielded pseudo-plug. Theflow divides into three regions: a plugged flow region (PF) where thesuperficial layer of fluid is not yielded, a lubrication zone (LZ) wherethe pseudo-plug arises, and the flow front (FF) where the dynamics isnot shallow.yield stress. Such a combination of shallow-flow analysis and computation hasproven effective in our previous work studying dambreak flows and their finalshapes [53, 54], illustrated in chapter 2 and 3. Here, we examine the extent towhich the theoretical solutions match the observations of Chambon et al. and de-termine the conditions for which the pseudo-plug of lubrication theory locks upinto a true plug.4.1 Formulation4.1.1 Model equationsAs sketched in figure 4.1, we consider a two-dimensional surge of incompressibleviscoplastic fluid flowing steadily down a plane that is inclined at an angle θ to thehorizontal. We use Cartesian coordinates aligned with the plane to describe thegeometry; in the frame of reference of the surge, the inclined plane travels upslope93with a speed ub. We model the rheology of the fluid using the Herschel-Bulkleyconstitutive law. The governing equations for the velocity u = (u,w), deviatoricstress tensor τ , and pressure p are then∇ · u = 0,ρ[∂u∂ t+(u ·∇)u]=−∇p+∇ · τ +ρg(sinθ−cosθ),(4.1)and γ˙ jk = 0, τI < τY ,τ jk =(κγ˙n−1+τYγ˙)γ˙ jk, τI > τY ,(4.2)where ρ is the density, g is gravity, τY is the yield stress, the plastic viscosityµ = κγ˙n−1 introduces the consistency κ and power-law index n as two furtherrheological parameters, and τI =√12 ∑ j,k τ2jk and γ˙ =√12 ∑ j,k γ˙2jk denote secondtensorial invariants, withγ˙ =(2ux uz +wxuz +wx 2wz). (4.3)Here, subscripts on the velocity components represent partial derivatives.For the boundary conditions, we assume that there is no slip over the inclinedplane, u(x,0) = (−ub,0), and that the upper surface is stress free, so that∂h∂ t+u∂h∂x= w and (τ − pI) ·(−hx1)=(00), (4.4)on z = h(x, t). The surge ends at a flow front, x = X(t), where h→ 0, and extendsback upstream to where the flow converges to a uniform sheet flow.944.1.2 The sheet-flow solutionIn the frame of reference in which the net flux vanishes, the steady, uniform,sheet-flow solution is given by(p,τxz) = (1, tanθ)(1− zH)ρgcosθ (4.5)andu = usheet(z) =nU tan1/n θn+1××[n2n+1Y2+1/n∞ − (Y∞− zH )1+1/n], 0< zH B,0=∫ h0u(x,z) dz,γ˙xx = 2εux, γ˙xz = uz+ ε2wx, γ˙ =√γ˙2xx + γ˙2xz,(4.11)with u(x,0) =−ub andτ +hx(p− εσ) = 0p+ εσ + ε2hxτ = 0}on z = h. (4.12)4.2.1 Plugged FlowIn the PF region, the flow is nearly uniform, so we introduce the asymptotic se-quences,h = 1+ εh1+ ..., p = 1− z+ ε p1+ ..., (4.13)τ = S(1− z)+ ετ1+ ..., u = u0+ εu1+ ... (4.14)and σ = σ0+ ... The leading-order terms from (4.11) recover the velocity profileof the sheet-flow solution, this time in dimensionless form:u0 =−ub + 12S×{z(2Y∞− z), 0< z z > Y . This demands that u1(x,Y ) = 0, orT =−12Y∞h1x. (4.23)Finally, we impose the flux constraint,∫ h0udz = 0 (4.24)which gives, at O(ε),ubh1 =12SY 2∞h1+16h1xY3∞ +12TY 2∞. (4.25)Hence,h1x = 2Sh1 and Y1 =Y∞h1, (4.26)and so the departure from the uniform sheet solution grows exponentially in thedownslope direction, with an exponent given by 2S.4.2.2 Breaking the plugReturning to (4.17), we now observe that, over the plug,h1x = 2σ0x + τ1z, (4.27)which can be integrated from z =Y∞ upto z = 1 to give∂∂x∫ 1Y∞σ0 dz =12Sh1(1−Y∞). (4.28)Hence ∫ 1Y∞[σ0(x,z)]x−∞dz1−Y∞ =14h1, (4.29)which, given that |h1| grows exponentially, implies that the net jump in extensionalstress across the plug must also increase towards the flow front. Moreover, since100h1 must eventually become O(ε−1) to curve the surge, the plug must inevitablybreak.A further constraint is provided by the unyielded condition of the plug (giventhe leading-order shear stress τ ∼ S(1− z)),−√B2−S2(1− z)2 < σ0 <√B2−S2(1− z)2, (4.30)which bounds the integral on the left of (4.29). In particular, that integral cannotexceed 12piB in absolute size. Therefore, the plug must have broken when|h1|> 2piB. (4.31)4.2.3 Lubrication ZoneThe preceding analysis indicates that the plug breaks when h−1= O(ε). Down-stream, the departure from the sheet-flow solution grows further as we enter thelubrication zone and h− 1 becomes O(1). With h = h(x) 6= 1, the leading-ordersolution of (4.11) now furnishes the standard lubrication result for the velocityprofile,u∼−ub + 12(S−hx)×{z(2Y − z), 0< z 4/X0; by contrast, thevertical edge of the rectangular block ensures the fluid always collapses in theshallow limit. The final state is given by Y → 0, which leads to (2.8).The shallow-layer theory is the leading order of an asymptotic expansionwhich, the current dimensionless scalings, corresponds to the limit B ≪ 1. Forthe final shape one can go further with the analysis and compute higher-order cor-rections in the effort to extend the accuracy of the approximation. In particular,following the analysis in [33] but bearing in mind that the slump comes to rest in135a state of horizontal expansion, we findp≈ h− z−√B2− v2, v =−hx(h− z), (A.11)τxx ≈√B2− v2− vτ1√B2− v2 , τxz ≈ v+ τ1, (A.12)whereτ1 =− ∂∂x[1hx(v√B2− v2+B2 sin−1 vB)]. (A.13)Imposing the lower boundary condition, τxz = B at z = 0, and integrating in x thengives12h2−h√B2−h2h2x−B2hxsin−1hhxB=C−Bx, (A.14)where C is an integration constant. Evaluating the higher-order corrections in(A.14) (i.e. the second and third terms on the left-hand side) using the leading-order approximation hhx =−B leads to (2.9) with C = BX∞−pi2B2/8.A.3 Slender columnsWhen the column of viscoplastic fluid remains slender throughout its collapse,we may use the thin-filament asymptotics outlined by [13]. The key detail is thatthe horizontal gradients are much larger than the vertical ones and the verticalvelocity greatly exceeds the horizontal speed. Moreover, because the sides of thecolumn are stress free, shear stresses must remain much smaller than the exten-sional stresses and the vertical velocity cannot develop significant horizontal shearand remains largely plug-like. These consideration indicate that (see [13])w≈W (z, t), u≈−xWz, (A.15)andp≈ τxx =−τzz ≈ B−2Wz, (A.16)136the latter of which follows from the leading-order horizontal force balance (whichis (τxx− p)x ≈ 0) and constitutive law (given γ˙ ≈ 2|Wz|). The width-averagedmass conservation equation and vertical force balance then imply [13]ξt +(ξW )z = 0 and ξ +2(ξ p)z = 0, (A.17)where x = ξ (z, t) is the local half-width.The equations in (A.17) can be solved analytically by transforming to La-grangian coordinates (a, t), where a denotes initial vertical position (i.e. themethod of characteristics; cf. [49]). For the rectangular or triangular blocks, wehave the initial condition ξ (a,0) = X0 or X0(1− 12a), respectively. The transfor-mation then indicates thatξ (a, t)∂ z∂a={X0X0(1− 12a)and ξt(a, t) =−ξWz. (A.18)Henceξ p = 12X0(1−a/a∗), (A.19)given that p = 0 at the top of the column where a = a∗ = 1 or 2. If the fluid isyielded, the constitutive law impliesp = B−2Wz = B+ 2ξtξ. (A.20)After a little algebra and the use of the bottom boundary condition z(a= 0, t)=0, we findξ (a, t)ξ (a,0)= E +1−E2B(1− aa∗), E = e−t/2B, (A.21)andz =2a∗B1−E log[1−E +2BE(1−E)(1−a/a∗)+2BE]. (A.22)The yield condition in this limit becomes p < −B, or a < a∗(1− 2B), which137translates toz < Z(t) =2a∗B(1− e−Bt/2) log[1− (1−2B)e−Bt/22B]. (A.23)The column does not therefore yield anywhere when B > 12. If fluid does yield,the base spreads out to a distance ξ (0, t)→ X∞ = X0/(2B), and the column fallsto a height H∞ = 2a∗B[1− log(2B)]. The yield condition and runout are the samefor both the rectangle and triangle because, in the slender limit, all that matters isthe weight of the overlying fluid. Sample solutions are shown in figure A.1.The main failing of the slender asymptotics is that the no-slip boundary con-dition is not imposed: the fluid slides freely over the base, leading to the collapsedcolumn being widest at z = 0, whereas the fluid actually rolls over the base in atank-treading motion (cf. figure 2.8). This failing must be remedied by adding aboundary layer at the bottom (with different asymptotic scalings). In any event,we attribute the lack of agreement between the slender asymptotics and the nu-merical results in figure 2.8 to this feature.A.4 Plasticity solutionsIn this appendix, we summarize slipline and bound computations based closely onexisting work in plasticity theory [20, 21], emphasizing some minor generaliza-tions for triangular initial shapes and circular failure curves. The reader is referredto [20, 21] for further details of the basic developments.A.4.1 Slipline fieldsThe slipline solutions of Chamberlain et al. [20] begin from the side free surfacewhere the stress field is specified and are constructed as follows: settingP≡ pB+zB& (τxx,τxz) = B(−sin2θ ,cos2θ), (A.24)1380 0.05 0.100.511.52x(a) B=0.10.02 0.04x(b) 0.3Figure A.1: Slender asymptotic solutions for (a) B = 0.1 and (b) B = 0.3,for X0 = 0.025, at times t = 0, 4, 9, 16, 36, 100 and 1000. Collaps-ing rectangles (triangles) are shown by solid (dotted) lines; the dotsindicate the final profiles.139the side boundary conditions implyP = 1+ zˇ & θ =3pi4+φ (A.25)on xˇ = −Szˇ, where (xˇ, zˇ) = B−1(x,z), and (S,φ) = (0,0) for the rectangle and(X0/2, tan−1 12X0) for the triangle. On the α−characteristics,P+2θ = constant,dzdx= tanθ ; (A.26)for the β−characteristics,P−2θ = constant, dzdx=−cotθ . (A.27)Beginning from the section of the side, 0≤ zˇ≤ zˇP, with zˇP a parameter, the char-acteristics can be continued into the fluid interior using a standard finite differencescheme to solve the characteristics equations in (A.26)-(A.27) [20, 64]. Below theresulting web, an expansion fan is then added that spreads out from the base point(x,z)= (X0,0)with θP≤ θ ≤ 3pi/4+φ , where θP is a second parameter (cf. figureA.2). The combined slipline field is then continued to x= 0, or xˇ=−X0/B. At thispoint, the two characteristics that bound the complete slipline field must cross andterminate with θ = 3pi/4, in view of the symmetry conditions there. This selectsthe two parameters zˇP and θP. Finally, along the uppermost α−characteristic, thetotal vertical force must match the weight of the overlying plug, which translatesto imposing the condition,X0− 12SB2zˇ2P = B2∫ zˇP0cos2θ dzˇ+B2∫ 0X0/B(P− sin2θ)dxˇ, (A.28)and determines the relation B = Bcrit(X0) (as plotted in figure 2.10). Examples ofthe slipline field are shown in figure A.2. Note that the sliplines of the two familiesbegin to cross over one another near the top of the expansion fan if X0 is increasedpast some threshold [20]. At that stage a curve of stress discontinuity must be1400 0.1 0.200.10.20.30.40.50.60.70.80.91yx(a)X0=0.20.1 0.2 0.3 0.4 0.5x(b) X0=0.50 0.5 100.20.40.60.811.21.41.61.82x(c) X0=1Figure A.2: Slipline solutions for a rectangle with (a) X0 = 0.2 and (b) X0 =0.5, and a triangle with X0 = 1. The dashed and dotted lines indicatethe lower bound and its improvement of A.4.2.introduced to render the slipline field single-valued. We avoid incorporatingthis detail here and only provide slipline solutions without a discontinuity, whichlimits the triangle data in figure 2.10(b) to X0 < 1. At still higher X0 > 2, theconstruction fails for the triangle altogether because the α−characteristic fromthe base of the free surface proceeds into z < 0.A.4.2 Simple failure modesFor the circular failure surface of a relatively wide initial rectangle, we refer thereader to existing literature (e.g. [24]). Here, we summarize the computation for141slender rectangles and triangles.We first consider the case where failure occurs on straight lines, as in [21].As illustrated in figure 2.11(b–c), we introduce two lines of failure, with slopestanγ and tanζ , that divide the initial block up into a lower stationary triangle, anintermediate triangle that slides out sideways, and the residual overlying materialthat falls vertically. When the downward speed of the top is W , the continuity ofthe normal velocity across each failure line demands that the intermediate triangleslides out parallel to the lower failure line with velocity U(cosγ,−sinγ), whereU = W secζ/(tanγ + tanζ ). Let AI and AII denote the areas of the intermediatetriangle and top block and LI and LII be the lengths of the lower and upper failurelines, both respectively. Equating the plastic dissipation across the failure lineswith the release of potential energy then furnishes (in our dimensionless notation)B(ULI +WLII secζtanγ + tanζ)= AIU sinγ +AIIW. (A.29)Geometric considerations allow us to express all the variables in terms of X0 andthe two angles γ and ζ . Equation (A.29) can therefore be formally written in thesuggestive form, B = B(γ,ζ ;X0). We then optimize the function B(γ,ζ ;X0) overall possible choices of the two angles (γ,ζ ) to arrive at the bound Bc(X0). It turnsout that Bc =12tanγ = 12tanζ = 12(√1+X20 −X0) for the rectangular block [21].In the case of the triangle, tanζ =√2(1+ tan2 γ)− tanγ , leaving a straightfor-ward algebraic problem to solve for the optimal γ and Bc, with solutions shownin figure 2.11(c) and 2.10(b). The failure lines of these bounds are compared tothe sample slipline fields in figure A.2, illustrating the manner in which the boundattempts to capture the actual plastic deformation.The streamlines of the numerically computed failure modes suggest that thepreceding bounds might be improved if the triangle at the side were allowed torotate out of position rather than slide linearly. In this situation, the failure sur-faces become circular arcs rather than straight lines which complicates the formof the power balance corresponding to (A.29) and the geometrical constraints.142Three optimization parameters are required to define the circular arcs; we use thelocal slopes at the bottom corner, sα , and midline, sβ and s, as illustrated in figure2.11(b,c). The optimization problem can then be continued through with the helpof the computer. We use the built-in function FMINSEARCH of Matlab to per-form the optimization of B(sα ,sβ ,s;X0) and improve the bounds on Bc(X0). Thecircular failure arcs corresponding to the three slipline solutions of figure A.2 areagain included in that picture.Note that the bounds for the triangle predict that sα falls to zero for X0 > 2.8with a straight failure surface and X0 > 1.2 for rotational failure. For wider initialstates, this parameter must then be removed from the optimization, which makesthe bounding procedure less effective. A more general and effective construction,that retains sα as a parameter, is to allow the lower circular failure arc to intersectthe base for x < X0, but not pass through that surface, and then continue beyond.That is, we allow the arc to proceed through a minimum at z= 0 and then intersectthe side surface at a finite height (cf. figure 2.13(b,c)). This extension permitscomputations of improved bounds for arbitrarily wide triangles and is plotted infigure 2.10(b). Figure 2.13(b,c) illustrates how the resulting arcs compare wellwith the computed failure modes for moderate width. Even for the widest trianglewith X0 = 8, the bound (Bc > 0.1635) is close to the computed value of Bc ≈0.1642. For X0≫ 1, the bound converges to Bc > 3/X0. By contrast, the shallow-layer asymptotics predict failure for Bc ∼ 4/X0 (see A.2), indicating that there isfurther room for improvement in this limit.143Appendix BAxisymmetric viscoplasticdambreaks and the slump testB.1 Numerical convergence studyFigure B.1 shows the results of a resolution study for the runout of a Newtonianslump with R = 1 and varying grid size. The first panel presents results for thePLIC scheme with the non-conservative correction of [53]; those with the con-servative correction are plotted in the second panel. Both are compared with theresults of lubrication theory [11, 47]. Also included in figure B.1(b)-(c) are corre-sponding results for the runout of a Bingham slump with B = 0.05 and a compar-ison of the fluid interfaces of the various solutions at t = 250, all for the case ofthe PLIC scheme with the conservative correction (and both values of B).Figure B.2 shows further details of the resolution studies with B = 0 and 0.05,plotting the runout and central depth of the slump at t = 250 for various compu-tations with differing resolution. Computations both with and without the correc-tion schemes are presented; the results for the original PLIC algorithm convergemuch more slowly than those for the corrected schemes, with the conservativecorrection scheme being superior. For B = 0.05, the three different treatments ofthe interface are shown for a computation using the regularization method to deal1441.52R(t) (a)1/401/801/1601/320Shallow0 50 100 150 200 25011.52tR(t)(b)rz(c)0 0.5 1 1.5 200.20.4Figure B.1: Computations with varying grid size (as indicated). Panels (a)and (b) show the flow front as a function of time for a Newtonianslump with R = 1 using the PLIC scheme with the non-conservativeand conservative corrections, respectively. Also inclued in (b) are cor-responding results for a Bingham slump with B = 0.05 (green curves).Panel (c) shows the interfaces at t = 250 of the computations in (b)(with inset showing magnifications near the flow fronts). The red cir-cles show leading-order shallow-layer solutions [9, 47].14540 80 160 3201.81.9Grid points pe r uni t l ength (a) R (250)B=0.0540 80 160 3200.460.48 (b) H (250)B=0.0540 80 160 3200.230.24 B=040 80 160 3202.42.6 B=0Figure B.2: (a) Runout and (b) central depth at t = 250 with varying gridsize for a Newtonian (B = 0; top) and Bingham slump (B = 0.05;bottom) using various interface-tracking schemes and Navier-Stokessolvers: regularization with the PLIC algorithm including the non-conservative (crosses), conservative (stars) or no correction (circles)scheme, or using the augmented-Lagrangian method with PLIC andthe conservative correction (squares). There is no difference betweenthe Navier-Stokes solvers for B = 0.146with the yield stress. Also shown is a series of computations using the augmented-Lagrangian method and the PLIC scheme with a conservative correction, whichillustrates how the two methods for the dealing with the yield stress yield similarresults.B.2 Higher-order shallow asymptoticsWhen flow becomes arrested, throughout the bulk of the fluid the constitutivemodel reduces to the plasticity law τi j = Bγ˙i j/γ˙ . In combination with the rescaledequations in (3.7) and (3.10), one may then deduce the leading-order stress com-ponents, τˇrz = Bˇ(1− z/h0),(τˇrr, τˇθθ ) = Bˇ(uχ ,χ−1u)√1− (1− z/h0)2√u2χ +u2/χ2+uuχ/χ,(B.1)where h0 =√1−χ/χ∞ and u are the leading-order profile and radial velocity,with χ∞ ≡ εr∞. If we instead neglect only the O(ε2) terms in (3.7)-(3.8), we findthat p+ ετˇzz = h− z+O(ε2). An integral of the radial force balance over thedepth of the fluid then furnishes the equationhhχ + Bˇ =3ε2ddχ[∫ h00(τˇrr + τˇθθ )dz]+ε2χ2ddχ[χ2∫ h00(τˇrr− τˇθθ )dz]. (B.2)To correct the profile, we evaluate the O(ε) terms on the right of (B.2) using theleading-order stress components.In order to accomplish this, however, we must first determine the leading-order radial velocity u in order to fix the indeterminacy of the axisymmetric stresscomponents. The convergence of the leading-order slump solution to rest hasbeen explored previously in [11], extending the work of Matson & Hogg [58] for2D dambreaks. It was found that, for t → ∞, the radial velocity u is plug-like1470 0.2 0.4 0.6 0.8 100.20.40.60.81r /R(τrr,τθθ,τzz)/B τ z zτ r rτ θ θ(a)0 0.2 0.4 0.6 0.8 100.20.40.60.81Plug speed(b)r /R NumericsAsymptoticsFigure B.3: (a) Scaled stress components (τrr,τθθ ,τzz)/B and (b) radialspeed u/umax plotted against r/R along the level where τrz = 0 atthe end of the computation of a simulation with B = 0.01 and R = 1.throughout the fluid depth and given by u = Q2ξ/(4ξ ), whereddξ[(1−ξ 2)Q2ξ]= 2λ (1−ξ 2)(1−Q), (B.3)with ξ =√1−χ/χ∞ and λ ≈ 23.3855. With the solution of this problem in hand,we may then compute τˇrr and τˇθθ to leading order. The results are compared withcomponents extracted from the numerical simulations in figure B.3(a) for a sampleslump near the arrest of motion. Plotted are the stress components along the levelwhere τrz = 0 (which is z = h0 in the leading-order asymptotics), after recastingthe asymptotic predictions in terms of the original variables. Except near thecentre and edge of the collapsed cylinder (where the shallow-layer asymptoticsbreaks down), the two agree qualitatively. Evidently, τzz is almost constant inradius, whereas τrr and τθθ match one another up to a small correction.The approximation τˇrr ≈ τˇθθ (which corresponds to a plug speed u ∝ r; cf.figure B.3(b)), offers a simpler analytical pathway to a corrected final profile: theyield condition now implies thatτˇrr + τˇθθ ≈ 2Bˇ√3√1−(1− zh0)2, (B.4)148and so the corrected profile now follows from (B.2) ashhχ + Bˇ∼ pi√34εBˇhχ , (B.5)which provides the solution quoted in §3.2.1.B.3 Shallow slippy flowsWhen the fluid slides more freely over the underlying surface, the lubricationmodel of [51] no longer captures the leading-order dynamics. In particular, theextensional stress components (τrr,τθθ ,τzz) become promoted to higher orderin comparison to the shear stress τrz because the reduced surface traction gen-erates little vertical shear and the horizontal flow becomes plug-like throughoutthe fluid layer. In this circumstance, the relevant asymptotic description is that ofa viscoplastic membrane model [4]: rather than rescaling as in §3.1, we set onlyτrz = ετˇrz. Then, ignoring inertia, the leading-order momentum equations are0=−pχ + 1χ(χτrr)χ + τˇrz,z− τθθχ+O(ε), (B.6)0=−pz + τzz,z−1+O(ε), (B.7)With the rescaled horizontal velocity u∼ ε−1uˇ(χ , t)≫ w(χ ,z, t), the constitutivelaw predicts that(τrr,τθθ)∼ 2(1+Bγ˙)(uˇχ , uˇ/χ), (B.8)γ˙ ∼ 2√uˇ2χ + uˇuˇχ/χ + uˇ2/χ2, (B.9)which are independent of z. The surface boundary conditions in (3.8) are replacedbyτˇrz+hχ(τrr− p)p+ τzz}= O(ε) on z = h. (B.10)149The integral of (B.7) now furnishes p+ τzz = h− z to leading order. From thedepth integrated continuity equation and (B.6), it then follows thatht +1χ(χhuˇ)χ = 0, (B.11)0=−hχ − τbh+1h[h(2τrr + τθθ )]χ +τrr− τθθχ. (B.12)Here τb = τˇrz(r,0, t) is the basal shear stress, which must be specified by an addi-tional model of slip.For a cylindrical dambreak with free slip (τb = 0), the equations admit a simplesolution in which the fluid retains the shape of a cylinder, with h = H (t) anduˇ = χU (t) for χ < εR(t), and R˙ ≡RU . The extensional stresses become τrr =τθθ = 2U +B/√3, and mass conservation demands that R2H = R2. Although(B.12) is now satisfied, the combined gravitational and extensional stress must stillvanish at the edge of the cylinder, which implies that −12H 2+H (2τrr + τθθ ) =0, or 12H = 6U +B√3. Thus,H˙ =−2H R˙R= 16R√H (B√3− 12H ). (B.13)Provided B < 1/(2√3), the cylinder therefore collapses to a final state given byH = 2B√3 and R =R(12B2)1/4. (B.14)B.4 Tall slender conesMost practical slump tests take a truncated cone as the initial shape with the radiusat the top equal half of that at the base. The slender analysis in this case has beenderived by Clayton et al. [25] using the Tresca yield condition. For von Mises,150the corresponding result iss = 1−h0−√3B ln( (1+α−1)3−1(1+α−1h0)3−1), (B.15)where h0 is the height of the unyielded region, given byB =α3√3(1+α−1h0)3−1(1+α−1h0)2and α = Rtop/(Rbase−Rtop).151Appendix CSurges of viscoplastic fluids onslopesC.1 Additional numerical detailsIn this appendix we provide further details of the numerical computations. We usethe dimensionless version of the problem in which lengths are scaled by H andvelocities by U as in §4.1.To find the steady surge states, we solve suites of initial-value problems inthe frame of reference of the surge, following the strategy outlined in [53, 54] orChapter 2 and 3. We begin from initial conditions in whichmotionless viscoplasticfluid is deposited on the inclined plane with a rectangular shape whose depthis set by the sheet-flow solution. The front of the rectangle is smoothed overa streamwise scale of order of a fraction of the fluid depth using a arc tangentfunction. The simulations do not appear to be sensitive to the initial condition,at least at the Reynolds numbers chosen for most of our simulations (§C.1.2),with no indication of multiple equilibria or unsteady states (but see the discussionsurrounding figure 4.9).The initial-value problem is solved exploiting the PLIC (Piecewise Linear In-terface Calculation) algorithm to evolve the interface within the volume-of-fluid152scheme. The essential details of the algorithm can be found in [57], although wemodify it slightly to surmount a numerical difficulty arising from an unresolvedlayer of the upper-layer fluid coating the inclined plane, as is described in [53, 54]or Chapter 2 and 3. In brief, the modification amounts to monitoring the vol-ume fraction adjacent to the inclined plane and adjusting the value of c there if itexceeds a threshold (set to 0.99). This replacement of ambient fluid with yield-stress material destroys mass conservation, which is restored by rescaling the flowheight uniformly over the length of the surge, incurring a further error, of order afraction of a grid spacing.The evolution observed in these initial-value problems suggest that the finalsteady states possess no contact line along the underlying plane. Rather, a contin-ually thinning finger of ambient fluid coats the plane as one progresses upstream.Therefore, although the adjustment to the PLIC algorithm applies an approxima-tion to allow the contact line of the initial condition to migrate with the plane tocreate this finger, once that feature is established, the scheme simply amounts toneglecting the small amount of ambient fluid within the finger once it becomesthinner than the lowest grid cell. This prevents interpolation errors in the veloc-ity field and shear stress of the finger from excessively lubricating the surge, andpermits the computation to otherwise remain resolved, as we now document.C.1.1 Resolution studyFigure C.1 shows the results of a resolution study for the profile of a Binghamsurge with θ = pi/18 and B = 0.0522 (Re = 1, ub = 0.041 and the total area offluid is 24). At the rear of the surge, a no-slip back wall is imposed, and the gridspacing is uniform with ∆x = 4∆z. Varying that grid spacing by a factor of 8 (asindicated) furnishes barely any discernible difference in the free surface profile.Indeed, the root-mean-square difference in the velocity field, defined as√∫∫ |u−u1|2c1 dxdy∫∫c1 dxdy,153between the coarsest and finest of these simulations is about 0.088ub, and de-creases to 0.025ub between the two finest simulations. Here, c1 and u1 denotethe reference solution, which is that for the finest simulation, interpolated ontothe grid of the coarser solution. Also shown are selected contours of constantinvariant τxz (which closely match those of τI below the plug), illustrating theconvergence with resolution, at least for stress levels sufficiently in excess of B.Contours closer to B show significantly less degree of convergence, as highlightedby the roughness of the yield surfaces plotted in figure 4.4, and the stress solutionremains sensitive to resolution for τI < B. These latter deficiencies are not prob-lematic as the solution is independent of the stress state over the plug as long asthe fluid there is not yielded. Note that the computations reported in §4 and theremainder of this appendix all use the finest grid of the resolution study.C.1.2 Inertial effectsFigure C.2 shows numerical simulations of Bingham surges with varying Reynoldsnumber Re= HU /κ from 0.1 to 10 (which span the range of all the simulationsreported in this study), with a back wall providing the left-hand boundary condi-tions. With the scalings of the problem outlined in §4.1, the dimensionless prob-lem retains the Reynolds number only as a factor in front of the inertial terms.Consequently, the flow profile becomes independent of Re in the inertialess limit.Indeed, the flow profiles and stress levels shown in figure C.2 closely collapse(save for the relatively rough yield surface) and the root-mean-square differencesin the velocity field and stress invariant are less than 4×10−3ub and 0.01B, respec-tively. Thus we conclude that inertial effects are not significant. The computationsreported in §4.1 are conducted with Re = 1.C.1.3 The effect of the back wallTo gauge the effect of the back wall on the computations, figure C.3 comparestwo simulations with different boundary conditions imposed along x = 0. In thefirst, the velocity field corresponding to the uniform sheet flow is imposed (so154Figure C.1: Simulations of a Bingham surge profile with θ = 10◦ and B =0.0522 (Y∞ = 0.7, fluid area of 24, ub = 0.041 and Re = 1) with thegrid sizes indicated. (a) shows h and Y ; (b) a magnification near theflow front; (c) contours of constant τxz/B, as indicated (with the surgeprofile of the finest resolution case shown by the lighter grey line);(d) the shear stress along z = 0.025 (the lowest grid cell of the coars-est computation). In (a), the inset shows the flow front X f , mean shearstress, τxz/B=∫ ∫(τxz/B) cdxdz, and τxz at the point (x,z)= (18,0.4),all plotted against the vertical grid spacing ∆z (∆x = 4∆z). The con-vergence of τxz is impeded by the need to resolve the sharply localizedregion of high stress underneath the flow front (cf. figure 4.3 and pan-els (c) and (d)).155Figure C.2: Simulations of Bingham surges for θ = 10◦, ub = 0.041 andB = 0.0522 (Y∞ = 0.7) for the Reynolds numbers indicated (the totalarea of fluid is 24). Plotted are h, Y (τxz = B) and the stress levelsτI/B = 1, 2 and 3.(u,w) = (usheet ,0)); for the second, we impose a no-slip condition, u = w = 0 Thefigure displays the stress invariant τI and plug regions; the solutions are much thesame except within a region near the back wall in which flow adjustments arisedue to the boundary condition there (and which changes the flow length slightly).In all our simulations, the extent of these flow adjustments was restricted to along-slope lengths of about 4H.156C.2 Improved lubrication solutionFor Herschel Bulkley fluid, the dimensionless system ispx = εσx + τz +S, pz =−εσz + ε2τx−1,0= τ(x,h, t)+ [p(x,h, t)− εσ(x,h, t)]hx,0= p(x,h, t)+ εσ(x,h, t)+O(ε2),(στ)=(γ˙n−1+Bγ˙)(2εuxuz + ε2wx)for 0< z < Y,B2 = σ2+ τ2 for Y < z < h,γ˙ =√4ε2u2x +(uz + ε2wx)2,ht =−(∫ h0udz)x.In the fully yielded region, uz ∼ O(1), σ ∼ O(ε) and γ˙ ∼ O(1), and droppingthe O(ε2) terms then givespz =−1 and px = τz+S (C.1)Hencep = h− z+Pτ = (S−hx−Px)(Y − z)+BY = h+T −BS−hx−Px(C.2)where P = P(x) and T = T (x), which further imply thatτ = B+unz ,uz = (S−hx−Px) 1n (Y − z) 1n ,u =nn+1(S−hx−Px)1n[Yn+1n − (Y − z) n+1n]−ub,(C.3)all to O(ε2).157In the pseudo-plug, u = up(x)−ub + εu1(x,z), and soγ˙ = ε√4u2px +u21z +O(ε) = εΓ.Hence, (στ)=BΓ(2upxu1z)+O(εn), (C.4)if n < 1. Thusσ =√B2− τ2+O(εn) (C.5)Force balance over this region demandspx = S+ τz+ εσx & pz + εσz =−1+O(ε2), (C.6)which, given the surface stress conditions, now providep = h− z− εσ +O(ε2),τ = (S−hx)(h− z)+2ε(∫ hzσdz)x+O(ε2),(C.7)or, given (C.5),τ = (S−hx)(h− z)+ 12εB2(2θ + sin2θS−hx)x+O(εn+1),whereθ = sin−1(S−hx)(h− z)B. (C.8)Next we observe that P ∼ O(ε2), because p = h− z+P and σ = O(ε) in thefully sheared region, but p = h− z− εσ +O(ε2) in the pseudo-plug. The matchof τ = (S−hx)(h− z)+T in z < Y with(S−hx)(h− z)+2ε(∫ hz√B2− τ2dz)x+O(εn+1)158for z >Y , then demands thatT = 2εB(∫ hY√1− (S−hx)2(h− z)2B2dz)x= 12εpiB2(1S−hx)x.(C.9)Now we match the velocity profile of the fully yielded region in (C.3) with that ofthe pseudo-plug u = up−ub + εu1, to find u1(x,Y, t) = 0 andup =nn+1(S−hx)1nYn+1n .Finally we compute u1 and the downslope flux: in the pseudo-plug,u1z2upx=τσ=(S−hx)(h− z)√B2− (S−hx)2(h− z)2+O(ε), (C.10)and so, given that (S−hx)(h−Y ) = B+O(ε),u1 = 2upx√B2− (S−hx)2(h− z)2S−hx +O(ε).The flux can then be computed as∫ h0udz =∫ Y0udz+∫ hY(up−ub + εu1)dz= up(h− n2n+1Y)−ubh+ 12εpiB2upx(S−hx)2 .The equations of the improved model quoted in the main text now follow, ontaking n = 1.159C.3 Improved slump profilesCuriously, the improved model for slumped shapes on slopes given by (4.38) im-plies thathhx ∼−12εpiB2(h−1x )x or h∼ [3εpiB2(x f − x)]1/3,for x → x f . This contrasts sharply with the solution of the iterated version ofthe model in (4.41) which has a finite depth at the edge. Evidently, the freedomafforded by the extra derivative allows us to reach the flow front with h → 0 andhx →−∞. This feature of the improved asymptotic theory was not appreciated inour earlier papers [53, 54] or Chapter 2 and 3, where the simpler iterated versionof the model was implemented.Figure C.4 compares the various asymptotic solutions, with or without a back-ground slope. The two versions of the improved model differ by O(ε2) away fromthe flow front; within a distance of O(ε2) of x = x f , however, the flow depths be-come O(ε) different, permitting the iterated version to terminate at finite depth.For comparison, figure C.4(a)-(b) also includes simulation data from figure 14 of[53] for dambreaks on a horizontal surface with either square or triangular initialconditions. The (red) corners visible in panel (b) are relics of square initial condi-tions, whereas the (blue) sharp spikes at the back evident in panel (a) are remnantsfrom triangular initial conditions (in the scaled coordinates the slumps have dif-ferent lengths). Further computations for the slump of a rectangular block on anincline (with unit height and an upstream back wall) are included in figure C.4(c)-(d). Both comparisons indicate that the two versions of the improved model out-perform the leading-order theory away from the flow front. Near that steep fea-ture, the smooth decline to zero thickness of the non-iterated model provides aslightly more satisfying comparison with simulations. However, the asymptotictheory is not valid at the flow front where the slope of the free surface diverges.Moreover, the surface in the numerical simulations eventually overturns to createmulti-valued profiles with a finite elevation at the leading edge. Consequently, itis not clear which version of the improved model is superior; the iterated model160(which we have employed previously, and continue to use in the main text) hasthe advantage of being the simpler.161Figure C.3: Simulations for θ = 10◦ and B = 0.0354 (Y∞ = 0.8) with dif-ferent left-hand boundary conditions: (a) (u,w) = (usheet ,0) and (b)(u,w) = (0,0). Shown is the second invariant τI as a density onthe (x,z)−plane. The (true) yield surfaces are indicated by the solid(green) lines. The simulation in (a) corresponds to the truncated so-lution shown in figure 4.3. The profiles and the true and fake yieldsurfaces are compared in (c).162Figure C.4: Slump profiles for (a)-(b) S = 0 and (b)-(c) S 6= 0. For (a)-(b), scaled variables are plotted, which eliminates any free param-eters; in (c)-(d) the profiles are shown for 12εpiB = 0.2. The darksolid lines show the solutions to (4.38), whereas the dotted lines showthe solutions to the iterated version in(4.40); the dashed lines showthe leading-order approximation. In (a)–(b) the lighter (red and blue)lines show the results of a series of simulations from [53] with ε = 1,B = 0.02, 0.03, ..., 0.1 and the flow fronts aligned. In (c)–(d) thelighter (red) lines show additional simulations of the slump of a rect-angular block on an incline with B = (0.1,0.2,0.3,0.4)/pi and θ = 5◦;the inset shows the aligned, but unscaled profiles.163