We consider the flow of a viscoplastic fluid on a horizontal or an inclined surface with a flat and
an asymmetric topography. A particular application of interest is the spread of a fixed mass – a block
– of material under its own weight. The rheology of the fluid is described by the Bingham model
which includes the effect of yield stress, i.e. a threshold stress which must be exceeded before flow
can occur. Both the plastic viscosity and the yield stress are modelled with temperature-dependent
parameters. The flow is described by the lubrication approximation, and the heat transfer by a depthaveraged energy conservation equation. Results show that for large values of the yield stress, only the
outer fraction of the fluid spreads outward, the inner fraction remaining unyielded. We also present
an analysis which predicts the threshold value of the yield stress for which partial slump occurs

Cusack D

Hewett J

Kennedy B

Moyers Gonzalez, Miguel

Sellier M

UC Research Repository

23rd Australasian Fluid Mechanics Conference - 23AFMCSydney, Australia Paper No: AFMC2022-2174 – 8 December 2022Thin-film flow of a Bingham fluid over topography with a temperaturedependent rheologyMiguel Moyers-González1∗, James Hewett2, Dale Cusack3, Ben Kennedy3 and Mathieu Sellier21School of Mathematics and Statistics, University of Canterbury, NZ. 2Department of Mechanical Engineering,University of Canterbury, NZ. 3School of Earth and Environment, University of Canterbury, NZ.∗mailto: miguel.moyersgonzalez@canterbury.ac.nzAbstractWe consider the flow of a viscoplastic fluid on a horizontal or an inclined surface with a flat andan asymmetric topography. A particular application of interest is the spread of a fixed mass – a block– of material under its own weight. The rheology of the fluid is described by the Bingham modelwhich includes the effect of yield stress, i.e. a threshold stress which must be exceeded before flowcan occur. Both the plastic viscosity and the yield stress are modelled with temperature-dependentparameters. The flow is described by the lubrication approximation, and the heat transfer by a depth-averaged energy conservation equation. Results show that for large values of the yield stress, only theouter fraction of the fluid spreads outward, the inner fraction remaining unyielded. We also presentan analysis which predicts the threshold value of the yield stress for which partial slump occurs.1 IntroductionAs the name suggests, a viscoplastic fluid can either behave as a viscous fluid, which deformscontinuously when subject to a sufficiently strong mechanical load, or as a solid which undergoes nodeformation. As an illustrative example, a volume of viscous fluid deposited on a flat surface willspread under its own weight until the gravity and capillary forces are at equilibrium. A viscoplasticfluid, on the other hand, will only start spreading if the gravity force is strong enough to overcome thefluid cohesive force, which explains why a lump of toothpaste rests on a surface without spreading.The flow/no flow boundary for a viscoplastic fluid is dictated by the yield stress, a rheological propertyof the fluid. Flow can only occur when the stress within the fluid exceeds the yield stress. Viscoplas-tic fluids appear in a wide range of contexts including food products (e.g. ketchup, chocolate, peanutbutter), healthcare products (e.g. toothpaste, creams), the petroleum industry, cosmetics, geophysicalflows, and many others. Lava is a geophysical fluid known to have a viscoplastic behaviour Griffiths(2000), Balmforth et al. (2014). A better understanding on the rheology of lava is the motivationbehind the work presented here. Specifically, by comparing model results and laboratory/field mea-surements, we aim to indirectly infer rheological properties of the lava, which would otherwise bedifficult if not impossible to measure using standard rheometry techniques. A challenge, however, isthat as the lava flows downhill from the eruption site, it cools down by radiation to the surroundingsand conduction to the ground and many rheological parameters are temperature-dependent. Buildingon the work of Bernabeu et al. (2016), we develop and implement, a mathematical model describingthe flow of a non-isothermal viscoplastic fluid either spreading under its own weight on a horizontalsurface or draining down an inclined plane. Of particular interest here is the effect of heat transfer onthe spreading and rest state of a viscoplastic slump.Copyright is held by the author(s) through the Creative Commons BY-NC 4.0 License2 Mathematical modelA two-dimensional, gravity driven, free surface slump of fluid described by the Bingham rheo-logical law was analysed under isothermal and non-isothermal conditions. Three-dimensional andtwo-dimensional isothermal thin-film flows of a viscoplastic material have been studied extensivelyin the past, with a vast range of applications and conditions for its validity, for example see Balm-forth & Craster (1999), Balmforth et al. (2000), Mei & Yuhi (2001), Balmforth et al. (2007), Hogg& Matson (2009), Mitsoulis & Tsamapolous (2017), Hinton & Hogg (2022) and references therein.Recently, Bernabeu et al. Bernabeu et al. (2016) considered the temperature-dependent thin-film flowof a Herschel-Bulkley fluid to model lava flow down an arbitrary topography. In order to average andreduce the heat-balance equation, a closure relation is necessary. They assumed that the temperaturein the vertical direction can be modelled by a polynomial function of degree three, leading to greaterfreedom for the temperature profile but also requiring more information to constrain the solution. Inthis work, we follow instead the original idea of López et al. López et al. (1996), where they assumea second degree polynomial dependence for the vertical temperature. The result of this assumptionappears as a linear reaction term in (2). We found the difference between these approaches to beminimal and controlled via the Péclet number. The resulting unsteady one-dimensional model isht +∂x[−h2c [3h−hc]6µ( fx +hx)]= ws, (1)h [θt + ūθx]−ws(1−θ)−1Pe[2a2hθ] = 0, (2)where h(x, t) is the height of the thin film, θ(x, t) the temperature, f (x) the topography, ws(x) thesource term, and Pe the Péclet number is defined asPe =LUρcpk, (3)where L and U are the characteristic length and velocity scales, ρ the density, cp the specific heatcapacity, and k the thermal conductivity. Motion only occurs below the critical height, hc(x, t), withhc = max(0,h− B| fx +hx|). (4)The rheological parameters are the viscosity, µ(θ), and Bingham number, B(θ), which are defined byµ = eαµ(1−θ), (5)B = BieαB(1−θ), (6)where αµ, αB and Bi are constant coefficients. The Bingham number is the ratio between yield stressand viscous stress, and is defined byB =τyHµpU, (7)where τy is the yield stress, H the characteristic height scale, and µp the plastic viscosity. The averagedvelocity, ū(x, t), is given byū =−h2c(3h−hc)6µh( fx +hx). (8)2The coefficient a2(x, t) in equation (2) is calculated by solving the following system of equations.[( f +h)2 +2( f +h)− f 2]a2(x, t)+ [h+1]a1(x, t) = 0,[( f +h)3 − f 33− f 2h]a2(x, t)+[( f +h)2 − f 22− f h]a1(x, t) = h. (9)We consider homogeneous Dirichlet boundary conditions at the ends of the block for both height,h, and temperature θ. The initial conditions for both the height, hinit, and temperature, θinit, wereprescribed byr(x) ={1 −0.5 ≤ x ≤ 0.50 otherwise.3 Numerical methodA first order implicit-explicit scheme was employed for time integration, second order derivativesin space were treated implicitly, and expressions (4) to (9) were evaluated explicitly. A second ordercentral difference scheme for equation (1) was used. Equation (2) is hyperbolic, and in order to obtaina stable numerical scheme, the upwind discretisation was applied for the first derivative. The nodespacing was ∆x = 10−2 and the time step was ∆t = 10−3 for all of the cases. The simulations wererun until the block length reached a steady state; determined by evaluating the difference between twosuccessive time steps and comparing with a specified tolerance of 10−7.4 ResultsThe rheological coefficients for the non-isothermal cases were αB = αµ = 0.5, whereas the tem-perature dependence was removed for the isothermal cases by setting αB = αµ = 0.We let ws(x) = 0,and fix Pe = 500, which indicates an advection dominated flow.4.1 Horizontal and inclined flat topographyIn this section we consider isothermal and nonisothermal flows in flat topographies.4.1.1 Isothermal steady state and numerical validationFirst we present steady state analytical solutions for the isothermal case for the horizontal surfacetopography. For Bingham fluids, it is known that an arrested state exists, i.e. when all of the materialexhibits a balance between the gravitational force, pressure gradient and yield stress. For an extensivediscussion about the arrested state see Balmforth et al. (2007), Hogg & Matson (2009). The arrestedstate is fully described by equation (4)h∞ | fx +hx|= Bi, (10)which is true for any domain or topography the material is on. Note that by satisfying the balance inequation (10), the yield surface, hc, is equal to zero over the whole domain and there is no flow.For the case of a flat horizontal and flat topography we have fx = 0 and the flow if symmetric. Twoscenarios exist in this configuration, either the material fully slumps or only partially deforms. In thelatter case, part of the material remains unyielded and a yield-point exists. We present the two casesbelow,31. For a full slump, the solution satisfiesh∞(x f ,∞) = 0,∫ x f ,∞0h∞(x)dx =12. (11)Thereforeh∞(x) =√2Bi(x f ,∞ − x), (12)withx f ,∞ =(325/2)2/3B−1/3i , (13)where x f ,∞ is the position of the front of the slump.2. For a partial slump we requireh∞(xc,∞) = 1,h∞(x f ,∞) = 0,∫ x f ,∞0h∞(x)dx =12, (14)thus,h∞(x) ={√1−2Bi(x− xc) for xc,∞ < x ≤ x f ,∞1 for 0 ≤ x ≤ xc,∞(15)withxc,∞ =12− 13Bi, x f ,∞ =12+16Bi. (16)Here, xc,∞ is steady state the position of the yield point and when it is equal to zero, it provides therequired condition to find a critical Bingham number Bi,c, that is the minimum Bingham number forwhich we have a partial slump, which in this case is Bi,c = 2/3.In order to validate our numerical method we set Bi = 1, run the isothermal case to steady stateand compare with (15). The results are presented in table 1. The method is first order accurate, andconverges to the analytical solution.∆x ||h−h∞||∞0.02 0.01530.01 0.00790.005 0.00410.0025 0.0021Table 1. Validation of numerical model, listing the difference between the solution of h(t = 300) with theanalytical steady state h∞ from equation (15) under isothermal conditions with Bi = 1.4.1.2 Evolution of height and temperature profilesIn this section we compare results of isothermal and nonisothermal spreading dynamics on a hori-zontal surface, figure 1, and an inclined surface, f = m(xmax − x) in figure 2. In figures 1(a)-1(d) wepresent results for the isothermal case for two Bingham numbers, Bi = 0.6, 1. As we can see arrestedstate was reached for both Bingham numbers, clearly hc → 0 as t → ∞, see figures 1(b) and 1(d). Note4-1 -0.5 0 0.5 1x00.20.40.60.81ha)-1 -0.5 0 0.5 1x00.10.20.30.4hcb)-1 -0.5 0 0.5 1x00.20.40.60.81hc)-1 -0.5 0 0.5 1x00.050.10.150.20.25hcd)-1 0 1x00.20.40.60.81he)-1 0 1x00.10.20.30.4hcf)-1 -0.5 0 0.5 1x00.20.40.60.81g)-1 0 1x00.20.40.60.81hh)-1 0 1x00.050.10.150.20.25hci)-1 -0.5 0 0.5 1x00.20.40.60.81j)Figure 1. Horizontal surface topography ( f = 0), showing the thin film height h(x), critical height hc(x), andtemperature θ(x) profiles for times t = 0, 0.3, 3, 30, 300. Isothermal case: Bi = 0.6 for (a) and (b), Bi = 1 for(c) and (d). Non-isothermal case: Bi = 0.6 for (e) to (g), Bi = 1 for (h) to (j).that for the smaller Bingham number (Bi = 0.6) the whole block slumps and spreads until reachingsteady state. Conversely, for Bi = 1, part of the block remains static for all time, i.e. does not deform.As shown in the previous section, there are two flow regimes: (i) no yield point exists and the wholeblock deforms, figure 1(a), and (ii) xc > 0 and there is a partial slump, figure 1(c). This is consistentto our previous results since Bi,c = 2/3.The non-isothermal horizontal profiles of h(x) and hc(x), shown in figures 1(e), 1(f), 1(h) and 1(i),exhibit similar behaviour to the isothermal case. Note that as hc → 0 the average velocity ū → 0.Therefore, we expect θ → 0 as t → ∞, which is clearly the case. The Péclet number provides anindication of the ratio between advective and diffusive time scales. The problem analysed in thiswork is treated as an advected dominated flow, hence a large Pe. The temperature distribution infigures 1(g) and 1(j) has not obtained a steady state even after the film has reached an arrested state.This behaviour is driven by the large but finite Péclet number.The results for an inclined surface with m = 0.2 are shown in figure 2, where there is no longera symmetric profile for any of the variables, however, the flow characteristics are similar to the hori-zontal topography. We confirm previous published results, see Hogg & Matson (2009), and show thatan arrested state also exists if Bi ̸= 0. The length of the slump, either full or partial, is larger than inthe horizontal problem; true for both the isothermal and non-isothermal cases.Now we are in a position to investigate the effects that a temperature-dependent rheology has onBi,c. The maximum height of the block occurs at x = 0, and the height h(0) at steady state, effectivelyh∞(0), is shown in figure 3(a) for Bingham numbers in the neighbourhood of Bi = 2/3. The non-isothermal critical Bingham number is less than 2/3. Therefore, the non-isothermal slumps shouldbe shorter in length than the isothermal ones. In order to corroborate this assertion, the length ofthe steady block as a function of Bingham number for the case of a horizontal surface is shown infigure 3(b). Finally, the influence that the slope m has on the length of the block was explored for both5-1 -0.5 0 0.5 1x00.511.5h+fa)-1 -0.5 0 0.5 1x00.20.40.6hc+fb)-1 -0.5 0 0.5 1x00.511.5h+fc)-1 -0.5 0 0.5 1x00.20.40.6hc+fd)-1 0 1x00.511.5h+fe)-1 0 1x00.20.40.6hc+ff)-1 -0.5 0 0.5 1x00.511.5+fg)-1 0 1x00.511.5h+fh)-1 0 1x00.20.40.6hc+fi)-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1x00.511.5+fj)Figure 2. Inclined surface topography ( f = m(xmax − x) with m = 0.2), showing the thin film height h(x),critical height hc(x), and temperature θ(x) profiles for times t = 0, 0.3, 3, 30, 300. Isothermal case: Bi = 0.6for (a) and (b), Bi = 1 for (c) and (d). Non-isothermal case: Bi = 0.6 for (e) to (g), Bi = 1 for (h) to (j).isothermal and non-isothermal problems by fixing Bi = 1, as shown in figure 3(c). Not surprisingly,as m grows, the length grows in what seems to be a power law manner. Further investigations areunderway and remain as future work.4.2 Asymmetric topographyIn this section we consider the case of an asymmetric topography, which consists of two equidistantbumps with different heights. Hence, f (x) is defined in the following way,f (x) = ale−(x−bl)2/c2 +are−(x−br)2/c2. (17)We have now centred the block at x = 0.5 and chosen al = 0.1, bl = −0.1, ar = 0.2, br = 1.1,and c = 0.1. The values for al and ar are carefully chosen such that the lubrication approximation inequation (1) remains valid, see Kalliadasis et al. (2000).In figure 4 we present the results of our simulations. Based on results from previous sections, it isclear that we are asymptotically close to the arrested state if t ≥ 300. Therefore, we plot height, h(t,x),and temperature, θ(t,x) only for t = 300. We fix Bi = 2/3 which is the critical Bingham number forthe isothermal case. As we can see in figure 4a the symmetry is broken and two distinct yield-pointsexist. As the block starts to slump, fx(x) has opposite sign than hx(t,x) at both fronts. This increasesyield-stress effects and part of the block never moves, i.e. the critical Bingham number, Bi,c is lessthan 2/3 now. Clearly, the “new” critical Bingham number will depend on the derivative of f (x).As the block spreads, it “climbs” the bumps and stops before reaching the flat ground. Now, eventhough the sign of fx(x) is the same as hx(t,x) at the front, is not enough to overcome the effects of60.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85Bi0.910.920.930.940.950.960.970.980.991h(0)a)0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Bi11.522.53block lengthb)0 0.5 1 1.5 2m1.41.61.82block lengthc)Figure 3. (a) Maximum block height (h(0)) at steady state for Bingham number close to Bi,c. (b) Steady stateblock length as a function of Bingham number. (c) Steady state block length for increasing slope and fixedBingham number Bi = 1. Isothermal problem (black line) and non-isothermal problem (red line).the yield-stress. Something similar happens for the nonisothermal case, presented in figures 4b and4c. The cooling effects increase the unyielded region, where h(t,x)≡ 1 for all t, and the spreading ofthe block slows. As we can see, the block now stops almost a the top of the bumps. Just as in previoussections, thermal effects have an effect in the flow dynamics. This is corroborated in figure 4d wherewe plot the time series of the length of the block for both cases. As we can see, spreading has reachedsteady state for both cases and the length of the block for the nonisothermal case is about 6% smallerthan the isothermal one.5 ConclusionsThe critical Bingham number for the temperature-dependent rheology was smaller than the isother-mal counterpart which has Bi,c = 2/3. The length of the viscoplastic block was shorter for higherBingham numbers and for the non-isothermal case, and longer for steeper inclined slopes. We alsoinvestigated the case of an asymmetric topography. As in previous results, thermal effects have thegreatest effect on the flow dynamics. The dependence of the arrested state on the non-isothermal be-haviour, topographical features, and Bingham numbers has implications for the spreading dynamicsof viscoplastic flows.AcknowledgementsThis work is part of the project “Indirect measurement of lava rheology”, Marsden Fund UOC1802.The funding is kindly acknowledged.ReferencesBalmforth, N.J. and Craster, R.V., A consistent thin-layer theory for Bingham fluids, J. Non-Newtonian FluidMech., 84, 1999, 65–81.Balmforth, N.J., Burbidge, A.S., Craster, R.V., Salzig, J. and Shen, A., Viscoplastic models of isothermal lavadomes, J. Fluid Mech., 403, 2000, 37–65.Balmforth, N.J., Craster, R.V., Perona, P., Rust, A.C. and Sassi, R., Viscoplastic dam breaks and the Bostwickconsistometer, J. Non-Newtonian Fluid Mech., 142, 2007, 63–78.7-2 -1 0 1 2 3x00.20.40.60.811.2h+fa)-2 -1 0 1 2 3x00.20.40.60.811.2h+fb)-2 -1 0 1 2 3x00.10.20.30.40.50.60.7 +fc)10-410-2100102t11.11.21.31.4block lenghtd)Figure 4. Asymmetric topography, f (x) as in equation 17. (a) and (b) h(x) for isotherman and nonisothermalcase, respectively. (c) temperature θ(x) profile. For all cases t = 300 and Bi = 2/3. (d) Time series for lengthof the block for both isothermal (black curve) and nonisothermal (red curve) models.Balmforth, N J., Frigaard I.A., and Ovarlez G., Yielding to stress: recent developments in viscoplastic fluidmechanics, Annual Review of Fluid Mechanics 46, 2014, 121–146.Bernabeu, N., Saramito, P. and Smutek C., Modelling lava flow advance using a shallow-dpeth approximationfor three-dimensional cooling of viscoplastic flows, Geological Society, London, Special Publications.,1, 2016, 404–423.Griffiths. R.W., The dynamics of lava flows, Annual Review of Fluid Mechanics 32, 2000, 477–518.Hinton, E.M. and Hogg, A.J., Flow of a yield-stress fluid past a topographical feature, J. Non-Newtonian FluidMech., 299, 2022, 104696.Hogg, A.J. and Matson, G.P., Slumps of viscoplastic fluids on slopes, J. Non-Newtonian Fluid Mech., 158,2009, 101–112.Kalliadasis, S. and Bielarz, C. and Homsy, G.M., Steady free-surface thin film flows over topography, Physicsof Fluids, 12(8), 2000, 1889–1898.López, P.G., Bankoff, S. and Miksis, M.J., Non-isothermal spreading of a thin liquid film on an inclined plane,J. Fluid Mech., 324, 1996, 261–286.Mei, C.C. and Yuhi, M., Slow flow of a Bingham fluid in a shallow channel of finite width, J. Fluid Mech.,431, 2001, 135–159.Mitsoulis, E., Tsamopoulos, J. Numerical simulations of complex yield-stress fluid flows. Rheol Acta 56, 2017,231–258.8

Thin-film flow of a Bingham fluid over topography with a temperature-dependent rheology

https://ir.canterbury.ac.nz/bitstream/10092/105094/2/Moyers-Gonzalez_et_al_2022.pdf

Thin-film flow of a Bingham fluid over topography with a temperature-dependent rheology

Abstract

Similar works

Full text

Available Versions

UC Research Repository