The Propagation of the GravityCurrent of Bingham FluidA Numerical InvestigationbyYe LiuB.Sc., Peking University, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2014c© Ye Liu 2014AbstractWe have studied the propagation of 2D unit block of viscoplastic fluid ofBingham type over a horizontal plane, underneath another Newtonian fluid.We numerically simulate the dynamics of a two-layer fluid in a rectangledomain, using the volume-of-fluid method to deal with the evolution of theinterface, and regularization scheme of the constitutive law, which replacesunyielded plugs with very viscous flow. We explore the final shape of theflow for varying yield stress, comparing the numerical results with the pre-dictions of the asymptotic theory, a plasticity model based on slipline theory,and other past results. Numerical difficulties with the moving contact linesare encountered during the numerical simulation. A slip boundary condi-tion is used to address this issue, the validity of which should be furtherinvestigated.iiPrefaceThis thesis is original, unpublished, independent work by the author, YeLiu.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . ixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Literature review . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Problem Description . . . . . . . . . . . . . . . . . . . . . . . . 42.1 Problem setting . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.3 Dimensionless form . . . . . . . . . . . . . . . . . . . . . . . 63 Asymptotic Results . . . . . . . . . . . . . . . . . . . . . . . . 73.1 A Newtonian two-layer model . . . . . . . . . . . . . . . . . 73.2 A Bingham-Newtonian two-layer model . . . . . . . . . . . . 93.3 Final shape from plasticity theory . . . . . . . . . . . . . . . 133.4 Relations of slump and aspect ratio . . . . . . . . . . . . . . 13ivTable of Contents4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . 164.1 Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . 164.1.1 Regularization method for Navier-Stokes equations . 174.1.2 Augmented Lagrangian method for Navier-Stokes e-quations . . . . . . . . . . . . . . . . . . . . . . . . . 194.1.3 Algorithm for advection-diffusion equation . . . . . . 204.2 Parameter settings . . . . . . . . . . . . . . . . . . . . . . . . 224.3 Numerical result for Newtonian fluid problem . . . . . . . . 234.4 The no-slip boundary condition . . . . . . . . . . . . . . . . 274.5 Discussion about the result . . . . . . . . . . . . . . . . . . . 364.6 Numerical result for Bingham fluid problem . . . . . . . . . 375 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48AppendicesA Code Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . 51A.1 The results are independent of tolerance . . . . . . . . . . . 51A.2 The results are independent of time step . . . . . . . . . . . 51A.3 The domain size . . . . . . . . . . . . . . . . . . . . . . . . . 52A.4 The density and viscosity ratio . . . . . . . . . . . . . . . . . 53A.5 Reynolds number . . . . . . . . . . . . . . . . . . . . . . . . 53A.6 Regularization parameter is small enough . . . . . . . . . . . 53A.7 Problems with the slip boundary condition . . . . . . . . . . 54vList of Figures2.1 problem setting . . . . . . . . . . . . . . . . . . . . . . . . . . 43.1 final shape for B=0.02 to 0.3 by 1-order shallow layer model . 143.2 final shape for B=0.02 to 0.2 by slipline model . . . . . . . . 144.1 upwind scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 214.2 MUSCL scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 224.3 flow profile for T = 10, 40, 90, 160, 250, µ2µ1 =ρ2ρ1= 10−3, Re =9.81× 10−4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 244.4 x vs t for different resolution, µ2µ1 =ρ2ρ1= 10−3, Re = 9.81×10−4 254.5 h vs t for different resolution, µ2µ1 =ρ2ρ1= 10−3, Re = 9.81×10−4 254.6 concentration profile at T = 250, µ2µ1 =ρ2ρ1= 10−3, Re =9.81× 10−4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264.7 details about the concentration layer on the bottom . . . . . 264.8 horizontal velocity at x = 3.5, T = 250, there is evidentdiscontinuity in velocity field . . . . . . . . . . . . . . . . . . 274.9 numerical scheme around the boundary, stress in the grid cellis interpolated by velocity at the vertex . . . . . . . . . . . . 284.10 no-slip condition problem, T = 10, 40, 90, 160, 250, µ2µ1 =ρ2ρ1=10−1, Re = 9.81× 10−4 . . . . . . . . . . . . . . . . . . . . . 294.11 Navier slip boundary condition, T = 10, 40, 90, 160, 250, µ2µ1 =ρ2ρ1= 10−1, Re = 9.81× 10−4 . . . . . . . . . . . . . . . . . . 304.12 horizontal velocity at x=3.5, from Navier slip condition . . . 314.13 flow profile comparison of no-slip condition, Navier slip con-dition, modified slip condition, at the same time step, T =250, µ2µ1 =ρ2ρ1= 10−1, Re = 9.81× 10−4 . . . . . . . . . . . . 32viList of Figures4.14 details about the finger, T = 250, µ2µ1 =ρ2ρ1= 10−1, Re =9.81× 10−4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.15 horizontal velocity at x=3.5, from modified slip condition, thevelocity is continuous and converges well . . . . . . . . . . . . 334.16 x vs t for modified slip model, µ2µ1 =ρ2ρ1= 10−1, Re = 9.81×10−4 344.17 h vs t for modified slip model, µ2µ1 =ρ2ρ1= 10−1, Re = 9.81×10−4 344.18 flow length vs time for different models, µ2µ1 =ρ2ρ1= 10−1, Re =9.81× 10−4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.19 flow height vs time for different models, µ2µ1 =ρ2ρ1= 10−1, Re =9.81× 10−4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.20 modified slip boundary condition, T = 10, 40, 90, 160, 250, µ2µ1 =ρ2ρ1= 10−1, Re = 9.81× 10−4, white curves are from the two-layer shallow layer model . . . . . . . . . . . . . . . . . . . . 364.21 final shapes for different B . . . . . . . . . . . . . . . . . . . . 394.22 equilibrium flow length for B in [0.02 0.1] . . . . . . . . . . . 414.23 equilibrium flow height for B in [0.02 0.1] . . . . . . . . . . . 414.24 The color contour map is the stress field of the equilibriumstate for B=0.02, the black curve is the 1-order shallow layermodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.25 equilibrium flow length for B in [0.1 0.2] . . . . . . . . . . . . 434.26 equilibrium flow height for B in [0.1 0.2] . . . . . . . . . . . . 444.27 The color contour map is the stress field of the equilibriumstate for B=0.2, the white curve is the slipline result . . . . . 45A.1 x vs t for different tolerance . . . . . . . . . . . . . . . . . . . 51A.2 h vs t for different tolerance . . . . . . . . . . . . . . . . . . . 51A.3 x vs t for different time step . . . . . . . . . . . . . . . . . . . 52A.4 h vs t for different time step . . . . . . . . . . . . . . . . . . . 52A.5 x vs t for different domain size . . . . . . . . . . . . . . . . . 52A.6 h vs t for different domain size . . . . . . . . . . . . . . . . . 52A.7 x vs t for different viscosity ratio . . . . . . . . . . . . . . . . 53A.8 h vs t for different viscosity ratio . . . . . . . . . . . . . . . . 53A.9 x vs t for different Reynolds number . . . . . . . . . . . . . . 54viiList of FiguresA.10 h vs t for different Reynolds number . . . . . . . . . . . . . . 54A.11 flow length vs time for different BCs . . . . . . . . . . . . . . 55viiiAcknowledgementsI would like to convey my gratitude to all people who gave me the possi-bility to complete this thesis. In the first place, I would like to express mysincere gratitude to my supervisor, Professor Neil Balmforth for the contin-uous support of my MSc study and research, for his patience, motivation,enthusiasm, and immense knowledge. His wide knowledge and his logicalway of thinking have been of great value to me. I am deeply grateful toProfessor Sarah Hormozi for her supervision, advice, and guidance from theearly stage of this research. Above all and the most needed, they providedme support and friendly help in various ways. I would like to thank MrAnthony Wachs, for assisting me with the validation of the numerical codesand Professor James Feng for useful advice. Most importantly, I would liketo thank my entire family. My immediate family to whom this dissertation isdedicated, has been a constant source of love, concern, support and strengthall these years. I would like to express my heart-felt gratitude to my parentswho have supported and encouraged me throughout this endeavor. Finally,financial support of the Natural Sciences and Engineering Research Councilof Canada (NSERC) is gratefully acknowledged.ixDedicationTo my wonderful parents, who stood by me and supported all my ideas anddreams. I love you dearly. To my best friends in all aspects of my life, whohave been by my side throughout difficult time. Your constant support andencouragement continue to help me reach my goals. I am forever thankfulto have you in my life.xChapter 1Introduction1.1 MotivationA Bingham plastic is a viscoplastic material that behaves as a rigid bodyat low stress but flows as a viscous fluid at high stress. Unlike Newtonianfluid in which the viscous stress τ is linearly proportional to the strain rateγ, with a viscosity coefficient µ, the Bingham fluid satisfies the followingconstitutive law: γ = 0 τ < τyτ =(µ+ τyγ)γ τ > τywhere τy is a scalar called the yield stress, τ =√12τ : τ is the secondinvariant of the stress tensor τ , and γ =√12γ : γ is the second invariant ofthe strain rate tensor γ. When τ < τy, the strain rate is zero, meaning thatthere is no local deformation. When τ > τy, there is a linear relationshipthat τ − τy = µγ.Bingham type plasticity is commonly used as a mathematical model forcomplex fluid in industry, such as drilling mud, concrete rheology, and foodprocessing[16]. In those industries, the dynamic of gravity-driven flow playsan important role in a number of processes, including methods designed tomeasure fluid rheology. For example, the slump test is a way of measuringthe consistency of fresh concrete. Furthermore, natural phenomenon likeavalanche and debris flows can be modeled as a gravity flow of Binghamplastic[5]. A common feature of those phenomenon is that due to the yieldcriteria, the fluid will stop spreading as the stress fall below the yield stress,and form certain steady shape, instead of just continually flattening itself11.2. Literature reviewlike Newtonian fluid on scales long enough that surface tension is negligible.Therefore it is interesting to know what is the final shape of the fluid, andits relation with varying yield stress.1.2 Literature reviewFor Newtonian fluid, the viscous gravity current of a single fluid has beenwell studied by Huppert[7] and others (see Bankoff et al [1]). Numericaltests for low Reynolds number flow, however, have rarely been implementedto validate the lubrication theory.For Non-Newtonian fluid, a thin layer theory for Bingham-type fluid hasbeen developed by Liu and Mei[12], based on Reynolds’ lubrication theory.This theory leads to a prediction of the final shape in the shallow limit (seeBalmforth, Craster and Sassi[21]). For non-shallow slumps, slipline theoryfrom plasticity theory [2] has been applied to predict the final shape, byDubash et al[20].For numerical techniques of computing Bingham type fluid, the maindifficulty is how to deal with the stress below the yield stress τ , as it isundetermined in the constitutive law. One way of addressing this problem isby considering the fluid as highly viscous when it is unyielded, thus makingthe strain rate close to 0. This is done by doing a regularization to theconstitutive law as below:τ =(µ+ τyγ + )γwhere is a small parameter. If τ is sufficiently bigger than τy, γ , theoriginal constitutive equation is recovered. A second way of addressing theundetermined stress field is by introducing a Lagrangian multiplier tensorλ, such that{γ = λ : γτ = µγ + τyλeverywhere in the fluid. Theories prove the existence of such λ, includingthe region γ = 0, so that the stress can be recovered everywhere. We21.3. Objectivesrefer to Dean[6] for a review of those methods. However, since the originalconstitutive law does not give any information for the stress in the unyieldedregion, thus recovering of unyielded stress field in the numerical schemes isartificial.Numerical simulations of gravity current of Bingham type fluid are rare.Vola et al[4] developed a numerical scheme of computing gravity currents ofNon-Newtonian fluid, using a Lagrangian technique. Roussel and Coussot[19]used a software FLOW-3D to simulate the slump of a cylinder of Binghamfluid. Staron et al[15] did a simulation of columns of Bingham plastics undergravity, using regularization methods.For experiments about this problem, there are no experimental resultswith the same context as ours now. Published papers with related exper-iments include Pashias et al[23] (measurement of the yield stress by theslump test), Cochard and Ancey[24] (experimental results related to thedam-break problem for viscoplastic fluids which is a 3D case), and Balm-forth et al[18] (experimental analysis of the dam break of a viscoplastic fluidin a horizontal channel, with Herschel-Bulkley constitutive law).1.3 ObjectivesThree theoretical models are given from the past literature, predicting thefinal shape of the slump of Bingham fluid in terms of the yield stress andthe initial aspect ratio. Currently there are not many experiments to verifythose models. Therefore, my thesis aims to use numerical simulations tocompare with the theoretical models, in order to find the final shapes ofa 2D unit block of Bingham fluid under gravity. I am using a Complexfluid solver on the platform Pelicans[25], developed by IRSN, to simulatethe evolution.3Chapter 2Problem Description2.1 Problem settingWe study the slumping of a unit block of Bingham fluid under gravity, in a2D case. We will consider a two-phase flow in a rectangular domain.Figure 2.1 illustrates the domain and initial condition. The domain isset to be Ω = [0, L]× [0, H]. There are two fluids in the domain. Fluid 1 iswhat we are interested in, with a higher density ρ1 and viscosity µ1, whichwe also refer to as the lower layer. Fluid 2 is an ambient fluid, with muchlower density ρ2 and viscosity µ2. We refer to this fluid as the upper layer.Since we want the lower layer to be a Bingham fluid, it is necessary todefine the constitutive law, which isγ = 0 τ1(u) < τyτ1(u) =(µ1 +τyγ)γ(u) τ1(u) > τyfluid 1fluid 2ρ1, µ1ρ2, µ2LHRRFigure 2.1: problem setting42.2. EquationsAnd the upper layer is Newtonian, so thatτ2(u) = µ2γ(u)where τ1(u), τ2(u) are the stress tensor, γ(u) = ∇u+(∇u)T is the strain ratetensor, and τ1(u),γ are the second invariants of τ1(u), γ(u), respectively:τ1(u) =√12τ1(u) : τ1(u), γ(u) =√12γ(u) : γ(u)τy is the yield stress of fluid 1. That is, when the stress invariant exceedsτy, the flow will behave like fluid, otherwise it is like solid, with no localdeformation.2.2 EquationsSuppose both fluids are incompressible, the equations are defined as in e-quation (2.1),lower layerρ1[∂u∂t + (u · ∇)u]= −∇p+∇ · τ1(u) + ρ1g∇ · u = 0upper layerρ2[∂u∂t + (u · ∇)u]= −∇p+∇ · τ2(u) + ρ2g∇ · u = 0(2.1)where u = (v, w) is the velocity vector, p is the pressure, and g is thegravity. A no-slip boundary condition is applied on the bottom and top,where u = (0, 0). A symmetry boundary condition is set on the left andright wall, where v = 0, ∂w∂x = 0. The velocity and stress must also becontinuous at the interface.52.3. Dimensionless form2.3 Dimensionless formSuppose the characteristic length and velocity are R and U , let x = Rxˆ, z =Rzˆ, u = Uuˆ, t = RU tˆ, where U =ρ1gR2µ1. Here, uˆ, tˆ, pˆ, τˆ all stand for thedimensionless form. The equations are reduced tolower layerRe[∂uˆ∂tˆ+ (uˆ · ∇)uˆ]= −∇pˆ+∇ · τˆ1(uˆ) + 1¯∇ · uˆ = 0upper layerρ2ρ1Re[∂uˆ∂tˆ+ (uˆ · ∇)uˆ]= −∇pˆ+ µ2µ1∇ · τˆ2(uˆ) +ρ2ρ11¯∇ · uˆ = 0(2.2)and the constitutive law is reduced toτˆ1(uˆ) =(1 + Bγˆ)γˆτˆ2(uˆ) = γˆ(2.3)where Re = ρ21gR3µ21is the Reynolds number, and B = τyρ1gR is the Binghamnumber, which is the dimensionless yield stress. For convenience, in the restof this thesis, we just use u, p, τ, γ to represent the dimensionless form.6Chapter 3Asymptotic ResultsWe will start by considering a Newtonian fluid for the lower layer, and thenextend to a Bingham fluid.3.1 A Newtonian two-layer modelA dimensionless model for the propagation of two-dimensional Newtonianfilm under gravity is discussed by Huppert[7].∂h∂t =13∂∂x(h3∂h∂x)It is based on the assumption that there is only one layer of fluid and theshear stress on the surface is 0. In our problem, there is an ambient flu-id above the lower layer. If both are Newtonian fluids, the Navier-Stokesequation is written as below.lower layerux + wz = 0ρ1(ut + uux + wuz) = −px + µ1(uxx + uzz)ρ1(wt + uwx + wwz) = −pz + µ1(wxx + wzz) + ρ1gupper layerux + wz = 0ρ2(ut + uux + wuz) = −px + µ2(uxx + uzz)ρ2(wt + uwx + wwz) = −pz + µ2(wxx + wzz) + ρ2g(3.1)where (u,w) is the velocity vector. Now suppose the upper layer is notignored, and the viscosity ratio and density ratio are S = µ2µ1 , Q =ρ2ρ1. Usinga similar non-dimensionalization methods as in [22], with different length73.1. A Newtonian two-layer modelscale for the horizontal and vertical length, letx = Lxˆ, z = Rzˆ, U = ρ1gR3µ1L, v = (u,w) =(Uuˆ, RLUwˆ), t = LU tˆ, p = ρgRpˆ, Hˆ = H/Rand suppose RL = 1, Re =ρ1UR2µ1L 1. Equation (3.1) is nondimension-alized aslower layerux + wz = 0Re(ut + uux + wuz) = −px + 2uxx + uzzRe(wt + uwx + wwz) = − 12 (pz + 1) + 2wxx + wzzupper layerux + wz = 0QSRe(ut + uux + wuz) = − 1Spx + 2uxx + uzzQSRe(wt + uwx + wwz) = − 12S (pz +Q) + 2wxx + wzz(3.2)Getting rid of O() term, we have{px = uzz, pz = −1 for the lower layerpx = Suzz, pz = −Q for the upper layer(3.3)Let h = h(x, t) be the height of the lower layer. In order to solve equation3.3, we need to include the following conditions.• no-slip condition on the upper and lower boundary, u = 0 at z = 0, 1.• u, p, τ continuous at the interface, here τ ≈ µuz.• zero flux across the domain,∫ Hˆ0 udz = 0.83.2. A Bingham-Newtonian two-layer modelCombining all the relations above, we obtainu = 0 z = 0uzz = px, pz = −1 0 < z < hu, p, τ continuous z = hSuzz = px, pz = −Q h < z < Hˆu = 0 z = Hˆ∫ Hˆ0udz = 0ht =∂∂x(∫ h0udz)(3.4)Then we derive a two-layer model for a shallow gravity flow.∂h∂t =13(1−Q)∂∂x((k + Sh)(k2 − Sh2)2 + 4SkhHˆ2h3k3∂h∂x)where k = Hˆ − h. As S → 0 and Q → 0, the equation is reduced toHuppert’s model. Notice that this model is only valid for large aspect ratio,where the horizontal length scale is much more than the vertical, so thatτ can be approximated as µuz. Therefore, it may not be valid at the flowfront. This is to be shown in the next chapter, in comparison with numericalsimulations.3.2 A Bingham-Newtonian two-layer modelFor shallow layer model for Bingham fluid, Liu and Mei[12] presented thelubrication model for two dimensional flow.∂h∂t =16∂∂x(Y 2(3h− Y )∂h∂x)(3.5)93.2. A Bingham-Newtonian two-layer modelY is the height of the yielded region, which satisfiesY = max{h− B|hx|, 0}This model assumes that there is only one layer of Bingham fluid, and theshear stress on the surface is 0. Therefore, there is always an unyieldedregion on the top of the fluid.Therefore we consider the problem with a Newtonian fluid above a Bing-ham fluid, so that the stress on the surface is balanced by the upper fluid,and give the following analysis. Using a similar nondimensionalization pro-cedure as for the Newtonian two-layer model, we get the reduced form.{px = ∂zτ, pz = −1, 0 < z < hpx = Suzz, pz = −Q, h < z < Hˆ(3.6)as in equation (3.3), the only difference is that uzz is replaced by ∂zτ . Sincepz is piecewise linear, we can integrate to get{p = P − (z − h), τ = τ0 + (Px + hx)z 0 < z < hp = P −Q(z − h), px = Suzz h < z < Hˆ(3.7)For the lower layer 0 < z < h, define τ0 = τ(x, 0), τh = τ(x, h), equation(3.7) gives a relationτh = τ0 + (Px + hx)hLet us define uz as a function of τ , uz = Γ(τ), using the constitutive law Asτ is a linear function of z.zh =τ − τ0τh − τ0differentiating the equation above givesdz = hτh − τ0dτ103.2. A Bingham-Newtonian two-layer modelThe velocity at the interface uh isuh =∫ h0uzdz =hτh − τ0∫ τhτ0Γ(τ)dτFor convenience, define I0 =∫ τhτ0Γ(τ)dτ , so thatuh =hI0τh − τ0And the flux of the lower layer is integrated as∫ h0udz = zu|h0 −∫ h0zuzdz= huh −∫ τhτ0h(τ − τ0)τh − τ0Γ(τ)dz hτh − τ0If we define I1 =∫ τhτ0τΓ(τ)dτ , we can rewrite the flux integral so that∫ h0udz = h2(τhI0 − I1)(τh − τ0)2For the upper layer h < z < Hˆ, since Px +Qhx = Suzz at and Suz|h = τh,we haveuz =1S [(Px +Qhx)(z − h) + τh]Introducing Px = τh−τ0h − hx furnishesSuz =1S[(τh − τ0h − hx +Qhx)(z − h) + τh]We can get the interface velocity uh again by integrating uz in the upperlayeruh = −∫ Hˆhuzdz= − 1S (Hˆ − h)[τh +Hˆ − h2(τh − τ0h − hx +Qhx)]113.2. A Bingham-Newtonian two-layer modelAnd the flux of the upper layer is∫ Hˆhudz = zu|Hˆh −∫ Hˆhzuzdz= −huh −∫ HˆhzS[τh +(τh − τ0h − hx +Qhx)(z − h)]dz= − 12S (Hˆ − h)2τh −13S (Hˆ − h)3(τh − τ0h − hx +Qhx)(3.8)To make the velocity continuous at the interface, and make zero flux acrossthe region∫ h0 udz +∫ Hˆh udz = 0, we haveS hI0τh − τ0= −(Hˆ − h)[τh +Hˆ − h2(τh − τ0h − hx +Qhx)]h2(τhI0 − I1)(τh − τ0)2= 12S (Hˆ − h)2τh +13S (Hˆ − h)3(τh − τ0h − hx +Qhx)(3.9)The evolution equation ht + (∫ h0 udz)x = 0, becomes∂h∂t +∂∂x[ h2(τh − τ0)2(τhI0 − I1)]= 0For Bingham fluid, the relation of uz and τ isuz = Γ(τ) ={τ −B |τ | > B0 |τ | < BUsing equation (3.9), τh, τ0 can be solved by a Newton iteration method ineach time step, as indicated in [9]. If S → 0, Q→ 0, the equation is reducedto equation (3.5).This equation indicates a final state of the flow ash(x) =√2B(X − x)Where X is the flow length of the final shape. In our case, since the totalvolume is 1, we can calculate the value of X by demanding∫ X0 h(x)dx = 1123.3. Final shape from plasticity theoryand get the final flow length and height as a function of B:X = 12( 9B)1/3, H0 = h(0) = (3B)1/33.3 Final shape from plasticity theoryDubash et al [20] provided a higher order solution to the steady state of theBingham slump. For B 1h =√(1 + pi2B)2( 13B +pi4)−2/3− 2Bx− pi2B( 13B +pi4)−1/3+O(B3)This model assumes the flow height is normalized to be 1. In our settings,as the total volume is 1, the final height and length areH0 =( 13B +pi4)−1/3X =( 12B +pi2)( 13B +pi4)−2/3(3.10)Figure 3.1 illustrates the final shape of the fluid for varying Bingham num-ber, predicted by equation 3.10. Dubash [20] also used a slipline methodof plasticity theory to compute the final shape, based on the assumptionthat the flow is yielded everywhere before reaching the steady state, so thatstress falls to the yielded stress τy everywhere. There is no explicit form ofthe final shape in terms of B, algorithms about construction of the sliplinefield can be found in the reference. Figure 3.2 illustrates the final shape ofthe fluid for varying Bingham number, predicted by the slipline theory.3.4 Relations of slump and aspect ratioOther than shallow-layer analysis, Pashias et al[23] derived a relation ofthe slump height h and yield stress B, as well as the initial aspect ratio133.4. Relations of slump and aspect ratio0 0.5 1 1.5 2 2.5 3 3.5 4 4.500.10.20.30.40.50.60.70.80.91xzfinal shape for B=0.02 to 0.20, 1-order shallow layer modelincreasing BFigure 3.1: final shape for B=0.02 to 0.3 by 1-order shallow layer model0 0.5 1 1.5 2 2.5 3 3.5 4 4.500.10.20.30.40.50.60.70.80.91xzfinal shape of B=0.02 to 0.2 from slipline modelincreasing BFigure 3.2: final shape for B=0.02 to 0.2 by slipline model143.4. Relations of slump and aspect ratioH0/R0 = a.h = 1− 2B[1− ln(2Ba)]However this is based on the assumption that the aspect ratio is large andthe stress variations in the horizontal direction are negligible compared tothose in the vertical direction, which is different from our assumption. Nev-ertheless, it is still worth comparing with their model for relatively greateraspect ratio.Staron et al[15] did numerical simulations in the same context as ours,except that the aspect ratio is varying from 0.2 to 19(in our problem it iskept unity). They get a relations of the slump height and the yield stress.h = 3.01B0.66This relation is based on the observation of a group of numerical tests withvarying B and initial aspect ratio.15Chapter 4Numerical Results4.1 Numerical schemeWe use a Volume-of-Fluid method to simulate the two-phase fluid. Werefer to Hirt and Nichols[8] for a description of this method. In brief, aconcentration field c(x, y, t) is introduced to distinguish the two fluids andsmooth out the interface between them: c = 1 for the lower layer, andc = 0 for the upper layer. This field is taken to be a material invariant.The problem therefore boils down to the Navier-Stokes equation and theadvection-diffusion equation:Re(c)[∂u∂t + (u · ∇)u]= −∇p+∇ · τ(u) +Ri(c)g∇ · u = 0∂c∂t +∇ · (uc) = 0(4.1)whereRe(c) :=[c+ (1− c)ρ2ρ1] ρ1URµ1Ri(c) :=[c+ (1− c)ρ2ρ1] ρ1gR2µ1Uτ(u) :=[c+ (1− c)µ2µ1]γ(u)In every time step, the Navier-Stokes equation is solved using a Galerkinfinite element method. And then the advection diffusion equation is solvedusing a finite volume method, with a MUSCL scheme in the advection term.We refer to Mikhlin [26] for a description of Galerkin method and Van164.1. Numerical schemeLeer[14] for a description of MUSCL scheme.4.1.1 Regularization method for Navier-Stokes equationsFor a general Navier-Stokes equation, with no-slip condition in ΓD, no-fluxcondition in ΓN , as belowρ(∂u∂t + (u · ∇)u)= −∇p+∇ · τ(u) + ρg in Ω∇ · u = 0 in Ωu = 0 in ΓD∇u · n = 0 in ΓN(4.2)Here τ(u) is given by the regularized constitutive law.τ(u) =[c(1 + Bγ + ) + (1− c)µ2µ1]γ(u)where is the regularization parameter. On a certain time tn, given thecurrent velocity field un, time step ∆t, the discrete scheme for solving un+1isρ∆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(4.3)The variational form is: Find (u, p) ∈ Su×V p so that for ∀(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(4.4)174.1. Numerical schemeHere Su, V u, V p is defined asSu = {v ∈ H1(Ω)|v = 0 in ΓD}V u = {v ∈ H1(Ω)|v = 0 in ΓD}V p = L2(Ω)(4.5)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(4.6)wherem(u, v) =∫Ωρu · vc(u; v, w) =∫ρ(u · ∇)v · wk(u, v) =∫Ωτ(u) : ∇vb(v, q) = −∫Ωq∇ · vf(v) =∫Ωρg · v − 1∆tm(un, v)(4.7)This variational equation is then solved using a Galerkin method. On afinite dimensional subspace V un ⊂ 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(4.8)After a finite element basis is chosen,Xh = span{ψuk | 0 ≤ k < Nudof}Mh = span{ψpk | 0 ≤ k < Npdof}184.1. Numerical schemethe variational equation is transformed to the linear system,1∆tMU +MRU +AU +BTP = FBU = 0(4.9)whereMIJ = m(ψuJ , ψuI ) 0 ≤ I, J < NudofMRIJ = mR(ψuJ , ψuI ) (I, J) ∈ ΓRAIJ = c(unh;ψuJ , ψuI ) + k(ψuJ , ψuI ) 0 ≤ I, J < NudofBIJ = b(ψuJ , ψpI ) 0 ≤ I < Npdof , 0 ≤ J < NudofFI = f(ψuI ) 0 ≤ I, J < Nudof(4.10)4.1.2 Augmented Lagrangian method for Navier-StokesequationsUnlike the regularization method, the augmented Lagrangian method doesnot require the constitutive laws to be regularized, and can produce trulyunyielded regions where γ = 0. Given the Navier-Stokes equationρ (∂tu+ (u · ∇)u) = ∇ · τ + f ∈ Ω∇ · u = 0 ∈ Ωu = uB(t) on Γτ = −pI +√2B D(u)|D(u)| + 2µD(u)194.1. Numerical schemewhere D(u) = 12γ = 12 [∇u+(∇u)T ], Duvaut and Lions derive the variationalinequality model: Find {u(u), p(t)} ∈ (H1(Ω))d × L2(Ω) such thatρ∫Ω∂tu(t) · (v − u(t))dx+ ρ∫Ω(u(t) · ∇)u(t) · (v − u(t))dx+µ∫Ω∇u(t) : ∇(v − u(t))dx+√2B(j(v)− j(u(t)))−∫Ωp(t)∇ · (v − u(t))dx ≥∫Ωf(t) · (v − u(t))dx, ∀v ∈ VB(t)∇ · u(t) = 0 in Ωj(v) =∫Ω|D(v)|dx, ∀v ∈ (H1(Ω))dVB(t) = {v|v ∈ (H1(Ω))d, v = uB(t) on Γ}We refer to Dean [6] for a review of the algorithm. In our numerical simula-tion with this method, the convergence rate for the linear system solution istoo slow for an adequate spacial resolution, thus we fail to get the numericalresults with large deformation. Nevertheless, for larger Bingham numberwhere there is little deformation, numerical results are achieved.4.1.3 Algorithm for advection-diffusion equationFor any region Ω in the domain∂C∂t +∇ · (uC) = 0can be integrated as∂∂t∫ΩCdx+∫∂ΩC(u · n)ds = 0Suppose the region is chosen to be a rectangular mesh cell, with 4 sides.Then the equation on this cell can be discretized asA∆tCn+1 +∑i∈{n,e,s,w}Li(ui · ni)Cn+1i =A∆tCn204.1. Numerical schemeCCsCnCeCw LeLnLwLsLe = Lw, Ls = Ln, A = LeLsnsnennnwFigure 4.1: upwind schemewhere ui is the approximated velocity at the middle point of side i, withoutward normal ni and length Li, and Ci is C if ui · ni > 0, or the adjacentcell of C with common side i if ui · ni < 0. This is known as the upwindscheme. A is the cell area. The details are shown in Figure 4.1. This isthe basic first order upwind scheme. The MUSCL scheme includes a secondorder correction term on the RHS, using a Van Leer limiter [14]:φ(r) = r + |r|1 + |r|which achieves more accurate solutions for large gradients in C, as is thecase in our problem. In order to understand the use of the flux limiter, werefer to Kuzmin[13]. Figure 4.2 illustrates how this correction is computed,suppose u · n > 0, then the correction term is computed ascorrection = −ud · L(u · n)φ(crglhg)lhg214.2. Parameter settingsud ddCuCuu Cduudlhg = Cu−Cuuuud crg =Cd−Cuud+ddnuLcorrection = −ud · L(u · n)φ( crglhg )lhgFigure 4.2: MUSCL schemeIt is added to the RHS of the equation for Cu, and its opposite value isadded to the RHS of the equation for Cd. This scheme is implemented bothon horizontal and vertical sides.All of the algorithms above are implemented in C++ as an applicationof PELICANS, an object oriented platform developed at IRSN, France, toprovide general frameworks and software components for the implementa-tion of PDE solvers. PELICANS is distributed under the CeCILL licenseagreement. We refer to Hormozi[25] and Wielage-Burchard[11] for a detaileddescription of this code. 14.2 Parameter settingsThe domain size is chosen to be [0, L] × [0, H] = [0, 5] × [0, 1.25], and thelower layer is initialized to be [0, 1]× [0, 1] on the left wall. The viscosity anddensity are normalized to be 1 for the lower layer , and 10−3 for the upperlayer, so as to minimize the effect of the upper layer. The Reynolds numberis set to be 9.81 × 10−4, in order to minimize the effect of inertia force.We choose a uniform mesh on the domain, with ∆x = ∆y varying from1PELICANS is distributed under the CeCILL license agreement(http://www.cecill.info/licences/Licence CeCILL V2-en.html ). PELICANS can bedownloaded from https://gforge.irsn.fr/gf/project/pelicans/ .224.3. Numerical result for Newtonian fluid problem0.005 to 0.025. For the computational reason, the maximum number of gridpoints we can achieve is 1000× 250 = 250000. The time step for solving theNavier-Stokes equation is set to be ∆t = 0.5, and the time step for solvingthe advection-diffusion equation is restricted by the CFL condition( u∆x +w∆y)dt < 12where (u,w) is the velocity vector. The choice of all parameters are discussedin Appendix A.4.3 Numerical result for Newtonian fluid problemBefore solving the Bingham fluid problem, it is necessary to test for a New-tonian problem to see if the code works well.This problem is solved using the shallow layer model, and by numericalsimulation. The results are compared below. Figure 4.3 is the flow profileat different time steps, the white curve is interface produced by the two-layer model. Figure 4.4 and 4.5 show the flow length and height as afunction of time for different resolutions and the asymptotic model(Noticethat as viscosity and density ratio are close to 0, the two-layer model is veryclose to the Huppert’s model). According to the graph, the flow height isindependent of resolution, while the flow length is converging slowly. Figure4.6 illustrates a snapshot of the concentration field at time step T = 250.The white curve shows the interface, which is defined to be the level curveof c = 0.5. According to the graph, there is a ‘finger’ of upper fluid alongthe bottom underneath the lower fluid. Figure 4.7 shows the details of thefinger, which is evidently poorly resolved. As a consequence of this resolutionproblem, the horizontal velocity from the bottom to the first layer of gridshows clear discontinuity, see Figure 4.8 for an example. In this graph, x-axis is the vertical position of the domain, and the y-axis is the horizontalvelocity, plotted at x = 3.5, T = 250.We attempted to resolve the finger by refining the resolution. However,even in our most refined resolution where the grid size is ∆y = 0.001, the234.3. Numerical result for Newtonian fluid problem0 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.5Figure 4.3: flow profile for T = 10, 40, 90, 160, 250, µ2µ1 =ρ2ρ1= 10−3, Re =9.81× 10−4244.3. Numerical result for Newtonian fluid problem0 100 200 300 400 500 600 70011.522.533.544.55txx vs t for different resolution 200x50300x75400x100500x125600x150700x175800x200900x225asymptoticsFigure 4.4: x vs t for different resolution, µ2µ1 =ρ2ρ1= 10−3, Re = 9.81×10−40 100 200 300 400 500 600 7000.20.30.40.50.60.70.80.91thh vs t for different resolution 200x50300x75400x100500x125600x150700x175800x200900x225asymptoticsFigure 4.5: h vs t for different resolution, µ2µ1 =ρ2ρ1= 10−3, Re = 9.81×10−4254.3. Numerical result for Newtonian fluid problemxzvolume fraction c(x,y) at T=250 0 1 2 3 4 500.050.10.150.20.250.30.350.40.10.20.30.40.50.60.70.80.9Figure 4.6: concentration profile at T = 250, µ2µ1 =ρ2ρ1= 10−3, Re =9.81× 10−4xzvolume fraction c(x,y) at T=250 0 0.5 1 1.5 2 2.5 3 3.5 400.010.020.030.040.050.060.070.080.090.10.10.20.30.40.50.60.70.80.9Figure 4.7: details about the concentration layer on the bottom264.4. The no-slip boundary condition0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.100.511.522.533.544.5x 10-3zUU at x= 3.5 for different grid resolution at time 250 200x50300x75400x100500x125600x150700x175800x200900x225Figure 4.8: horizontal velocity at x = 3.5, T = 250, there is evident discon-tinuity in velocity field‘finger’ is still there and velocity is still effectively discontinuous. Figure4.9 shows a plot of the numerical scheme. To make the shear stress µ∂u∂zcontinuous around the boundary, we have[c1 + (1− c1)10−3] u2 − u1∆y ≈[c0 + (1− c0)10−3] u1 − u0∆yGetting rid of all small terms, and notice that the finger on the bottomimplies c0 ≈ 0.1, c1 ≈ 1, we will haveu2 − u1u1 − u0= u2 − u1∆y ×(u1 − u0∆y)−1≈ c0c1≈ 0.1This difference in velocity gradient explains why velocity field is not resolved.4.4 The no-slip boundary conditionThe finger is believed to be a direct result of imposing a no-slip boundarycondition. Since the boundary velocity is 0, the fluid near the bottom sticks274.4. The no-slip boundary conditionC0C1u0 = 0u1u2∆y∆yFigure 4.9: numerical scheme around the boundary, stress in the grid cell isinterpolated by velocity at the vertexto it. This is a common problem in dealing with the moving contact line.The same finger problem can be found in Hartel et al[2]. We may have alook at another case about the finger, which leads to more severe problem.Suppose the viscosity ratio µ2µ1 and density ratioρ2ρ1are all reset to be 0.1,which means the upper layer is heavier and more viscous than before, thenthere is an obvious finger on the bottom. See Figure 4.10, the flow profilefor different time steps for this case.In this case, a layer of the lighter fluid is underneath the heavier one.We do not know if this makes sense physically, as there are no related exper-iments. According to our analysis, if the finger is physically true, then theonly way to resolve it is to refine the resolution along the bottom. However,due to our computation ability, we are not able to resolve the finger. If thefinger is a numerical artefact, then we should think of ways to remove it,so that there is no resolution issue caused by the finger. In this situation,we can only deal with the case that the finger is not physically true, andwe choose to remove the finger by introducing a slip condition on the bot-tom. From a physical viewpoint, there are several studies showing that thereexists slip at some micro-scale surface, see Chang-Hwan Choi et al[3] for ex-ample. The slip length is of a µm scale, which is far from the macroscopicdescription of our problem. The slip condition we are trying to implement284.4. The no-slip boundary condition0 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.5Figure 4.10: no-slip condition problem, T = 10, 40, 90, 160, 250, µ2µ1 =ρ2ρ1=10−1, Re = 9.81× 10−4is based on a numerical artefact, not physical reason. This is commonlyused in removing the singularity in force. Liu et al[10], for example, solvesfor the spreading of droplet with a slip boundary condition. We also referto Renardy et al[17] for a discussion about the slip boundary condition inVOF methods. As there is slip velocity along the bottom, the finger can beremoved as the flux at the boundary is evidently increased. We first applyu = λ∂u∂ywhich is known as the Navier slip condition. λ is usually chosen to be on theorder of the mesh size ∆y. So we choose λ = ∆y, and solve for the problemwith µ2µ1 =ρ2ρ1= 0.1. Figure 4.11 is the flow profile at T = 10, 40, 90, 160, 250.Compared with the no-slip condition result 4.10, the finger is partly removedand the flow profile totally changes, which we think is reasonable, based onthe assumption that the finger is numerical artefact. We also modified the294.4. The no-slip boundary condition0 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.5Figure 4.11: Navier slip boundary condition, T = 10, 40, 90, 160, 250, µ2µ1 =ρ2ρ1= 10−1, Re = 9.81× 10−4Navier slip condition to provide another model:u = 1− cc ∆y∂u∂ywhere c = c(x, y) is the volume fraction, this is because the Navier slipcondition does not fully remove the finger and there is still discontinuityin velocity field, see Figure 4.12, for example. And we add a term 1−cc ,considering that, for c = 1 we have u = 0, meaning that the lower layersatisfies the no-slip condition, and for c = 0 we have ∂u∂y = 0, meaning thatthe upper layer satisfies the free slip condition, while for 0 < c < 1 thetwo conditions are mixed. In this way, the finger is effectively removed,see Figure 4.13 and 4.14 for example, in those graphs the flow profile atthe same time step for three different boundary conditions are compared,among which the finger is removed mostly by the modified slip model.Furthermore, the velocity field is continuous for the slip model, see Figure4.15 as an example. Therefore, we choose the modified slip model304.4. The no-slip boundary conditionFigure 4.12: horizontal velocity at x=3.5, from Navier slip condition314.4. The no-slip boundary condition0 0.5 1 1.5 2 2.5 3 3.5 400.20.40 0.5 1 1.5 2 2.5 3 3.5 400.20.40 0.5 1 1.5 2 2.5 3 3.5 400.20.4Figure 4.13: flow profile comparison of no-slip condition, Navier slip condi-tion, modified slip condition, at the same time step, T = 250, µ2µ1 =ρ2ρ1=10−1, Re = 9.81× 10−4 0 0.5 1 1.5 2 2.5 3 3.5 400.050.100.51 0 0.5 1 1.5 2 2.5 3 3.5 400.050.100.51 0 0.5 1 1.5 2 2.5 3 3.5 400.050.100.51Figure 4.14: details about the finger, T = 250, µ2µ1 =ρ2ρ1= 10−1, Re =9.81× 10−4324.4. The no-slip boundary condition0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.100.511.522.533.544.5 x 10-3zUU at x= 3.5 for different grid resolution at time 250 400x100500x125600x150700x175800x200900x2251000x250Figure 4.15: horizontal velocity at x=3.5, from modified slip condition, thevelocity is continuous and converges wellu = 1− cc ∆y∂u∂yas an alternate for the no-slip condition, in order to remove the finger andavoid the resolution issue. Figure 4.17 and 4.16 show better convergenceof flow height and length versus time, than the no-slip condition. Thus wecan assume that this model is independent of resolution, even though theboundary condition depends on the mesh size and the boundary volumefraction. Then we compare the numerical results with the asymptoticmodels, Figure 4.18 and 4.19 give the flow length and height vs time inthe case µ2µ1 =ρ2ρ1= 0.1. In principle, it is not reasonable to compare anumerical simulation with slip condition, with an asymptotic model, withno-slip condition. However, as the slip model depends on the mesh size,we cannot impose it into asymptotic analysis which does nothing aboutmesh. However, we want to argue that the slip model effectively reduces thepropagation speed from the no-slip condition. That is, with a slip condition,the flow is expected to spread faster than no-slip condition, however, it is the334.4. The no-slip boundary condition0 100 200 300 400 500 600 70011.522.533.544.55txx vs t for different resolution 400x100500x125600x150700x175800x200900x225asymptoticsFigure 4.16: x vs t for modified slip model, µ2µ1 =ρ2ρ1= 10−1, Re = 9.81×10−40 100 200 300 400 500 600 7000.20.30.40.50.60.70.80.91thh vs t for different resolution 400x100500x125600x150700x175800x200900x225asymptoticsFigure 4.17: h vs t for modified slip model, µ2µ1 =ρ2ρ1= 10−1, Re = 9.81 ×10−4344.4. The no-slip boundary condition0 100 200 300 400 50011.522.533.544.55 flow length vs timetx Huppert modeltwo-layer modelno-slip BCour slip BCNavier-slip BCFigure 4.18: flow length vs time for different models, µ2µ1 =ρ2ρ1= 10−1, Re =9.81× 10−40 100 200 300 400 5000.40.50.60.70.80.91 flow height vs timeth Huppert modeltwo-layer modelno-slip BCour slip BCNavier-slip BCFigure 4.19: flow height vs time for different models, µ2µ1 =ρ2ρ1= 10−1, Re =9.81× 10−4354.5. Discussion about the result0 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.50 1 2 3 4 500.5Figure 4.20: modified slip boundary condition, T = 10, 40, 90, 160, 250, µ2µ1 =ρ2ρ1= 10−1, Re = 9.81 × 10−4, white curves are from the two-layer shallowlayer modelopposite. This might indicate that: in numerical simulation, imposing a slipcondition might be more effective than a no-slip condition, in solving for amoving contact line problem. Figure 4.20 is the flow profile for different timesteps, from the modified slip model, compared with the two-layer shallowlayer theory. If the shallow layer theory is correct, that implies that thenumerical results with a slip condition is similar to the asymptotic modelwith no-slip condition. In this case, we can think of the slip condition as anumerical artefact, aimed at recovering the resolution issue caused by theno-slip condition.4.5 Discussion about the resultIn our numerical simulation, we first impose a no-slip boundary condition onthe bottom of the domain. However, the computation is not resolved, due tothe existence of a finger of the upper fluid along the bottom. It is possible364.6. Numerical result for Bingham fluid problemthat a highly refined resolution may help resolve the finger, however, ourmost refined resolution does not give a well-resolved result. Therefore, wechange the boundary condition to be slip, in order to remove the finger. Inour modified slip model u = 1−cc ∆y ∂u∂y , it is still a no-slip condition as c = 1,but a free-slip condition as c = 0. This implies that for the lower layer, theboundary condition is still no-slip, but for the upper layer, the boundarycondition is free-slip. Since we are not interested in the dynamics of theupper layer, imposing a free-slip condition will not cause severe problem.Then we make a resolution study of the modified slip model, the result isencouraging, because the finger is removed, and the velocity field near theboundary is resolved. Moreover, for the same problem considered, the slipmodel spreads more slowly than the no-slip model, meaning that numerical-ly, the slip model is better to simulate the problem than the no-slip model.All these results are based on the assumption that the finger is a numericaleffect. However, if the finger exists in reality, we have to have much refinedresolution to resolve it. Also it is possible that the thickness of the finger isso small that its nature is beyond the reach of macro-scale fluid dynamics.Anyway, as the modified slip model is the best one so far in simulating ourproblem, we will impose it in the problem for Bingham fluid case.4.6 Numerical result for Bingham fluid problemUnlike Newtonian fluid, Bingham fluid will stop at a final state. In theprevious chapter, we gave several theoretical models predicting the lengthand height of the final state, as well as the equilibrium flow profile. Now wewill use numerical simulation to get the final shape, for varying yield stress.The numerical settings for the standard problem areρ2ρ1= 10−3, µ2µ1= 10−3, Re = 9.81× 10−4for varying Bingham number B from 0.02 to 0.6. The resolution is set to∆x = ∆y = 0.005. Numerically, the regularized model τ = (1 + Bγ+)γassumes that the viscosity is O(B ) for small value of γ, thus the fluid will374.6. Numerical result for Bingham fluid problemnever stop spreading. Therefore, it is necessary to define a stopping criteria.Remember in the numerical scheme, the constitutive law isτ(u) =[c(1 + Bγ + )+ (1− c)µ2µ1]γ(u)Getting rid of small terms when γ 1, and take the second invariant,the fluid reaches the equilibrium as τ ≈ Bc. Therefore in our problem thecomputation is considered as stopped when τ < 1.01Bc everywhere in thefluid.Figure 4.21 is the numerical final shape for B varying from 0.02 to 0.3,with slip and no-slip condition. For larger B, the fluid is unyielded fromthe very beginning, as the stress everywhere is below the yield stress. AsB < 0.28, the fluid starts to spread, with a horn-like unyielded region on theupper right corner, which is not predicted by any asymptotic theories. Andas B goes even lower, the unyielded region is less evident and flattened, andthe flow is yielded everywhere, and it approaches certain final shape. Theresults of the no-slip condition has a finger on the bottom, the fluid spreadsfarther than the slip condition.We will compare the numerical results with the following models frompast literature:• The 0-order shallow layer modelX = 12( 9B)1/3, H = (3B)1/3This assumes that B 1 so that the aspect ratio of the fluid is small.• The 1-order shallow layer modelX =( 12B +pi2)( 13B +pi4)−2/3, H =( 13B +pi4)−1/3from Dubash et al[20]. This again assumes B 1.• slipline theory model from Dubash et al[20], there is no explicit form of384.6. Numerical result for Bingham fluid problemxzfinal shape for B from 0.02 to 0.30, slip condition0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 500.20.40.60.811.2(a) slip conditionxzfinal shape for B from 0.02 to 0.30, no slip condition0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 500.20.40.60.811.2(b) no slip conditionFigure 4.21: final shapes for different B394.6. Numerical result for Bingham fluid problemthe equation for the interface, the algorithm of recovering the interfaceis in this paper. It does not require B to be small, but assume thefluid yields everywhere at the beginning.• Pashias et al[23] gives a relation of the equilibrium flow height andBingham numberH = 1− 2B[(1− ln(2B)]It is based on the assumption of a much larger aspect ratio, while inour problem it is normalized to be 1. We include their result just tosee how different the results are under different assumptions.• Staron et al[15] gives a relation of the equilibrium flow height andBingham numberH = 3.01B0.66They have generalized the relations of the equilibrium length and theBingham number, but it is also dependent on viscosity in that model.As our methods of nondimensionalization are different, we do not havethe viscosity available to compare, thus we only use the equilibriumheight model.For B in [0.02 0.1], the equilibrium spreading length and height areshown in Figure 4.22 and 4.23. Here we also include the numerical result-s with the no-slip condition to see the problem caused by the finger. In4.22, the numerical no-slip results are still not resolved, leading to a furtherpropagation than the numerical slip results and other theoretical models,for B < 0.05. Thus the numerical no-slip results are not believable, andthe slip model makes the computation resolved, and it fits better with theasymptotic models than the no-slip results. As is said before, the modifiedslip model is the best one so far in simulating our problem, this comparisonfurther confirms our expectation about the application of the slip model.We can also see from the graphs that the numerical results agree well with1-order shallow layer theory and the slipline theory for smaller B, this indi-cates that the shallow layer model and the slipline model from Dubash [20]404.6. Numerical result for Bingham fluid problem0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.0911.522.533.544.55BXequilibrium flow length for different B numerical no-slipnumerical slipslipline0-order shallow layer1-order shallow layerFigure 4.22: equilibrium flow length for B in [0.02 0.1]0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.0900.10.20.30.40.50.6BHequilibrium flow height for different B numerical no-slipnumerical slipslipline0-order shallow layer1-order shallow layerStaron et alPashias et alFigure 4.23: equilibrium flow height for B in [0.02 0.1]414.6. Numerical result for Bingham fluid problem0 1 2 3 4 500.050.10.150.20.250.30.350.40.450.5xzsecond invariant of stress tensor at B=0.02 equilibrium state 0.0020.0040.0060.0080.010.0120.0140.0160.0180.020.022Figure 4.24: The color contour map is the stress field of the equilibriumstate for B=0.02, the black curve is the 1-order shallow layer modelare valid for B 1. Figure 4.24 is the equilibrium state for B = 0.02, theblack curve is the 1-order shallow layer model. They correspond well witheach other. We can also see from the stress field that the fluid yields every-where and the stress falls to the dimensionless yield stress B as it reachesthe equilibrium. From Figure 4.23, both Staron and Pashias’ models are notin a good comparison with the numerical results or the lubrication theory.Pashias’ model is different mainly because they assume the equilibrium statehas a larger aspect ratio, while in our problem it is small. Staron’s model isdifferent because they assume that the slump height is linearly dependenton the initial aspect ratio, for a constant B, and derive that the equilibriumflow height is on the order of O(B2/3), while this is not true in the shallowlayer model(O(B1/3)), which is more close to the numerical results. For Bin [0.1 0.2], the equilibrium spreading length and height are shown in Figure4.25 and 4.26. There is no big difference between the no-slip model and theslip model. And also the numerical results deviate from both the shallowlayer model and other models. This is because B is not small enough, and424.6. Numerical result for Bingham fluid problem0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.211.522.533.5BXequilibrium flow length for different B numerical no-slipnumerical slipslipline0-order shallow layer1-order shallow layerFigure 4.25: equilibrium flow length for B in [0.1 0.2]some region is not yielded. For example, Figure 4.27 is the stress field ofB = 0.2 at equilibrium state, the white curve is the final shape predicted bythe slipline model, we can see clearly that, at the right-upper corner, thereis an unyielded region, which makes the flow spread less than expected.434.6. Numerical result for Bingham fluid problem0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.200.20.40.60.81BHequilibrium flow height for different B numerical no-slipnumerical slipslipline0-order shallow layer1-order shallow layerStaron et alPashias et alFigure 4.26: equilibrium flow height for B in [0.1 0.2]444.6. Numerical result for Bingham fluid problemτ of B=0.2 equilibrium state vs slipline modelxz 0 0.5 1 1.5 2 2.500.20.40.60.810.050.10.150.2Figure 4.27: The color contour map is the stress field of the equilibriumstate for B=0.2, the white curve is the slipline result45Chapter 5ConclusionWe have performed numerical simulations for a two-layer fluid problem,using a VOF method for the evolution of the fluid and a regularizationmethod for the constitutive law of Bingham fluid, in order to simulate thepropagation of the Bingham fluid gravity current.A Newtonian fluid problem is solved first to validate the code. A res-olution problem caused by the no-slip boundary condition is found, due tothe appearance of a thin finger of upper layer fluid that is over-ridden bythe lower fluid. We are not sure if the finger is physically reasonable ornot. If it really exists, it probably requires much higher resolution to berecovered, or demands more than macro-scale fluid dynamics to be under-stood. Therefore, we consider the possibility that the finger is a numericaleffect from the no-slip condition and needs to be removed. A slip boundarycondition is shown to remove the finger formed by the upper fluid, therebyavoiding numerical resolution issues. However, we are not convinced thatthe numerical simulation with a slip boundary condition, is still our originalproblem(with no-slip condition). So far as we achieved, the slip model isproved to make the result independent of resolution, and it performs betterthan the no-slip condition results, in comparison with the asymptotic mod-els. Those observations, may imply that the slip condition is an alternativeof no-slip condition in numerical simulation of two-layer fluid moving on thewall. Further investigation is still needed.The slip boundary condition is then used for the Bingham fluid prob-lem. The numerical results for the Bingham fluid problem agree with theasymptotic models for B < 0.05, for B > 0.1, they become different, asthere is genuine unyielded region on the upper-right corner of the fluid. andthe final shapes are totally different. As B > 0.28, the fluid is always static.46Chapter 5. ConclusionTherefore, we numerically confirm the validity of the shallow layer theoryin certain parameter regimes.Although the numerical simulations agree with the asymptotic results,we are not very confident with the numerical result because the mechanismbehind imposing a slip condition is not yet understood, and we are notsure if the modified slip condition u = 1−cc ∆y ∂u∂y is a proper choice(Thereis a problem caused by imposing the slip-condition in the Bingham fluidproblem. Details can be found in Appendix.A.). In fact, it is a commonproblem that no-slip condition poses a problem in viscous flow theory atcontact lines: places where an interface between two fluids meets a solidboundary. Here, the no-slip boundary condition implies that the position ofthe contact line does not move, which is not observed in reality. Analysis ofa moving contact line with the no slip condition demands infinite stresses.Our next step includes: extending our simulation to an axisymmetriccontext, and varying the initial shape and aspect ratio and application of theaugmented Lagrangian method. We hope to make use of both experimentsand large number of numerical simulations to address the finger problem. Ifthe finger is not observed in experiments, that means it is simply a numericalartefact, then the slip model can be validated by comparison with no-slipasymptotic model and experiments, in different contexts.47Bibliography[1] S.G.Bankoff A.Oron, S.H.Davis. Long scale evolution of thin liquidfilms. Reviews of Modern Physics, 69(3):931, 1997.[2] F. Necker C. Hartel, E. Meiburg. Analysis and direct numerical simu-lation of the flow at a gravity-current head. part 1. flow topology andfront speed for slip and no-slip boundaries. J. Fluid Mech., 2000.[3] Kenneth S. Breuer Chang-Hwan Choi, K. Johan A. Westin. Apparentslip flows in hydrophilic and hydrophobic microchannels. Physics ofFluids, 15(10), 2003.[4] J.C. Latche D. Vola, L. Boscardin. Laminar unsteady flows of binghamfluids: a numerical strategy and some benchmark results. Journal ofComputational Physics, 2003.[5] L. Preziosi E. Bovet, B. Chiaia. A new model for snow avalanche dy-namics based on bingham fluids. Meccanica, 45:753–765, 2010.[6] Giovanna Guidoboni Edward J. Dean, Roland Glowinski. On the nu-merical simulation of bingham visco-plastic flow: Old and new results.J. Non-Newtonian Fluid Mech, 142:36–62, 2007.[7] H.E.Huppert. The propagation of two-dimensional and axisymmetricviscous gravity currents over a rigid horizontal surface. J. Fluid Mech.,121:43–58, 1982.[8] C. W. Hirt and B. D. Nichols. Volume of fluid (vof) method forthe dynamics of free boundaries. JOURNAL OF COMPUTATION-AL PHYSICS, 39:201–225, 1981.48Bibliography[9] N.J.Balmforth I.J.Hewitt. Viscoplastic lubrication theory with appli-cation to bearings and the washboard instability of a planing plate. J.Non-Newtonian Fluid Mech., 2012.[10] Y.F.Yap J.Liu, N.T.Nguyen. Numerical studies of sessile droplet shapewith moving contact lines. Micro and Nanosystems, 2011.[11] I.A. Frigaard K. Wielage-Burchard. Static wall layers in plane channeldisplacement flows. Journal of Non-Newtonian Fluid Mechanics, 2011.[12] C.C.Mei K.F.Liu. Approximate equations for the slow spreading of athin sheet of bingham plastic fluid. Physics of Fluids A: Fluid Dynam-ics, 2:30, 1990.[13] D. Kuzmin. On the design of general-purpose flux limiters for finite el-ement schemes. i. scalar convection. Journal of Computational Physics,219:513–531, 2006.[14] B.van Leer. Towards the ultimate conservative difference scheme. v.asecond-order sequel to godunov’s method. Journal of ComputationalPhysics, 142(1):101–136, 1979.[15] P.Ray S.Popinet L.Staron, P.Lagree. Scaling laws for the slumping ofa bingham plastic fluid. J. Rheol., 2013.[16] Evan Mitsoulis. Flows of viscoplastic materials: models and computa-tions. Rheology Reviews, pages 135–178, 2007.[17] J.Li M.Renardy, Y.Renardy. Numerical simulation of moving contactline problems using a volume-of-fluid method. Journal of Computation-al Physics, 171:243–263, 2001.[18] P. Perona A. C. Rust R. Sassi N. J. Balmforth, R. V. Craster. Vis-coplastic dam breaks and the bostwick consistometer. Journal of Non-Newtonian Fluid Mechanics, 142:63–78, 2007.49[19] P. Coussot N. Roussel. fifty-cent rheometer for yield stress measure-ments: From slump to spreading flow. Journal of Rheology, 49:705,2005.[20] A.C.Slim S.Cochard N.Dubash, N.J.Balmforth. What is the final shapeof a viscoplastic slump? J. Non-Newtonian Fluid Mech., 158:91–100,2009.[21] R. Sassi N.J. Balmforth, R.V. Craster. Shallow viscoplastic flow on aninclined plane. J. Fluid Mech., 470:1–29, 2002.[22] A.C.Rust R.Sassi N.J.Balmforth, R.V.Craster. Viscoplastic flow overan inclined surface. J.Non-Newtonian Fluid Mech., 142:219–243, 2007.[23] J.Summers D.J.Glenister N.Pashias, D.V.Boger. A fifty cent rheometerfor yield stress measurement. Journal of Rheology, 1996.[24] C. Ancey S. Cocharda. Experimental investigation of the spreading ofviscoplastic fluids on inclined planes. Journal of Non-Newtonian FluidMechanics, 158:73–84, 2009.[25] I. A. Frigaard S. Hormozi, K.Wielage-Burchard. Entry and start upeffects in visco-plastically lubricated viscous shear flow in pipe. Journalof Fluid Mechanics, 673:432–467, 2011.[26] S.G.Mikhlin. Variational methods in Mathematical Physics. PergamonPress, 1964.50Appendix ACode ValidationA.1 The results are independent of toleranceThe convergence criteria for solving the Navier-Stokes equations is whentolerance tol > max(||∇ · u||, ||un+1 − un||). We test the tol = 10−6, 10−8to see if the results are different, Figure A.1 and A.2 are the plots for theevolution of flow height and length, for the two tolerance(with all othersettings the same). It shows that tol = 10−6 is enough for convergence.A.2 The results are independent of time stepIt is necessary to show that the time step ∆T for each iterations of Navier-Stokes solver is small enough so that the velocity during any time is wellcaptured. The time step here is set to be 0.5, 0.25, 0.125, with all othersettings the same. Figure A.3 and A.4 are the results. It shows that ∆T =0.5 is small enough.0 50 100 150 200 25011.522.533.54 flow length vs timetx tol=1e-6tol=1e-8Figure A.1: x vs t for different toler-ance0 50 100 150 200 2500.40.50.60.70.80.91 flow height vs timeth tol=1e-6tol=1e-8Figure A.2: h vs t for different toler-ance51A.3. The domain size0 50 100 150 200 25011.522.533.54 flow length vs timetx dt=0.5dt=0.25dt=0.125Figure A.3: x vs t for different timestep0 50 100 150 200 2500.40.50.60.70.80.91 flow height vs timeth dt=0.5dt=0.25dt=0.125Figure A.4: h vs t for different timestep0 50 100 150 200 25011.522.533.54 flow length vs timetx LxH=5x1.25LxH=6x1.5Figure A.5: x vs t for different domainsize0 50 100 150 200 2500.40.50.60.70.80.91 flow height vs timeth LxH=5x1.25LxH=6x1.5Figure A.6: h vs t for different domainsizeA.3 The domain sizeSince we want to make the effect of the upper layer as small as possible, itis necessary to make sure that the change in domain size will not lead toa big difference in the result. We choose L ×H = 5 × 1.25, 6 × 1.5, whilethe initial size of the lower layer is still 1 × 1. Figure A.5 and A.6 are theresults. It shows that the effect of the domain size is small enough, thus wechoose L×H = 5× 1.25 as the standard setting in our problem.52A.4. The density and viscosity ratio0 50 100 150 200 25011.522.533.54 flow length vs timetx µ2/µ1=1e-3µ2/µ1=5e-4Figure A.7: x vs t for different viscos-ity ratio0 50 100 150 200 2500.40.50.60.70.80.91 flow height vs timeth µ2/µ1=1e-3µ2/µ1=5e-4Figure A.8: h vs t for different viscos-ity ratioA.4 The density and viscosity ratioFor the same reason as the domain size test, it is necessary to show thatthe change in density and viscosity ratio will not leading to big difference inthe results. We make the viscosity ratio µ2µ1 = 10−3, 5× 10−4, with all othersettings the same. Figure A.7 and A.8 show the results. It implies thatµ2µ1= 10−3 is small enough to make the results independent of the upperlayer.A.5 Reynolds numberFor Bingham fluid, the inertial force sometimes can change the equilibriumshape. To avoid the effect of inertia, it is necessary to make sure that theReynolds number is small enough, so that the Bingham fluid propagation isnot affected by inertia. We do computations for Re = 9.81 × 10−4, 7.85 ×10−3, Figure A.9 and A.10 are the results. It shows that the inertia can beneglected.A.6 Regularization parameter is small enoughThe regularization parameter of the model τ = (1 + Bγ+)γ needs to besmall enough so that it has no effect on the propagation of Bingham fluid53A.7. Problems with the slip boundary condition0 50 100 150 200 25011.522.533.54tXflow length vs time for diffrent Reynolds number Re=9.81e-4Re=7.85e-3Figure A.9: x vs t for differentReynolds number0 50 100 150 200 2500.40.50.60.70.80.91tHflow length vs time for different Reynolds number Re=9.81e-4Re=7.85e-3Figure A.10: h vs t for differentReynolds numberbefore stopping. We did a test by setting B = 0.04, and = 10−4 to 10−8,with all other settings the same. It indicates that as ≤ 10−6, there is nodifference in the results. To make the effect small enough, is always chosento be 10−8 in our tests.A.7 Problems with the slip boundary conditionFigure A.11 compares the flow evolution of the two types of boundary con-ditions for B = 0.1. We can see that, before reaching the equilibrium, theslip model spreads more slowly than the no-slip one. However, around theequilibrium, the no-slip model stops right away, while the slip model contin-ues to spread. The observation above implies that the modified slip modelis invalid around the equilibrium. It needs to be corrected to address thisissue. Therefore, the no-slip model is a relatively more valid model for Bin [0.1 0.2], as there is not much of a finger. This observation reflects thedeficiency of slip condition in Bingham fluid problem, as we cannot tell ifthe flow really stops at the equilibrium, which implies that the slip modelis not a proper choice in this Bingham fluid problem. One way of resolvingthis is to use the stopping time from the no-slip results.54A.7. Problems with the slip boundary condition0 100 200 300 400 500 600 700 800 900 10000.811.21.41.61.822.2tXflow length vs time for B=0.1 numerical slipnumerical no-slipFigure A.11: flow length vs time for different BCs55
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The propagation of the gravity current of Bingham fluid...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The propagation of the gravity current of Bingham fluid : a numerical investigation Liu, Ye 2014
pdf
Page Metadata
Item Metadata
Title | The propagation of the gravity current of Bingham fluid : a numerical investigation |
Creator |
Liu, Ye |
Publisher | University of British Columbia |
Date Issued | 2014 |
Description | We have studied the propagation of 2D unit block of viscoplastic fluid of Bingham type over a horizontal plane, underneath another Newtonian fluid. We numerically simulate the dynamics of a two-layer fluid in a rectangle domain, using the volume-of-fluid method to deal with the evolution of the interface, and regularization scheme of the constitutive law, which replaces unyielded plugs with very viscous flow. We explore the final shape of the flow for varying yield stress, comparing the numerical results with the predictions of the asymptotic theory, a plasticity model based on slipline theory, and other past results. Numerical difficulties with the moving contact lines are encountered during the numerical simulation. A slip boundary condition is used to address this issue, the validity of which should be further investigated. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2014-12-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0167092 |
URI | http://hdl.handle.net/2429/51637 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2015-02 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2015_february_liu_ye.pdf [ 2.64MB ]
- Metadata
- JSON: 24-1.0167092.json
- JSON-LD: 24-1.0167092-ld.json
- RDF/XML (Pretty): 24-1.0167092-rdf.xml
- RDF/JSON: 24-1.0167092-rdf.json
- Turtle: 24-1.0167092-turtle.txt
- N-Triples: 24-1.0167092-rdf-ntriples.txt
- Original Record: 24-1.0167092-source.json
- Full Text
- 24-1.0167092-fulltext.txt
- Citation
- 24-1.0167092.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0167092/manifest