A finite volume method for the simulation of elastoviscoplastic flows and its application to the lid-driven cavity case
AA Finite Volume method for the simulation of elastoviscoplastic flows and itsapplication to the lid-driven cavity case
Alexandros Syrakos a, ∗ , Yannis Dimakopoulos a , John Tsamopoulos a a Laboratory of Fluid Mechanics and Rheology, Department of Chemical Engineering, University of Patras, Patras, Greece
Abstract
We propose a Finite Volume Method for the simulation of elastoviscoplastic flows, modelled after the ex-tension to the Herschel-Bulkley model by Saramito [J. Non-Newton. Fluid Mech. 158 (2009) 154–161]. Themethod is akin to methods for viscoelastic flows. It is applicable to cell-centred grids, both structuredand unstructured, and includes a novel pressure stabilisation technique of the “momentum interpolation”type. Stabilisation of the velocity and stresses is achieved through a “both sides diffusion” technique andthe CUBISTA convection scheme, respectively. A second-order accurate temporal discretisation schemewith adaptive time step is employed. The method is used to obtain benchmark results of lid-driven cavityflow, with the model parameters chosen so as to represent Carbopol. The results are compared againstthose obtained with the classic Herschel-Bulkley model. Simulations are performed for various lid velocities,with slip and no-slip boundary conditions, and with different initial conditions for stress. Furthermore, weinvestigate the cessation of the flow, once the lid is suddenly halted.
Keywords: elastoviscoplastic flow; Finite Volume method; Carbopol; lid-driven cavity; benchmarkproblem
1. Introduction
Viscoplastic (VP) fluids are a class of materials whose distinctive property is that they flow as fluids ifsubjected to large enough stresses but behave as solids if the applied stress is below a critical value, termedthe yield stress. Although solids can also undergo plastic deformation, viscoplastic fluids are characterisedby reversibility of the structural changes caused during plastic flow once the flow ceases [1, 2]. The class ofVP fluids includes a variety of materials such as foams, emulsions, colloids and physical gels, with possiblydifferent microscopic mechanisms being responsible for the emergence of a yield stress in each case [3].Viscoplastic flows are of major relevance in many industries (oil, construction, cosmetics, foodstuffs, etc.)and many natural flows can be classified as such (mud and lava flows, landslides, avalanches etc.).Viscoplastic fluids were first studied in depth by Eugene Bingham [1] who proposed the renowned consti-tutive equation named after him to describe their behaviour [4]. A short time later, Herschel and Bulkley [5]extended the Bingham model to describe also the shear-thinning (or shear-thickening) post-yield behaviourthat most of these materials exhibit, by assuming a power-law dependence of the viscosity on the shear rate.These models were originally proposed in scalar form, but it was not long until a full tensorial form wasproposed, employing the von Mises yield criterion [6]. The empirical Herschel-Bulkley (HB) constitutiveequation has been found to represent well the behaviour of yield stress fluids under steady shear flow, andis arguably the most popular model for VP fluids. It is commonly given in the following form: τ ≤ τ y ⇒ ˙ γ = 0 τ > τ y ⇒ τ = (cid:18) τ y ˙ γ + k ˙ γ n − (cid:19) ˙ γ (1)where τ is the deviatoric stress tensor; ˙ γ = ∇ u + ( ∇ u ) T is the rate-of-strain tensor, u being the fluid velocityvector and T denoting the tensor transpose; and τ ≡ ( τ : τ ) and ˙ γ ≡ ( ˙ γ : ˙ γ ) are their magnitudes. ∗ Corresponding author
Email addresses: [email protected], [email protected] (Alexandros Syrakos), [email protected] (Yannis Dimakopoulos), [email protected] (John Tsamopoulos)
Preprint submitted to Journal of Non-Newtonian Fluid Mechanics April 23, 2019 a r X i v : . [ phy s i c s . f l u - dyn ] A p r he HB model includes three parameters, the yield stress τ y , the consistency index k , and the exponent n ;the latter is usually in the range 0 . − . τ < τ y ) or fluid ( τ > τ y ),separated by the yield surface ( τ = τ y ). In order to express the von Mises yield criterion, τ in Eq. (1) isdefined to be the deviatoric component of the total stress tensor σ ; such a decomposition of σ into deviatoric(traceless), σ d , and isotropic, σ i , components is useful for solids: σ = σ −
13 tr( σ ) I (cid:124) (cid:123)(cid:122) (cid:125) σ d + 13 tr( σ ) I (cid:124) (cid:123)(cid:122) (cid:125) σ i (2)where tr( σ ) ≡ (cid:80) i σ ii is the trace of σ . Therefore, in the common form (1) of the HB equation, τ identifieswith σ d of Eq. (2). However, from a constitutive equation for a fluid one expects a decomposition of σ intothe extra-stress tensor (for which the symbol τ is usually reserved) which expresses the forces that arise dueto deformation of the fluid, and an isotropic pressure component which is responsible for enforcing continuityin incompressible fluids: σ = τ − pI (3)where I is the identity tensor. Equation (1) does not suffice to define the HB extra-stress tensor because itallows the possibility that this tensor has an isotropic part, for which no information is given. This ambiguityis removed by the tacit assumption that the HB extra-stress tensor is traceless, so that decompositions (2)and (3) become equivalent and τ of Eq. (1) also happens to equal the extra-stress tensor (hence the symbol τ instead of σ d in (1)). This makes the HB fluid a generalised Newtonial one, for which the extra anddeviatoric stresses are identical since the former is traceless due to the incompressible continuity equation.When elasticity is introduced into the HB equation in Sec. 2, the extra-stress τ will no longer be tracelessand thus can no longer identify with σ d .In the solid region, the concepts of extra-stress and pressure are not meaningful, since the HB solid doesnot deform, and the continuity equation is already enforced by the rigidity condition without a need forpressure. The symbols τ and p are still used there, but it must be kept in mind that they simply refer tothe deviatoric and isotropic components of σ . In other words, in the HB solid the stress decomposition iswritten as in Eq. (3), but what is really meant is the decomposition (2).With some manipulations, noticing that in the fluid branch the stress tensor and the rate-of-strain tensorare parallel and thus the unit tensors ˙ γ/ ˙ γ and τ /τ are equal, the solid and fluid branches of Eq. (1) can becombined into a single expression for ˙ γ [7]:˙ γ = (cid:18) max(0 , τ − τ y ) k (cid:19) n τ τ (4)Viscoplastic constitutive equations such as the HB have been studied extensively during the past decades.Their discontinuity at the yield surfaces and their inherent indeterminacy of the stress tensor within theunyielded regions require specialised numerical techniques for performing flow simulations [8–10]. Further-more, there is the question of whether their assumptions of a completely rigid (inelastic) solid phase and apurely viscous fluid phase are physically realistic. Allowing for some deformation of the solid phase understress seems more natural and indeed experimental studies on Carbopol, a prototypical material often usedin experimental studies on viscoplasticity, have shown that prior to yielding it exhibits elastic deformationunder stress [11, 12]. Elastic effects can be observed also in the fluid phase. For example, bubbles risingin Carbopol solutions usually acquire the shape of an inverted teardrop, with a cusp at their leeward side[13, 14], but the classical Bingham and HB viscoplastic models fail to predict such behaviour [15, 16]; on theother hand, such shapes are observed also in viscoelastic fluids, and are correctly captured by viscoelasticconstitutive equations [17]. Similarly, for the settling of spherical particles in Carbopol, classical VP modelscannot predict phenomena such as the loss of fore-aft symmetry under creeping flow conditions and the for-mation of a negative wake behind the sphere, but these phenomena are predicted if elasticity is incorporatedinto the constitutive modelling [18].Therefore, recently the focus has been shifting towards constitutive equations that incorporate bothplasticity and elasticity, usually called elastoviscoplastic (EVP) constitutive equations. Actually, EVP con-stitutive modelling dates back to the beginning of the previous century – a nice historical overview can be2ound in [19]. Although several EVP models have been proposed (see e.g. the literature reviews in [19, 20]),many of them appear only in scalar form. To be applicable in complex two- and three-dimensional flowsimulations (2D/3D), a full tensorial form is required; some such models are compared in [21].Complex simulations also require that these models be accompanied by appropriate numerical solverscharacterised by accuracy, robustness and efficiency. Usually, FEM solvers are employed in EVP flow sim-ulations (e.g. [22, 23, 21]). An alternative, very popular discretisation / solution method in ComputationalFluid Dynamics is the Finite Volume Method (FVM); it has been successfully applied to viscoplastic (e.g.[24, 25]) and viscoelastic (e.g. [26–29]) flows individually, but not to EVP flows, to the best of our knowl-edge (a hybrid FE/FV method was used in [30]). In this paper a FVM for the simulation of EVP flows isdescribed. The EVP constitutive equation chosen is that of Saramito [31] which introduces elasticity intothe classic HB model, to which it reduces in the limit of inelastic behaviour. This model shall be referred toas the Saramito-Herschel-Bulkley (SHB) model. We chose this model because of its simplicity, its potentialas revealed in a number of recent studies on materials such as foams [32] and Carbopol [33], and because ofthe popularity of the classic HB model. Nevertheless, the general framework of the presented FVM shouldbe applicable to a range of other EVP models, particularly those that can be regarded as modifications /extensions of viscoelastic constitutive equations.The presentation of the method in Sections 3 and 4 and its validation in Sec. 5 are followed by applicationof the method to simulate EVP flow in a lid-driven cavity. The lid-driven cavity test case is arguably themost popular benchmark test case for new numerical methods for flow simulations. As such, it has beenused also as a benchmark problem for viscoplastic [24, 34] and viscoelastic [35] flows; there exist also theEVP flow studies of [36, 23], but with a different EVP model that incorporates a kind of regularisation.In the present study the parameters of the SHB model are chosen so as to represent Carbopol, which isregarded as a simple VP fluid (more complex behaviour such as thixotropy [37] and kinematic hardening [38]are not considered, but may be incorporated into the model in the future). The lid-driven cavity problemconstitutes a convenient “playground” for testing the numerical method, testing the behaviour of the SHBmodel under conditions of complex 2D flow, comparing its predictions against those of the classic HB model,and providing benchmark results. The tests include varying the lid velocity to vary the flow character asquantified by dimensionless numbers, flow cessation (which occurs in finite time for VP flows), and varyingthe initial conditions to investigate the issue of multiplicity of solutions of the SHB model.
2. Governing equations
The flow is governed by the continuity, momentum, and constitutive equations. The first two are: ∇ · ( ρu ) = 0 (5) ∂ ( ρu ) ∂t + ∇ · ( ρuu ) = −∇ p + ∇ · τ (6)where t is time, and ρ is the density of the material, assumed to be constant; the rest of the variables havebeen defined in Sec. 1. The right-hand side of Eq. (6) can be written collectively as ∇ · σ .Closure of the system of governing equations requires that the extra stress tensor be related to the flowkinematics through a constitutive equation. In the present work we use the EVP equation of Saramito[31], which assumes that the total deformation of the material is equal to the sum of an elastic strain γ e (applicable in both the solid and fluid phases) and the provisional (if the yield stress is exceeded) viscousdeformation γ v predicted by the HB model (4). This behaviour is depicted schematically in Fig. 1 by amechanical analogue. In terms of the rate-of-strain this is written as:1 G (cid:79) τ (cid:124)(cid:123)(cid:122)(cid:125) ˙ γ e + (cid:18) max(0 , τ d − τ y ) k (cid:19) n τ d τ (cid:124) (cid:123)(cid:122) (cid:125) ˙ γ v = ˙ γ (7)where G is the elastic modulus, the triangle denotes the upper-convected time derivative (cid:79) τ ≡ ∂τ∂t + u · ∇ τ − ( ∇ u ) T · τ − τ · ∇ u (8)3 igure 1: A mechanical analogue of the Saramito-Herschel-Bulkley (SHB) model [31]. and τ d ≡ ( τ d : τ d ) is the magnitude of the deviatoric part of the EVP stress tensor, τ d ≡ τ −
13 tr( τ ) I (9)Note that because τ and σ differ only by an isotropic quantity ( − pI , Eq. (3)) it holds that τ d = σ d . Animportant difference with the classic HB model (4) is that, due to the last two terms of the upper convectedderivative (8), the trace of τ is now not necessarily zero and thus, in general, τ d (cid:54) = τ , necessitating explicituse of τ d inside the “max” term of Eq. (7) in order to express the von Mises yield criterion. The full 3Dformula (9) must be used for the deviatoric stress even in 1D or 2D simulations: a 2D stress state where all τ ij are zero except τ = τ (cid:54) = 0 is not isotropic ( τ (cid:54) = τ , τ ) and τ d is not zero.Other qualitative differences with the HB model exist. For example, the SHB model allows non-zerorate of strain in the unyielded regions, arising from elastic deformations of the solid phase. Conversely, it istheoretically possible that the rate of strain is zero in yielded material because the stresses there have nothad enough time to relax. In particular, in the SHB model, contrary to the HB model, the stresses do notrespond instantaneously to changes in the rate of strain; rather, the EVP stress tensor is also a function ofall past states of ˙ γ , as in viscoelastic fluids.Furthermore, as noted in [39, 32], the combination of plastic and elastic terms in the SHB Eq. (7) canproduce complex behaviour not observed in either purely viscoplastic or purely viscoelastic flows. One suchfeature is that the extra stress tensor and the velocity gradient can vary discontinuously across yield surfaces[39]. Another feature is that flows are inherently transient, in the sense that there exist an infinitude ofsteady states which depend in a continuous manner on the initial conditions. Even if it is just the steadystate that is sought, one cannot simply discard the time derivatives in Eqs. (6) and (7); the steady statemust be obtained by a transient simulation with appropriate initial conditions. It is not only the steadystate stresses that depend on the initial conditions, but also the steady state velocity. This issue was studiedin [32] in the context of cylindrical Couette flow. In fact, actual EVP materials do behave this way inexperiments, with the steady state depending on the residual stresses present in the unyielded, stationarymaterial at the start of the experiment [32]. The residual stresses are stresses that are “trapped” insideunyielded material because there is no relaxation mechanism there. E.g. for stationary ( u = 0) unyieldedmaterial, Eq. (7) predicts ∂τ /∂t = 0. In experiments, residual stresses do develop during the preparation ofthe material and are difficult to eliminate.This contrasts the behaviour of classic VP models such as the HB, where there is a stress indeterminacyin the unyielded regions: there exist infinite σ fields within these regions that have the same divergence ∇ · σ which yields ˙ γ = 0 when substituted in the momentum equation (6). Each of these fields is a valid solution.The HB steady-state does not depend on the initial conditions, and can be obtained directly, by droppingthe time derivative from Eq. (6). The HB model makes no connection between this indeterminacy and theinitial conditions, and in fact it even allows that stresses in the unyielded regions vary discontinuously intime. No such indeterminacy is exhibited by the SHB model, the stresses in both the yielded and unyieldedregions being precisely determined for given initial conditions.In the limit of very high elastic modulus, G → ∞ , the SHB material becomes so stiff that it behaves as aninelastic fluid or solid; the importance of the first term on the left-hand side of the SHB constitutive equation(7) diminishes and that equation tends to become identical to the classic HB equation (4). However, it mustbe kept in mind that the difference in character between the SHB and HB models is in some respects retained4o matter how large the value of G is. This concerns mostly the stress state in the unyielded regions, whichfor the SHB model remains uniquely determined by the initial conditions whereas that of the HB model isindeterminate. Also, as discussed in Appendix A, in the HB model τ was defined to equal σ d whereas forthe SHB model the identification of τ with σ d as G → ∞ is not guaranteed.Figure 1 shows that the complete SHB model contains an additional viscous component labelled κ notdiscussed thus far. This is a Newtonian component, of viscosity κ , which makes the unyielded phase behaveas a Kelvin-Voigt solid. The entire extra stress tensor is then τ e = τ + κ ˙ γ (10)The SHB model then reduces to the Oldroyd-B viscoelastic model when τ y = 0 and n = 1. However, in themain results of Sec. 6 we will use κ = 0.Finally, concerning the boundary conditions, we consider solid wall boundaries. We will mostly employthe no-slip condition, but since (elasto)viscoplastic materials are usually slippery we will also employ a Navierslip condition, according to which the relative velocity between the fluid and the wall, in the tangentialdirection, is proportional to the tangential stress. For two-dimensional flows this is expressed as follows: Let n be the unit vector normal to the wall, and s be the unit vector tangential to the wall within the plane inwhich the equations are solved. Let also u and u w be the fluid and wall velocities, respectively. Then,( u − u w ) · s = β (cid:0) n · τ (cid:1) · s (11)where the parameter β is called the slip coefficient. In the present work we will solve the dimensional governing equations. Nevertheless, nondimensionali-sation of the equations reveals a number of dimensionless parameters that characterise the flow. So, let L and U be characteristic length and velocity scales of the problem, with T = L/U the associated time scale,and choose the following characteristic value of stress, S (the Newtonian viscosity κ , Eq. (10), is omitted): S ≡ τ y + k (cid:18) UL (cid:19) n (12)Then, we define the dimensionless variables, denoted with a tilde (˜), as x i = ˜ x i L , t = ˜ tT , u = ˜ uU , p = ˜ pS , and τ = ˜ τ S . Substituting these into the governing equations (5), (6), and (7)-(8), we obtain theirdimensionless forms: ˜ ∇ · ˜ u = 0 (13) Re (cid:18) ∂ ˜ u∂ ˜ t + ˜ ∇ · (˜ u ˜ u ) (cid:19) = − ˜ ∇ ˜ p + ˜ ∇ · ˜ τ (14) Wi (cid:18) ∂ ˜ τ∂ ˜ t + ˜ u · ˜ ∇ ˜ τ − ( ˜ ∇ ˜ u ) T · ˜ τ − ˜ τ · ˜ ∇ ˜ u (cid:19) + (cid:18) max (0 , ˜ τ d − Bn )1 − Bn (cid:19) n τ d ˜ τ = ˜˙ γ (15)where ˜ ∇ = (1 /L ) ∇ and ˜˙ γ = ˙ γ/ ( U/L ). The following dimensionless numbers appear:Reynolds number Re ≡ ρU S (16)Weissenberg number Wi ≡ SG (17)Bingham number Bn ≡ τ y S (18)A more standard choice for S would account only for viscous effects, omitting plasticity: S = k ( U/L ) n .This would lead to more standard definitions of the Reynolds and Bingham numbers: Re (cid:48) ≡ ρU − n L n k = Re − Bn (19)5 n (cid:48) ≡ τ y L n kU n = Bn − Bn (20)However, it was shown in [40] (see also the discussion in [41]) that in viscoplastic flows Re suffices as astandalone indicator of inertial effects in the flow, whereas Re (cid:48) does not (the inertial character must beinferred from the values of Re (cid:48) and Bn (cid:48) combined). Hence the definition (12) was preferred. Concerning theBingham number, Bn ∈ [0 ,
1] and Bn (cid:48) ∈ [0 , ∞ ) are simply related through Eq. (20) and carry exactly thesame information. Because of familiarity with Bn (cid:48) that is almost universally used in the literature, we willuse both Bn and Bn (cid:48) in this study.The Weissenberg number Wi is an indicator of elastic effects in the flow. Unlike the Re and Bn numbers,it is not a ratio between stresses of different nature or momentum fluxes. In fact, as seen in the mechanicalanalogue of Fig. 1, omitting the Newtonian component κ , the viscoplastic and elastic components areconnected in series and thus carry the same load. Therefore, the representative stress S , although definedby considerations pertaining to the viscoplastic component of the material behaviour (Eq. (12)), is bornealso by the elastic component of the material structure and thus Wi ≡ S/G is a typical elastic deformation.In the literature, the standard form for the Weissenberg number is Wi = λU/L , where λ is a relaxationtime, which becomes equivalent to our definition if we define λ ≡ ηG where η ≡ SU/L (21)The relaxation time λ is proportional to the apparent viscosity η (the more viscous the flow is, the slowerthe recovery from elastic deformations) and inversely proportional to the elastic modulus G (the stiffer thematerial is, the faster it recovers). Note that λ and η as defined by (21) are not material constants butreflect also the influence of the flow, as they depend on U and L in addition to the material parameters τ y , G , k and n (Eq. (12)). The definition Wi = λU/L is also interpreted as a typical elastic strain: it isthe product of the strain rate U/L by the time period λ during which the material has not yet significantlyrelaxed.Finally, an important dimensionless number is the yield strain, γ y = τ y G = Bn · Wi (22)which depends only on the material parameters and not on kinematic or geometric parameters of the flowsuch as U and L . The fact that it is equal to the product of Bn and Wi means that these two numbersare not independent but their product is constant for a given material. Since Bn ∈ [0 , Wi ∈ [ γ y , ∞ ). That Wi is bounded from below by γ y follows also directly from the definitions (17) and (12),and derives from the characteristic stress definition (12) which assumes the existence of unyielded regions;in such regions the strain, of which Wi is a measure, is higher than the yield strain γ y .
3. Discretisation of the governing equations
In this section we propose a Finite Volume Method (FVM) for the discretisation of the governing equa-tions. The method will be described for two-dimensional problems, but extension to three dimensions isstraightforward. The first step is the tessellation of the domain Ω into a number of non-overlapping volumes,or cells . Each cell is bounded by straight faces, each of which separates it from another single cell or fromthe exterior of the domain (the latter are called boundary faces). Figure 2 shows such a cell P along withits faces and neighbours, and the associated nomenclature. Our FVM is applicable to grids of arbitrarypolygonal cells, although for the results of Sec. 6 we will employ Cartesian grids due to the regularity of theproblem geometry. Grids will be labelled after a characteristic cell length h .The discretisation procedure will convert the governing equations into a large set of algebraic equationsinvolving only (approximate) values of the dependent variables at the cell centroids. These values will bedenoted using the cell index as a subscript, i.e. φ P is the value of the variable φ at the centroid of cell P ; ifthe variable name already includes a subscript then the cell index will be separated by a comma (e.g. τ d,P is the deviatoric stress at cell P ). We aim for a second-order accurate method, i.e. one whose discretisation6 (cid:48) n d c P N N N N c m Figure 2:
Cell P and its neighbouring cells, each having a single common face with P . Its faces and neighbours arenumbered in anticlockwise order, with face f separating P from its neighbour N f . The shaded area lies outside thedomain, so face 5 is a boundary face. The geometric characteristics of face 1 are displayed. The position vectors ofthe centroids of cells P and N f are denoted as P and N f ; c f is the centroid of face f and c (cid:48) f is its closest point on theline connecting P and N f ; m f is the midpoint between P and N f ; n f is the unit vector normal to face f , pointingoutwards of P , and d f is the unit vector pointing from P towards N f . The volume of cell P is denoted as Ω P . (cid:78)(cid:78) (cid:78)(cid:78)(cid:78)(cid:78) (cid:78)(cid:78)(cid:72)(cid:72) (cid:72)(cid:72)(cid:72)(cid:72) (cid:72)(cid:72) Figure 3:
A checkerboard distribution of values on a Cartesian grid. Different values are stored at points (cid:78) and (cid:72) . error scales as h . With the present governing equations, a difficulty arises from the fact that they all containonly first spatial derivatives because this allows the development of spurious oscillations in the solution. Forexample, assume that on the grid of Fig. 3, a quantity φ has the value zero at the boundary points ( • ),the value +1 at points ( (cid:78) ), and the value − (cid:72) ). FVM integration of first derivatives of φ overeach cell eventually comes down to a calculation involving values of φ at face centres ( ◦ ) and ( • ), due toapplication of the Gauss theorem. If we use linear interpolation to approximate the values of φ at the facecentroids ( ◦ ) from the values at the cell centres on either side of each face, then we obtain φ = 0 at everyface centre ( ◦ ). This ultimately results in obtaining a zero integral of the first derivative of φ over any cell.Thus, interpolation filters out the oscillating cell-centre field and leaves a smooth field over the face centres,which in turn produces an image of the discrete operator that varies smoothly from one cell to the next.To see how this creates a problem, consider the FVM discretisation on grid h of the following PDE f ( φ ) = b (23)where b is a known function. The FVM produces a system of algebraic equations that can be written as F h ( φ h ) = B h (24)where φ h is the vector of unknown values of φ at the cell centroids, F h is the discrete operator representingthe integral of the differential operator f over each cell, and B h is the vector of integrals of the knownfunction b over each cell. The aforementioned filtering action of the face interpolation scheme means that F h will filter out any oscillations in φ h , producing a smooth image F h ( φ h ); therefore, the solution φ h to the7ystem (24) can be oscillatory even if the right-hand side B h is smooth. This is not just a remote possibility,but it does occur in practice, a notorious example being the spurious pressure oscillations produced bythe FVM solution of the Navier-Stokes equations (see [42] for a demonstration). But, if the discretisationschemes incorporated in the operator F h are designed so as to reflect any oscillations in the vector φ h ontothe image F h ( φ h ), then for a smooth right-hand side B h the system (24) can only have a smooth solution φ h . For the Navier-Stokes equations, the main concern is the pressure oscillations because any velocityoscillations will be reflected onto the second derivatives of velocity present in the momentum equations andare therefore inhibited. But in our case, all three of the governing equations, (5), (6) and (7)–(8), onlycontain first derivatives and thus we have to be concerned about the possibility of spurious oscillations inall three variables u , p , and τ . In the following subsections the adopted measures will be described.In what follows we will need a “characteristic viscosity”, a quantity with units of viscosity somewhatcharacteristic of the flow’s viscous character. Our first choice is the following: η a ≡ κ + SU/L (25)where S is given by (12). This value is constant throughout the domain. A second option we tried is similarto the coefficient of the DAVSS-G technique proposed in [43] and varies throughout the domain, dependingon the ratio between the magnitudes of the stress and rate-of-strain tensors: η a ≡ κ + 12 SU/L + 12 S + τ P + τ N f U/L + ˙ γ P + ˙ γ N f (26)The above formula gives η a at face f of cell P . With definition (26), η a never falls below κ + 0 . S/ ( U/L ),tends to κ + S/ ( U/L ) when both τ and ˙ γ are small, and tends to κ + 0 . S/ ( U/L ) + τ / ˙ γ when both τ and˙ γ are large. Also, it does not tend to infinity when ˙ γ →
0, although it has no upper bound.The discretisation of several terms (e.g. the calculation of ˙ γ ≡ ∇ u +( ∇ u ) T and its magnitude) requires theuse of a discrete gradient operator which approximates ∇ φ at the cell centres. We employ the least-squaresoperator described in detail in [44], denoted here as ∇ qh , the superscript q being the exponent employed;in [44] it was shown that on smooth structured grids the choice q = 1 . ∇ . h φ = ∇ φ + O ( h )) while with other q values the accuracy degrades to first-order at boundary cells. Onirregular (unstructured) grids all choices of q result in first order accuracy everywhere. Nevertheless, first-order accurate gradients suffice for second-order accuracy of the differentiated variables, i.e. are compatiblewith second-order accurate FVMs [44].Finally, linear interpolation for calculating face centre values on unstructured grids is performed as¯ φ c f = ¯ φ c (cid:48) f + ¯ ∇ qh φ c (cid:48) f · ( c f − c (cid:48) f ) (27)where the overbar denotes an interpolated value, and in the right-hand side both φ and its gradient areinterpolated linearly to point c (cid:48) f (the closest point to the face centre lying on the line joining P to N f , Fig.2) according to the formula ¯ φ c (cid:48) f ≡ (cid:107) c (cid:48) f − N f (cid:107)(cid:107) N f − P (cid:107) φ P + (cid:107) c (cid:48) f − P (cid:107)(cid:107) N f − P (cid:107) φ N f (28) Integrating Eq. (5) over cell P and applying the divergence theorem, we get the mass flux balance forthat cell: (cid:88) f (cid:90) s f n · ( ρu ) d s = 0 (29)where s f is the surface of face f , n is the outward (of P ) unit vector normal to that face, and d s is aninfinitesimal surface element. The integrals summed are the outward mass fluxes through the respectivefaces of cell P . They are approximated by midpoint rule integration with additional stabilising terms: (cid:90) s f n · ( ρu ) d s ≈ ρ s f (cid:104) ¯ u c f · n f + u p + f − u p − f (cid:105) ≡ ˙ M f (30)8here u p + f ≡ a mif (cid:0) p P − p N f (cid:1) (31) u p − f ≡ a mif ¯ ∇ qh p m f · ( P − N f ) , ¯ ∇ qh p m f ≡ (cid:0) ∇ qh p P + ∇ qh p N f (cid:1) (32) a mif ≡ ρ (cid:18) (cid:107) u N f − u P (cid:107) + h f ∆ t (cid:19) + 2 η a h f , h f ≡ (cid:20) (cid:0) Ω P + Ω N f (cid:1)(cid:21) D (33)where ∆ t is the current time step (see Sec. 3.5), Ω P and Ω N f are the cell volumes, and D equals either 2or 3, for 2D and 3D problems, respectively. With the definition (30), the discrete version of the continuityequation for a cell P is (cid:88) f ˙ M f = 0 (34)To expound the above scheme, we make the following remarks: Remark 1: If u p + f and u p − f are omitted from Eq. (30) then we are left with simple midpoint ruleintegration (¯ u c f is obtained using the scheme (27)). The part of ˙ M f due to u p + f and u p − f can be viewed asbelonging to the truncation error of the continuity equation of cell P :1Ω P ρ s f (cid:16) u p + f − u p − f (cid:17) = ρ s f Ω P a mif (cid:2)(cid:0) p P − p N f (cid:1) − ¯ ∇ qh p m f · ( P − N f ) (cid:3) (35)(the truncation error is defined per unit volume, hence we divide by the cell volume Ω P ). It is easy to show,by expanding the pressure in a Taylor series about point m f (Fig. 2), that ( p P − p N f ) − ∇ p ( m f ) · ( P − N f ) = O ( h ). However, since we use only an approximation to the pressure gradient, ∇ qh p = ∇ p + O ( h ), the termin square brackets in Eq. (35) is O ( h ) (because O ( h ) · ( P − N f ) = O ( h )). Given also that Ω P = O ( h ), s f = O ( h ) and a mif = O ( h ), the whole term (35) is O ( h ), which is compatible with a second-order accuratemethod such as the present one. Remark 2:
The term (35) inhibits spurious pressure oscillations by reflecting them on the image of thediscrete continuity operator (see Eq. (24) and associated discussion). Pressure oscillations cause the term u p + f to oscillate from face to face, and are thus reflected on the continuity image. On the other hand, suchoscillations are filtered out by the gradient operator ∇ qh , so that the term u p − f varies smoothly from face toface (it does not pass the oscillations on to the image). The u p − f term is used simply for counterbalancingthe truncation error introduced by u p + f . Remark 3:
The coefficient a mif is chosen so that the terms u p + f and u p − f of the mass flux expression(30), which have units of velocity, are neither too small to have a stabilising effect nor so large that theydominate the mass flux. The “velocities” u p + f and u p − f are functions of local pressure variations, and attemptto quantify the contributions of these pressure variations to the velocity field. Pressure and velocity areconnected through the momentum equation, so consider this equation in the following non-conservative form,where we assume that the stress tensor can be approximated through the use of a characteristic viscositysuch as (25) or (26) as τ ≈ η ( ∇ u + ( ∇ u ) T ): ρ (cid:18) ∂u∂t + u · ∇ u (cid:19) = −∇ p + ∇ · ( η ∇ u ) (36)where we have neglected the term ∇ · ( η ( ∇ u ) T ), assuming that it is small (it is zero if η is constant). Weare free to make these sorts of approximations because all we want is a rough estimate of the effect that thepressure gradient has on velocity. So, consider the simple uniform grid, of spacing h , shown in Fig. 4, where u denotes the velocity component normal to the face separating cells P and N . We will employ a simpleFV discretisation of Eq. (36) in order to relate the velocity u c at the face centre c to the pressures at thecentres of the adjacent cells P and N . The momentum conservation, Eq. (36), in the direction x normal tothe face, for the imaginary cell drawn in dashed line surrounding that face, can be discretised as ρ u c − u oldc ∆ t h + ρu c d u d x (cid:12)(cid:12)(cid:12)(cid:12) c h = − p N − p P h h + 1 h (cid:18) η N u c N − u c h − η P u c − u c P h (cid:19) h (37)9 P P N N Nh u c cc P c N Figure 4:
A row of grid cells. where u old c was the velocity at time ∆ t ago. Assuming that ( η N + η P ) u c ≈ η c u c we can solve this for u c : u c = 1 ρ (cid:18) d u d x (cid:12)(cid:12)(cid:12)(cid:12) c h + h ∆ t (cid:19) + 2 η c h ( p P − p N ) + · · · (38)where the dots ( · · · ) denote the terms that are not related to pressure. The above equation provides a quan-tification of the local effect of pressure gradient on velocity. It was derived from a simplistic one-dimensionalconsideration; for more general flows, the coefficient a mif (33) can be seen to be a generalisation of the coeffi-cient multiplying the pressure difference in Eq. (38). It should be noted that the one-dimensional momentumequation (37) accounts for the viscous force due to the velocity variation in the direction perpendicular tothe face, but omits that due to the velocity variation in the direction parallel to the face; had the latter alsobeen accounted for, the viscous term in the denominator of (38), and of a mif (33), would have been 4 η/h (or6 η/h in 3D) instead of 2 η/h , which would have reduced the magnitude of the stabilisation terms in the massflux scheme (30), but would also likely somewhat increase the accuracy, according to the results of Sec. 5.We did not investigate this issue further, but the choice is not crucial as it affects neither the stabilisationability of the technique nor the O ( h ) magnitude of the error (35) that it introduces. Remark 4:
An explanation of how the scheme retains its oscillation-inhibiting effect when h → u p + f , u p − f tending to zero is given in Appendix B.The scheme (30) – (33) is a variant of the popular technique known as momentum interpolation , whichwas originally proposed in [45]. Ever since, many variants of this technique have been proposed (see [46]and references therein) almost all of which are intertwined with the SIMPLE algebraic solver (exceptionsinclude [47, 42]). Although this connection with SIMPLE can be useful in some respects [48], and SIMPLEis employed in the present work, we prefer an independent method such as the presently proposed becauseit is more general, transparent, and easily adaptable. It can be used with other algebraic solvers such asNewton’s method, and it does not lead to dependence of the solution on underrelaxation factors. As for the continuity equation, the FVM discretisation of the momentum equation (6) begins by inte-grating it over a cell P and applying the divergence theorem, to get (cid:90) Ω P ∂ ( ρu ) ∂t dΩ + (cid:88) f (cid:90) s f n · ( ρuu ) d s = − (cid:88) f (cid:90) s f n p d s + (cid:88) f (cid:90) s f n · τ d s (39)Using midpoint rule integration, the above equation is approximated by ∂ ( ρu ) ∂t (cid:12)(cid:12)(cid:12)(cid:12) aP Ω P + (cid:88) f ˙ M f ¯ u c f = − (cid:88) f ¯ p c f s f n f + (cid:88) f F τf (40)where ∂/∂t | aP , the approximate time derivative at P , will be defined in Sec. 3.5, ˙ M f is the mass outfluxthrough face f defined by Eq. (30), ¯ u c f and ¯ p c f are approximated from cell-centre values with the interpola-tion scheme (27), and F τf is the approximation of the force on face f due to the EVP stress tensor τ . Notethat Eq. (40) is a vector equation, with two (in 2D) or three (in 3D) scalar components. The force F τf iscalculated as (cid:90) s f n · τ d s ≈ F τf ≡ s f (cid:104) n f · ¯ τ c f + D τ + f − D τ − f (cid:105) (41)10gain, the stress tensor at the face centroid, ¯ τ c f , is calculated via the linear interpolation scheme (27)(applied component-wise). The viscous pseudo-stresses D τ + f and D τ − f are stabilisation terms employed tosuppress spurious velocity oscillations. They are vectors, whose i -th components are D τ + f,i ≡ η a u i,N f − u i,P (cid:107) N f − P (cid:107) (42) D τ − f,i ≡ η a ¯ ∇ qh u i,m f · d f , ¯ ∇ qh u i,m f ≡ (cid:0) ∇ qh u i,P + ∇ qh u i,N f (cid:1) (43)where u i is the i -th velocity component ( u = ( u , u ) in 2D), and d f = ( N f − P ) / (cid:107) N f − P (cid:107) is the unit vectorpointing from P to N f (Fig. 2). It can be seen that both D τ + f,i and D τ − f,i equal a characteristic viscosity, givenby (25) or (26), times the velocity gradient in the direction d f , albeit calculated differently: the gradient ascomputed in D τ + f,i is sensitive to velocity oscillations whereas that computed in D τ − f,i is not. The mechanismof oscillation suppression is similar to that for momentum interpolation: in the presence of spurious velocityoscillations, the smooth part of D τ + f,i is cancelled out by D τ − f,i , but the oscillatory part produces oscillationsin the operator image (in F h ( φ h ), in the terminology of Eq. (24)).In the context of colocated FVMs for viscoelastic flows, this technique was first proposed in [26], in-spired from the corresponding “momentum interpolation” technique for pressure. A question concerning themethod is the appropriate choice of the parameter η a . The original method [26] as well as some subsequentvariants (e.g. [49, 50]), similarly to the original momentum interpolation [45], derived the coefficient fromthe SIMPLE matrix of the linearised discrete constitutive equation, arriving at complicated expressions.The present simpler approach is essentially equivalent to that adopted by [51, 52] who, for viscoelastic flows,set the coefficient η a equal to the polymeric viscosity. The aim is that D τ + f and D τ − f are of the same order ofmagnitude as the EVP stress acting on the cell face, and this can be achieved using a characteristic viscosity η a defined as the ratio of a typical stress to a typical rate of strain for the given problem.The present technique can also be interpreted as a “both-sides diffusion” technique [52] where themomentum equation discretised by the FVM is not (6) but an equivalent equation where the same diffusionterm ∇ · ( η a ∇ u ) has been subtracted from both sides: ∂ ( ρu ) ∂t + ∇ · ( ρuu ) − ∇ · ( η a ∇ u ) = −∇ p + ∇ · τ − ∇ · ( η a ∇ u ) (44)The left-hand side such term is not discretised in the same way as the right-hand side term; in particular,the component D τ + f in (41) comes from the discretisation of the left-hand side term, whereas the component D τ − f comes from the discretisation of the right-hand side term. The other components of these diffusionterms are discretised in exactly the same way for both of them and cancel out, leaving only D τ + f and D τ − f . Since both of these diffusion terms are discretised by central differences they contribute O (1) to thetruncation error [44] but O ( h ) to the discretisation error. In particular, consider that Eq. (44) correspondsto Eq. (23), and its FVM discretisation with the present schemes leads to the algebraic system (24), withthe dependent variable φ ≡ u being the velocity. If φ e is the exact velocity (the solution of Eq. (23)) and φ is the approximate velocity obtained by solving the system (24), then the discretisation error is ε = φ e − φ .The latter is related to the truncation error α by F (cid:48) h ( φ h ) · ε h ≈ − α h , where F (cid:48) h ( φ h ) is the Jacobian matrix ofthe operator F h evaluated at the approximate solution φ h (see e.g. [53] for a derivation). Let us focus on the ∇ · ( η a ∇ u ) term of the left-hand side of (44), which is discretised with a compact stencil giving rise to the D τ + f term of the scheme (41) (a similar analysis of the term on the right-hand side can be made). Becausethe diffusion operator is linear, its contribution to the truncation error F (cid:48) h ( φ h ) · ε h is equal to its contributionto F h ( ε h ), i.e. it can be calculated by replacing the velocities by their discretisation errors in (42) (we neglectthe rest of the contributions, which cancel out with the corresponding ones of the diffusion term of the rightside of Eq. (44)). So, plugging ε into Eq. (42) instead of u i we get a contribution to the truncation error ofcell P of ( s f / Ω P ) η a ( ε N f − ε P ) / (cid:107) N f − P (cid:107) = O ( h − )( ε N f − ε P ). If our method is second-order accurate then ε = O ( h ). On unstructured grids, the variation of discretisation error in space has a random component[53] so that ε N f − ε P = O ( h ) as well, which results in a O (1) contribution to the truncation error. Onsmooth structured grids the discretisation error varies smoothly [53] so that ε N f − ε P = O ( h ) and each facecontributes O ( h ) to the truncation error (on such grids there is also a cancellation between opposite facesfor each cell which reduces the net contribution to the truncation error to O ( h )).11nother benefit that comes from the incorporation of the stabilisation terms in the scheme (41) concernsthe algebraic solver, SIMPLE [54]. Within each SIMPLE iteration a succession of linearised systems ofequations are solved, one for each dependent variable. The systems for the velocity components u i are derivedfrom linearisation of the (discretised) momentum equation (6). The latter, unlike in (generalised) Newtonianflows, contains no diffusion terms and the velocity appears directly only in the convection terms of the left-hand side. In low Reynolds number flows, these terms play only a minor role and the momentum balanceinvolves mainly the pressure and stress. Constructing the linear systems for the velocity components onlyfrom these inertial terms would lead to bad convergence of the SIMPLE algorithm. Nevertheless, velocityplays an important indirect role in the momentum equation through the EVP stress tensor, to which itis related via the constitutive equation (7). Therefore, either the momentum and constitutive equationshave to be solved in a coupled manner, or the effect of velocity on the stress tensor has to somehow bedirectly quantified in the momentum equation, which is precisely what the stabilisation term D τ + f achieves;in particular, the diffusion term on the left-hand side of Eq. (44) contributes to the matrix of coefficients ofthe linear systems for u i , making this matrix very similar to those for (generalised) Newtonian flows.Finally, the momentum balance for cells that lie at the wall boundaries involves the pressure and thestress tensor at these boundaries. The pressure is linearly extrapolated from the interior, as is a commonpractice [55]. For stress two possible options are linear extrapolation, as for pressure, and the imposition ofa zero-gradient condition which is roughly equivalent to setting the stress at the boundary face equal to thatat the owner cell centre. In [56], the former approach was found, as expected, to be more accurate, whilethe latter was found to be more robust. In the present work we employed either linear extrapolation or amodified extrapolation which is akin to the scheme (41). So, if b is a boundary face of cell P , the linearlyextrapolated value of stress at its centre is calculated as¯ τ c b = τ P + ∇ h τ P · ( c b − P ) (45)where the least-squares stress gradient ∇ h τ P is calculated with q = 1 and using only the stress values atthe centres of P and its neighbour cells. This scheme is used also for extrapolation of pressure, but when itcomes to stress we can go a step further: the force due to the EVP stress tensor on boundary face b can becalculated by a scheme that has precisely the form (41) only that the interpolated stress ¯ τ c f (27) is replacedby the extrapolated value ¯ τ c b (45) and the D terms are replaced by D τ + b,i ≡ η a u i,c b − u i,P (cid:107) c b − P (cid:107) (46) D τ − b,i ≡ η a ¯ ∇ qh u i,m b · d b (47)where d b is the unit vector pointing from P to c b and the gradient in (47) is extrapolated to point m b ≡ ( P + c b ) / c b replaced by m b . The terms (46) – (47) do not have an apparentstabilising effect, and in fact our experience has shown that they sometimes can cause convergence difficultiesto the algebraic solver. Such is the case with wall slip examined in Sec. 6.4, where we had to set these termsto zero and use simple linear extrapolation of stresses at the walls. Nevertheless, we also found that theseterms can increase the accuracy of the solution (see Sec. 5) and therefore we employed them wheneverpossible (they were used for all the main results except those of Sec. 6.4). Again, by integrating the SHB equation (7) – (8) and applying the divergence theorem we obtain (cid:90) Ω P ∂τ∂t dΩ + (cid:88) f (cid:90) s f n · uτ d s = G (cid:90) Ω P ˙ γ dΩ − G (cid:90) Ω P (cid:18) max(0 , τ d − τ y ) k (cid:19) n τ d τ dΩ + (cid:90) Ω P (cid:0) ( ∇ u ) T · τ + τ · ∇ u (cid:1) dΩ (48)The surface integral on the left-hand side derives from application of the divergence theorem to the secondterm on the right-hand side of Eq. (8), using the identity u ·∇ τ = ∇· ( uτ ) − ( ∇· u ) τ and the incompressibilitycondition ∇ · u = 0. Equation (48) has the form of a conservation equation for stress. In the left-hand side,12 U Dc = c (cid:48) d Figure 5:
On a uniform structured grid, CUBISTA interpolates the value at face centre c from the values at thedownwind ( D ), upwind ( P ), and farther-upwind ( U ) cells (assuming the mass flux is directed from cell P to cell D ). there is the rate of change of stress in the volume plus the outflux of stress, which equal the rate of changeof stress in a fixed mass of fluid moving with the flow. In the right-hand side, there are three source termswhich express, respectively, stress generation (the ˙ γ term), relaxation (the viscous term), and transformation(the upper convective derivative term). Equation (48) is then discretised as: ∂τ∂t (cid:12)(cid:12)(cid:12)(cid:12) aP Ω P + 1 ρ (cid:88) f ˙ M f ¯ τ C c f = G Ω P ˙ γ P − G Ω P (cid:18) max(0 , τ d,P − τ y ) k (cid:19) n τ d,P τ P + Ω P (cid:0) ( ∇ qh u ) T P · ¯ τ P + ¯ τ P · ∇ qh u P (cid:1) (49)where ∂/∂t | aP is the approximate time derivative at P (Sec. 3.5), ˙ M f is the mass flux (30), ˙ γ P = ∇ qh u P +( ∇ qh u ) T P is the discrete value of ˙ γ at P , and τ d,P is the value of τ d at P , calculated from τ P via Eq. (9).In the convection term of Eq. (49) (the second term on the left-hand side), ¯ τ C c f is the value of τ inter-polated at the face centre c f , but not by the linear interpolation (27). Since the constitutive equation lacksdiffusion terms, there is a danger of spurious stress oscillations, and a common preventive measure is to useso-called high resolution schemes [57] for the discretisation of the convection term. In viscoelastic flows, theCUBISTA scheme [58] has proved to be effective and is widely adopted (e.g. [59–61, 29]). In the presentwork, we adapted it for use on unstructured grids as follows. To account for skewness, a scheme similar to(27) is used, but the value at point c (cid:48) f is interpolated with CUBISTA:¯ φ C c f = ¯ φ C c (cid:48) f + ¯ ∇ qh φ c (cid:48) f · ( c f − c (cid:48) f ) (50)Note that the second term of the right-hand side is exactly the same as for the scheme (27) (the gradientis interpolated linearly). The first term of the right-hand side is the value at point c (cid:48) f interpolated viaCUBISTA. CUBISTA is a high-resolution scheme based on the Normalised Variable Formulation [62], andon structured grids such schemes interpolate the value at a face from the values at two upwind cell centresand one downwind cell centre; these are the two cells on either side of the face, labelled P and D in Fig.5, and a farther cell U . But on an unstructured grid the farther-upwind cell U is not straightforward (oreven not possible) to identify. To overcome this hurdle, we follow an approach [63, 57] which is based onthe observation that, on a uniform structured grid such as that shown in Fig. 5, it holds that φ U = φ D − ∇ qh φ P · ( D − U ) (51)because the least-squares gradient is such that ∇ qh φ P · ( D − U ) = φ D − φ U . On unstructured grids, following[63, 57] we define the location U ≡ P − ( D − P ) which lies at the diametrically opposite position of D relativeto P ; since the grid is unstructured, it is unlikely that U coincides with an actual cell centre, but still wecan use Eq. (51) to estimate a value there, φ U . Information from the upwind side of cell P is implicitlyincorporated into φ U through the gradient ∇ qh φ P , and if the scheme is used on a uniform structured gridthen it automatically retrieves the value at the centroid of the actual cell U .So now we have a set of three co-linear and equidistant points, U , P and D , and three correspondingvalues, φ U , φ P and φ D , and CUBISTA can be employed to calculate the value at point c (cid:48) f , which lies on the13ame line as these three points. Defining the normalised variables as ξ = (cid:107) c (cid:48) − P (cid:107)(cid:107) D − P (cid:107) (52)˜ φ P = φ P − φ U φ D − φ U (53)˜ φ c (cid:48) f ( ξ, ˜ φ P ) = ¯ φ C c (cid:48) f − φ U φ D − φ U (54)(note that ξ is just the linear interpolation factor multiplying φ N f in Eq. (28)), CUBISTA blends QuadraticUpwind Interpolation (QUICK) and first-order upwinding (UDS) as follows:˜ φ c (cid:48) f ( ξ, ˜ φ P ) = ˜ φ P if ˜ φ P ≤ (1 + ξ )(3 + ξ ) ˜ φ P if 0 < ˜ φ P < (Transition 1)(1 − ξ ) ˜ φ P + ξ (1 + ξ ) if < ˜ φ P ≤ (QUICK)(1 − ξ ) ˜ φ P + ξ (2 − ξ ) if < ˜ φ P ≤ φ P if ˜ φ P > φ c (cid:48) f the un-normalised value ¯ φ C c (cid:48) f can be recovered via (54) for use in Eq. (50). However, the schemecan be conveniently reformulated in terms of un-normalised φ values using the median function [64] as:¯ φ C ∗ c (cid:48) f = median (cid:18) φ P , φ U + (1 + ξ )(3 + ξ )3 ( φ P − φ U ) , φ D − (1 − ξ ) ( φ D − φ P ) (cid:19) (56)¯ φ C c (cid:48) f = median (cid:18) φ P , ¯ φ C ∗ c (cid:48) f , φ P + φ D − φ U ξ + 12 ( φ D − φ P + φ U ) ξ (cid:19) (57)which also automatically sets ¯ φ C c (cid:48) f = φ P (UDS) when φ D = φ U . The scheme (55) is based on QUICK,but switches to UDS when there is a local minimum or maximum at P ( ˜ φ P < φ P > φ P ∈ [3 / , / φ U , φ P and φ D is not far from linear, is considered to be “safe” and QUICK is applied unreservedly there.The upper bound of this region, which for non-uniform structured grids is given to be a ξ -dependent valuein [58], is chosen here to be the fixed number ˜ φ P = 3 / U , P and D is monotonic (noovershoots / undershoots). In terms of normalised variables, this condition requires that ˜ φ c (cid:48) f ( ξ, ˜ φ P ), asgiven by the QUICK branch of Eq. (55), increases monotonically towards the value 1 as ξ →
1, which isensured by the condition ∂ ˜ φ c (cid:48) f ( ξ, ˜ φ P ) /∂ξ | ξ =1 ≥ ⇒ ˜ φ P ≤ /
4. We also note that, if CUBISTA achieves itsgoal of producing smooth φ distributions, then grid refinement, which brings the points U , P and D closertogether, will cause φ to vary more linearly in the neighbourhood of the three points, i.e. ˜ φ P → . φ .Finally, one can notice an interpolated value ¯ τ P in the stress transformation terms in the right-handside of Eq. (49). An obvious choice is to set ¯ τ P = τ P , but it was found in [65] that using instead a weightedaverage of the stress at P and its neighbours can mitigate the high- Wi problem. In the original scheme [65]the weighting was based on the mass fluxes through the cell’s faces. However, we noticed that this causes anoticeable accuracy degradation and instead used the following scheme (written here for 2D problems):( ∇ qh u ) T P · ¯ τ P + ¯ τ P · ∇ qh u P = τ ,P ∂u ∂x + 2 τ ,P ∂u ∂x τ ,P ∂u ∂x + τ ,P ( ∂u ∂x + ∂u ∂x ) + τ ,P ∂u ∂x τ ,P ∂u ∂x + 2 τ ,P ∂u ∂x (1,1) component(1,2) component(2,2) component (58) It can be implemented as median( a, b, c ) = max( min( a, b ) , min(max( a, b ) , c ) ). ∇ qh u P and τ ij,P = 1 D Ω P F (cid:88) f =1 ¯ τ C ij,c f s f ( c f − P ) · n f (59)where D = 2 or 3, for 2D or 3D problems, respectively. The interpolation scheme (59) is 2nd-order accurate(Appendix D). In the (1,2) component of (58), the middle term, although equal to zero in the limit of infinitegrid fineness due to the continuity equation, has been retained for numerical stability. Compared to simplysetting τ ij,P = τ ij,P , this scheme was found, in the numerical tests of Sec. 5, to allow an increase of about40% in the maximum solvable Wi number and to provide a slight increase of accuracy. At wall boundaries itwas found that using linearly extrapolated stress values in (59) leads to convergence difficulties. Therefore,wall boundary values were omitted in the calculation, weighting the values at the remaining faces by their s f ( c f − P ) · n f factors, i.e. by replacing D Ω P in the denominator of (59) by D Ω (cid:48) P = (cid:88) f / ∈{ W P } s f (cid:0) c f − P (cid:1) · n f = D Ω P − (cid:88) f ∈{ W P } s f (cid:0) c f − P (cid:1) · n f (60)where { W P } is the set of wall boundary faces of cell P . This is a 1st-order accurate interpolation.In the simulations of Sec. 6 we did not encounter the high- Wi problem, but should it arise one cantransform the constitutive equation in terms of a more weakly varying function of the stress tensor such asthe logarithm of the conformation tensor [66, 67]. This has been applied to the Saramito model in [68]. The approximate time derivatives in the momentum and constitutive equations are calculated with a2nd-order accurate, three-time level implicit scheme with variable time step. Suppose we enter time step i and let φ iP be the current, yet unknown, value of φ at cell P . Fitting a quadratic Lagrange polynomialbetween the three points ( t i , φ iP ), ( t i − , φ i − P ) and ( t i − , φ i − P ) and differentiating it at t = t i gives thefollowing approximate time derivative:d φ d t (cid:12)(cid:12)(cid:12)(cid:12) aP ( t i ) = l (cid:48) i,i ( t ) φ iP + l (cid:48) i,i − ( t ) φ i − P + l (cid:48) i,i − ( t ) φ i − P (61)where l (cid:48) i,i ( t ) = 2 t i − t i − − t i − ( t i − t i − )( t i − t i − ) , l (cid:48) i,i − ( t ) = 2 t i − t i − t i − ( t i − − t i )( t i − − t i − ) , l (cid:48) i,i − ( t ) = 2 t i − t i − t i − ( t i − − t i )( t i − − t i − ) (62)All other terms of the governing equations are evaluated at the current time (the scheme is fully implicit).After solving the equations to obtain φ iP and deciding (as will be described shortly) the next time stepsize ∆ t i +1 = t i +1 − t i , to facilitate the solution at the new time t i +1 we can make a prediction by evaluatingthe aforementioned Lagrange polynomial at t = t i +1 to obtain a provisional value φ i +1 ∗ P . This value servesas an initial guess for solving the equations at the new time step i + 1, offering significant acceleration.We perform such a prediction even for pressure, whose time derivative does not appear in the governingequations, as this was found to accelerate the solution of the continuity equation.The subsequent solution at time t i +1 will give us a value φ i +1 P (cid:54) = φ i +1 ∗ P , in general. In order to obtain anestimate of the O ( h ) truncation error of scheme (61) we can augment our set of three points by the additionof the point ( t i +1 , φ i +1 P ) and fit a cubic polynomial through the four points; the difference between ∂φ/∂t | aP and ∂φ/∂t | cP , the derivatives predicted at time t i by the quadratic and cubic polynomials, respectively, isan approximation of the truncation error associated with ∂φ/∂t | aP . Omitting the details, it turns out that ∂φ∂t (cid:12)(cid:12)(cid:12)(cid:12) cP ( t i ) − ∂φ∂t (cid:12)(cid:12)(cid:12)(cid:12) aP ( t i ) = c φ i +1 P − φ i +1 ∗ P ∆ t i +1 (63)where the factor c = O (1) depends on the relative magnitudes of ∆ t i − , ∆ t i and ∆ t i +1 ( c = 1/3 for∆ t i − = ∆ t i = ∆ t i +1 ). This allows us to keep an approximately constant level of truncation error by15djusting the time step sizes so as to keep the right-hand side of Eq. (63) constant. We set the followinggoal, which accounts for all grid cells at once: g it ≡ (cid:118)(cid:117)(cid:117)(cid:116) (cid:88) P Ω P (cid:18) φ iP − φ i ∗ P ∆ t i (cid:19) = ˜ g t Φ T (64)I.e. we want g it , the L norm of ( φ iP − φ i ∗ P ) / ∆ t i , for any time step i , to equal a pre-selected non-dimensionaltarget value ˜ g t times Φ /T where Φ is either U (for the momentum equation, where φ ≡ u j ) or S (for theconstitutive equation, where φ ≡ τ jk ). Suppose then that at some time step i we are somewhat off target,i.e. Eq. (64) is not satisfied; we can try to amend this at the new time step i +1 by noting that the truncationerror, and the associated metric g t , are of order O (∆ t ), which means that g i +1 t /g it = (∆ t i +1 / ∆ t i ) . Thisrelation allows us to choose ∆ t i +1 so that g i +1 t will equal ˜ g t Φ /T (Eq. (64)):∆ t i +1 = r it ∆ t i where r it = (cid:113) ˜ g t / ˜ g it , ˜ g it = g it Φ /T (65)In practice, we calculate the adjustment ratios r it for all variables whose time derivatives appear in thegoverning equations (one per velocity component and per stress component) and choose the smallest amongthem. We also limit the allowable values of r it to a range r it ∈ [ r min t , r max t ], and also the overall time stepsize to ∆ t i ∈ [∆ t min , ∆ t max ]. Typical values used for these parameters in the simulations of Sec. 6 are ˜ g t =10 − , r min t = 0 . r max t = 1 . t min =5 × − s, ∆ t max = 2 × − s (these may vary slightly betweensimulations). The scheme automatically chooses small time steps when the flow evolves rapidly, such asduring the early stages of the flow, and large time steps when it evolves slowly, such as when the flow isnearing its steady state. A similar scheme is used in [69].
4. Solution of the algebraic system
Discretisation results in one large nonlinear algebraic system per time step. These systems are solvedusing the SIMPLE algorithm, with multigrid acceleration, where the algebraic equations are arranged intogroups, one for each momentum and constitutive equation component, each of which through linearisationproduces a linear system for one of the dependent variables. The algorithm is applied in the same manneras for viscoelastic flows [26, 28]; in each SIMPLE iteration, we solve successively: the linear systems foreach stress component, the linear systems for the velocity components, and the pressure correction systemto enforce continuity. The systems are solved with preconditioned conjugate gradient (pressure correction)or GMRES (all other variables) solvers.In the velocity systems, derived from the components of the momentum eq. (40), the matrix of coef-ficients is constructed via a contribution from the time derivative, a UDS discretisation of the convectiveterm, and the D τ + f,i (Eq. (42)) part of the EVP force F τf . The remaining terms are evaluated from theircurrently estimated values and moved to the right-hand side vector, as is the difference between the UDSand CDS convection schemes (deferred correction). For the stress systems we follow the established practicein viscoelastic flows of constructing the matrix of coefficients with contributions only from the terms of theleft-hand side of Eq. (49), breaking the CUBISTA flux into a UDS component and a remainder which ismoved to the right-hand side vector (deferred correction). This poor representation of the constitutive equa-tion by the matrix of coefficients may be responsible for the observed convergence difficulties of SIMPLE atlarge time steps; for small time steps the role of the time derivative in the constitutive equation becomesdominant and convergence is fast.It was noticed that if the residual reduction within each time step is not at least a couple orders ofmagnitude, then the temporal prediction step (Sec. 3.5) may not produce good enough initial guesses andthe solution may not be smooth in time. To avoid this possibility we set a maximum effort per time step ofabout 20 single-grid SIMPLE iterations followed by about 5 W(4,4) multigrid cycles [70, 24], and if furtheriterations are needed to accomplish the required residual reduction then the time step is reduced by a factorof r min t (Sec. 3.5). Setting ∆ t max to an appropriate value avoids the need for this.16 able 1: ˜ x − and ˜ y − coordinates (˜ x c , ˜ y c ) of the centre of the main vortex for Oldroyd-B flow in a square cavity of side L , and associated value of the streamfunction there, ˜ ψ c . The top wall moves in the positive x − direction with variablevelocity u ( x ) = 16 U ( x/L ) (1 − x/L ) . The coordinates are normalised by the cavity side L and the streamfunctionby U L . The flow is steady, the Reynolds number is zero, and the solvent and polymeric viscosities are equal ( κ = k ). Wi = 0 . Wi = 1 . x c ˜ y c ˜ ψ c ˜ x c ˜ y c ˜ ψ c Pan et al. [71] 0 .
469 0 . − . .
439 0 . − . .
429 0 . − . .
467 0 . − . .
434 0 . − . .
468 0 .
798 0 .
430 0 . .
468 0 . − . .
434 0 . − .
5. Validation and testing of the method in Oldroyd-B flow
For τ y = 0 and n = 1, the SHB model reduces to the Oldroyd-B viscoelastic model, for which benchmarkresults for square lid-driven cavity flow can be found in the literature. We apply our method to this problemto validate it and to evaluate alternative options for the FVM components. The problem is solved as steady-state, omitting the time derivatives from the governing equations. Table 1 lists the computed location andstrength of the main vortex, along with values reported in the literature, with which there is very goodagreement.To estimate discretisation errors, we obtained accurate estimates of the Wi = 0 . Wi = 1 . × × (cid:15) u ≡ (cid:48) (cid:88) P ∈ Ω (cid:48) (cid:107) u eP − u P (cid:107) / (cid:107) u e (cid:107) avg (66) (cid:15) τ ≡ (cid:48) (cid:88) P ∈ Ω (cid:48) (cid:107) τ eP − τ P (cid:107) / (cid:107) τ e (cid:107) avg (67)where the superscript e denotes the “exact” solution, calculated by Richardson extrapolation. The sum-mation is over all cells P whose centres lie inside the area Ω (cid:48) = [0 . L, . L ] × [0 . L, . L ] ( (cid:48) is thetotal number of such cells), L being the cavity length, i.e. a strip of width 0 . L touching the boundary isexcluded, because the stress magnitude ( τ xx in particular) appears to grow exponentially over part of thetop boundary and this inflates the discretisation errors there. The vector and tensor magnitudes in (66)and (67) are (cid:107) u (cid:107) = ( (cid:80) i u i ) / and (cid:107) τ (cid:107) = ( (cid:80) ij τ ij ) / . The definitions (66) and (67) incorporate a normali-sation of the errors by the average velocity and stress magnitudes over Ω (cid:48) : (cid:107) u e (cid:107) avg ≡ (1 / Ω (cid:48) ) (cid:82) Ω (cid:48) (cid:107) u e (cid:107) dΩ and (cid:107) τ e (cid:107) avg ≡ (1 / Ω (cid:48) ) (cid:82) Ω (cid:48) (cid:107) τ e (cid:107) dΩ. In Fig. 7 the discretisation errors are plotted with respect to the average gridspacing h = (Ω / / which equals L/N for a N × N grid of any of the kinds shown in Fig. 6.To ensure that our method is applicable to unstructured grids, although the geometry favours discreti-sation by uniform (Fig. 6a, “CU”) or non-uniform (Fig. 6b, “CN”) Cartesian grids, we also employed gridsobtained from the CU ones by randomly perturbing their vertices as described in [44]. In particular, froma N × N CU grid (Fig. 6a) a N × N distorted grid is constructed (Fig. 6c) by moving each interior node( x, y ) to a location ( x (cid:48) , y (cid:48) ) = ( x + δx, y + δy ) where δx and δy are randomly selected for each node underthe restriction that | δx | , | δy | < h/
4. This restriction ensures that all resulting grid cells are simple convexquadrilaterals. This procedure produces a series of grids whose skewness, unevenness and non-orthogonalitydo not diminish with refinement [44], as is typical of unstructured grids. Checking the convergence of themethod on this sort of grids is important because in [44] it was shown that there are popular FVM dis-cretisations, widely regarded as second-order accurate, which actually do not converge to the exact solution17 a) Cartesian uniform (CU) (b)
Cartesian non-uniform (CN) (c)
Distorted (D)
Figure 6: ×
32 grids of different kinds, used for the Oldroyd-B benchmark problems. with refinement on unstructured grids. So, we employed three series of grids, one for each of the grid typesshown in Fig. 6, having 32 ×
32, 64 ×
64, 128 × ×
256 and 512 ×
512 cells. The CN grids wereconstructed as follows: in the 512 ×
512 such grid, the grid spacing at the walls equals L/ ×
256 grid, and so on for coarser grids.The numerical tests led to the following observations. Concerning the oscillations issue, the stabilisationstrategies achieved their goal as all dependent variables varied smoothly across the domain without spuriousoscillations; Fig. 8 shows smooth pressure and τ fields and streamlines obtained on the 512 ×
512 distortedgrid. Concerning the accuracy of the method, Fig. 7 shows convergence to the exact solution with gridrefinement on any type of grid, including the randomly distorted type. The order of accuracy is close to 1on coarse grids and increases towards 2 (the design value) as the grid is refined, but for the Wi = 1 case(Figs. 7b, 7e) it is still quite far from 2 even on the finest (512 × n ≥
1, but not for n < Wi <
1, which, combined with the fact thattheir behaviour is described by n <
1, avoids the “high- Wi ” problem. Nevertheless, in the tests of thepresent Section we also tried the log-conformation technique [66, 67] which is implemented as an optionin our code [29], and, interestingly, Fig. 7 shows that it can be beneficial even at low Wi as it providesenhanced accuracy. Interpolation (59) also appears to improve the accuracy compared to simply setting τ ij,P = τ ij,P , although rate of convergence of the latter method appears to be higher for Wi = 0 . Wi compared to setting τ ij,P = τ ij,P (up to Wi ≈ . Wi = 1). For the main results of Sec. 6 the scheme (59) was selected, given also that it is cheaper thanthe log-conformation method.The scheme (26) performs poorly in terms of accuracy (Figs. 7a,7d and 7b,7e, dash-dot lines) comparedto the simpler scheme (25). Concerning the latter, it is noteworthy that the larger the value of η a (adjustedby replacing S with larger / smaller values in (25)) the less accurate the results (Figs. 7a,7d, lines withlong or short dashes). Thus, stabilisation should be exercised with parsimony. Henceforth we abandon thescheme (26) and employ scheme (25).Fig. 7 shows that packing the grid lines close to the walls – i.e. using the CN grids – achieves a verysignificant increase of accuracy. This is likely related to the aforementioned singular behaviour near thelid. As expected, the discretisation errors are largest on the distorted grids (Figs. 7c,7f, red lines), but notby a great margin compared to the uniform Cartesian grids. The rates of error convergence towards zeroare similar on all grid types. Another disadvantage of the distorted grids is that the high- Wi problem isintensified: for Wi = 1, we only managed to obtain a solution up to grid 64 ×
64 with the scheme (25),and up to grid 128 ×
128 with the scheme (26). Finally, we note that Figs. 7c,7f show that the stress18 a) Wi = 0 .
5, CU grids: (cid:15) u (b) Wi = 1 . (cid:15) u (c) Wi = 0 .
5; CU,CN,D grids: (cid:15) u (d) Wi = 0 .
5, CU grids: (cid:15) τ (e) Wi = 1 . (cid:15) τ (f ) Wi = 0 .
5; CU,CN,D grids: (cid:15) τ Figure 7:
Velocity (top row) and stress (bottom row) discretisation errors (Eqs. (66) and (67)) of Oldroyd-B flow asa function of the grid spacing h , for various cases. Unless otherwise stated, the ∇ . h gradient is used and the D ’s atboundary faces are given by (46) – (47). Figs. (a), (d): Wi = 0 .
5, CU grids (Fig. 6a). Green, red and blue linescorrespond to setting τ ij,P = τ ij,P in (58), using the log-conformation technique, and defining τ ij,P by (59) in (58),respectively. In the latter case, the solid and dash-dot lines correspond to the use of Eqs. (25) and (26), respectively.Dashed lines with short or long dashes correspond, respectively, to replacing S by S/ S in (25). Figs. (b), (e): as for Figs. (a), (d), but for Wi = 1 .
0. The dashed line now corresponds to CN grids (Fig. 6b), using (25) and (12).
Figs. (c), (f): Wi = 0 .
5, use of (25) and (12). Green, blue, and red lines: CU (Fig. 6a), CN (Fig. 6b), and D (Fig.6c) grids, respectively. Dashed lines correspond to use of ∇ . h and simple linear extrapolation of stresses to the walls(i.e. setting D τ + b,i = D τ − b,i = 0 instead of (46) – (47)). extrapolation scheme at the boundaries which uses (46) – (47), in combination with the gradient ∇ . h ,which retains 2nd-order accuracy at the boundaries unlike the more popular ∇ h [44] , leads to a noticeableaccuracy improvement compared to simple linear extrapolation of stresses.
6. EVP flow in a lid-driven cavity
We now turn our attention to EVP flow in a lid-driven cavity. The SHB model parameters were chosenso as to represent the behaviour of Carbopol. Carbopol gels are used very frequently as prototypicalviscoplastic fluids in experimental studies [11, 74]. In [33], the SHB model was fitted to experimental datafor a Carbopol gel of concentration 0.2% in weight, and we rounded those parameters to arrive at our ownchosen parameters, listed in Table 2. This material has a yield strain (Eq. (22)) of γ y = 0 . L = 0.1 m and the flow is driven by horizontal motion of the top wall (lid) towardsthe right. We employed a grid of 384 ×
384 cells, of uniform size in the x direction but packed near the lid sothat the vertical size of the cells touching the lid is L/ This is true only for variables with Dirichlet boundary conditions, i.e. the velocity components. For pressure and stress, ∇ . h also deteriorates to first-order accuracy. Hence we use the more standard gradient ∇ h in the extrapolation scheme (45). a) Pressure (b)
Stress τ (c) Stress τ Figure 8:
The Oldroyd-B, Wi = 0 . ×
512 distorted grid. (a) Pressure contours(colour, δp = 4 S with S given by Eq. (12)) and streamlines (white lines). (b) Contours of τ ( δτ = 2 S ) andstreamlines. (c) Close-up view of (b) near the upper-right corner. Table 2:
Properties of the fluid used in the EVP lid-driven cavity simulations. ρ / m τ y
70 Pa G
400 Pa k
20 Pa s n n L/ ×
48 grid of similar packing is shown in Fig. 9.
We start with a case where the enclosed material is initially at rest and fully relaxed ( τ = 0 everywhere).Starting from rest, the lid accelerates towards the right until it reaches a maximum velocity of U = 0.1 m / sat time T = L/U = 1 s. The lid velocity remains constant thereafter: u lid = U sin (cid:18) min( t, T ) T π (cid:19) (68)The dimensionless numbers for this case are listed in the “ U = 0.100 m / s” column of Table 3.Figure 10 shows the evolution in time of the normalised kinetic energy of the fluid and of the normalisedaverage absolute value of the trace of the stress tensor, calculated asN . K . E . = 1 U Ω (cid:88) P Ω P (cid:107) u P (cid:107) (69) Table 3:
Values of the dimensionless numbers, and of the relaxation time (21), at different lid velocities. U [m / s] 0.025 0.100 0.400 Bn (cid:48) Bn Wi Re λ [s] 0.815 0.255 0.066 igure 9: A coarse, 48 ×
48 cell grid, showing the packing of the cells near the lid (a coarse grid is shown for clarity;the EVP lid-driven cavity simulations were performed on a similar grid of 384 ×
384 cells). (a)
N.K.E. (b) tr(˜ τ ) avg Figure 10: (a) Normalised kinetic energy (Eq. (69)) and (b) normalised average trace of τ (Eq. (70)), versus time,for different lid velocities. In (a) time is normalised by L/U . tr(˜ τ ) avg = 1 S Ω (cid:88) P Ω P | tr( τ P ) | (70)where Ω = L is the volume of the cavity. It is immediately obvious from Fig. 10 that the kinetic energyassumes a constant value very quickly, but the average stress trace evolves much more slowly. To investigatethis, the simulation was carried on until a time of t = 210 s.The flow field is visualised in Fig. 11. It resembles that for pure viscoplastic flow [24, 34], in that thereare two unyielded zones ( τ d < τ y ), one at the bottom of the cavity touching the walls, containing stationaryfluid, and one near the lid which is rotating with the flow and does not touch the walls, called a plug zone.Whether the material at a point is in a yielded or unyielded state is, of course, determined by whether τ d islarger or smaller, respectively, than the yield stress τ y there (Eq. (7)). Figure 11b shows that the streamlinesdo not cross into the lower unyielded zone, which therefore always consists of the same material, but theycross into and out of the plug zone. Thus, at every instant in time, liquefied particles are entering the plugzone, solidifying upon entry, while other, solidified particles, are exiting the zone, liquefying upon exit.A striking feature of the flow evolution is the slowness with which the stationary unyielded zone tends toobtain its final shape. Figure 11a shows that although the plug zone has reached its steady state already from t = 30 s, the stationary zone continues to expand even at t = 210 s. The shape of the yield line delimitingthis zone is quite irregular. Furthermore, Fig. 11b shows that this zone is surrounded by an amount of fluid21 a) Unyielded areas (b) (cid:107) u (cid:107) [m / s] Figure 11: (a) Yielded (white) and unyielded (colour) regions for the base case (Sec. 6.1). Different colours denotethe extent of the unyielded regions at t = 30, 60, 90, 120, 150, 180 and 210 s (the lower unyielded region is expandingwith time). (b) Colour contours of the magnitude of the velocity vector at t = 210 s. The white lines are the yieldlines ( τ d = τ y ) and the dashed black lines are streamlines, both at t = 210 s. The streamlines are plotted at equalstreamfunction intervals. that is practically stationary and yet yielded. It is likely that this fluid will eventually become part of thestationary unyielded zone as t → ∞ , i.e. that the yield line will eventually lie somewhere close to the nearbystreamline drawn in Fig. 11b, which separates the fluid with near-zero velocity from that whose velocity islarger. However, to ascertain this the simulation would have to be prolonged to prohibitively long times;already the expansion of the lower unyielded zone from t = 180 s to t = 210 s, marked in black colour inFig. 11a, is very small.A related feature is that the magnitude of the deviatoric stress tensor, τ d , is very close to the yield stressthroughout the aforementioned near-zero velocity region into which the lower unyielded zone is expanding.Thus, this region appears as an almost completely flat surface in the three-dimensional plot of τ d in Fig. 12(the surface outlined in dashed line contains fluid where τ d is within ± τ y ). Due to its distinctivefeatures, we will refer to this region as a “transition zone”. In it, the fluid practically behaves as unyielded,although some of it can be formally yielded. The plug zone possesses no transition zone.In order to explain the transition zone behaviour, we consider the constitutive equation (7) – (8) in thecase that the fluid velocity is zero and τ d is slightly above τ y . It becomes: ∂τ∂t = − Gk n ( τ d − τ y ) n τ d τ ≡ − r ( τ d ) τ (71)Since the function r ( τ d ) assumes only positive values, the minus sign on the right-hand side of Eq. (71) drivesthe components of τ towards zero as time passes, with a rate proportional to r ( τ d ) (and to the magnitudeof those components themselves). Hence τ d decreases towards zero as well (Eq. (9)). However, the rate atwhich they decrease, r ( τ d ), diminishes as τ d → τ y and in fact r ( τ y ) = 0. Therefore, τ d actually convergestowards τ y and not to zero. This is what we observe happening inside the transition zone.It is useful to consider also a scalar version of Eq. (71) ( τ d ≡ τ in this case) :d τ d t = − Gk n ( τ − τ y ) n (72) Eq. (72) can be viewed as describing the behaviour of the mechanical system of Fig. 1 (with κ = 0) in the case that atsome stress τ = τ > τ y we pin the right end of the spring to a fixed location, so that the total length γ v + γ e remains constanthenceforth, and we leave the system to relax. The spring will try to recover its equilibrium length by contracting or expanding γ e , which will cause an opposite expansion or contraction of γ v , resisted by the viscous and friction elements. Once the tensionin the spring drops to τ y , the spring cannot overcome the friction τ y of the friction element and motion stops, without the springhaving attained its equilibrium length. igure 12:
3D contour plot of the base case results (Sec. 6.1) at t = 210 s. Both the colour contours and the plotheight ( z axis) represent τ d . In addition, the dashed line encloses an area where τ d ∈ [69 ,
71] Pa.
For n = 1, the solution to the above equation is τ − τ y = ( τ − τ y ) e − t/λ (cid:48) (73)where τ is the value of τ at t = 0 and λ (cid:48) ≡ k/G is a relaxation time. This behaviour is similar to that ofa Maxwell viscoelastic fluid, only that now τ decays exponentially towards τ y instead of towards zero. For n (cid:54) = 1, Eq. (72) can be written as ( λ (cid:48) ≡ k /n /G now has units of s / Pa − / n ):d( τ − τ y )d t = − λ (cid:48) ( τ − τ y ) (cid:124) (cid:123)(cid:122) (cid:125) rate of decay for n =1 ( τ − τ y ) − nn (74)so that the rate of decay of τ towards τ y equals that for n = 1 multiplied by a factor of ( τ − τ y ) (1 − n ) /n . If n <
1, as is the present choice but also the most common case, then (1 − n ) /n > τ → τ y , the extra factor ( τ − τ y ) (1 − n ) /n tends to zero, and the rate of decay becomes progressivelysmaller than that of the n = 1 case, eventually becoming infinitesimal compared to it. This behaviour isexplained physically by the fact that the relaxation time (Eq. (21)) is proportional to the fluid viscosity,which resists recovery from deformation (relaxation), and inversely proportional to the elastic modulus,which drives towards such recovery. For shear-thinning fluids ( n < increases as τ becomessmaller and, in the SHB and HB models, tends to infinity as τ → τ y . Thus, for these fluids the relaxationtime tends to infinity as τ → τ y . The opposite happens in the less common case of n > Re , because Re is very small (Table 3), it becomes apparent thatthe fluid particle accelerations (the left-hand side of that equation) are large and the velocity adjusts veryquickly to force changes (the inertial time scales are very small). Thus, the velocity should vary at the samerate as these forces, i.e. the velocity and stress variations should go hand-in-hand (the boundary conditionsare not a source of variation as the lid velocity remains constant for t > T = 1 s). However, in Fig. 10the velocity field appears to reach a steady state much faster than the stress field. This can be explainedby observing that the stress actually also evolves quickly over most of the domain, the relaxation time λ =0.255 s (Table 3) being quite small compared to the time scale of stress evolution in Fig. 10b, except in the In the language of the mechanical analogue of Fig. 1 this means that the damper component resisting the relaxation of thespring becomes stiffer as this relaxation proceeds. τ y , so its local slow evolution has a noticeableimpact on the average stress trace plotted in Fig. 10. The velocity does follow the stress evolution, but sinceit is almost zero in this zone its local variations caused by the local stress evolution have a negligible impacton the overall kinetic energy plotted in Fig. 10. We noticed that between t = 30 s and t = 210 s the changein the velocity components at any point in the domain is of the order of 10 − m / s, except in the transitionzone where it can be about five times larger.Figure 13 shows vertical profiles of most flow variables along the centreline ( x = L/
2, Figs. 13a and 13b)and close to the right wall ( x = 0 . L , Figs. 13c and 13d). Two profiles are drawn for each variable, oneat time t = 30 s and one at t = 210 s. These profiles are identical inside the yielded and plug zones, andonly deviate very slightly inside the bottom unyielded zone. Thus the steady state has largely been reachedalready at t = 30 s. The vertical line at x = 0 . L cuts through a significant width of the transition zone,which can be recognised by the vertical line segment where τ d ≈ τ y in Fig. 13c (approximately from y =0.041 m to y = 0.051 m). A close-up view of the variation of τ d in the transition zone is shown in the insetof the same figure, where one can see that as time passes, τ d decreases towards τ y throughout the zone.Interestingly, Fig. 13a shows that the normal stress component τ tends to vary discontinuously acrossthe boundary of the lower unyielded zone as time passes. The possibility of the SHB model producingsolutions with discontinuous stress components and/or velocity gradients was noted and discussed in [39, 22],where it was attributed to the combination of the upper convective derivative and the viscoplastic “max”term. Nevertheless, as noted in [39], stress discontinuities must be such that the components of the force ∇ · σ = ∇ · ( τ − pI ) remain bounded (discontinuities in τ that lead to infinite derivatives in ∇ · τ canexist if they are counterbalanced by opposite discontinuities in − pI ). Otherwise, at a point of discontinuitythere will result a finite (i.e. non-zero) force acting on an infinitesimal mass, producing infinite accelerationand making the velocity field discontinuous, but this violates the continuity equation for an incompressiblemedium. In our case though, the variation of τ in the x direction does not cause any force on the fluid(there is no ∂τ /∂x derivative in ∇ · τ ) and is therefore allowed to be discontinuous. Next we examine the flow driven by higher ( U = 0.4 m / s) and lower ( U = 0.025 m / s) lid velocities, whichaffects the dimensionless numbers as shown in Table 3. Lowering the velocity increases the Bingham numberand decreases the Weissenberg number, i.e. it accentuates the plastic character of the flow at the expenseof its elastic character. The Re values listed in Table 3 suggest that inertia effects should be negligible for U = 0.025 and 0.100 m / s, but should be slightly noticeable for U = 0.400 m / s. The Table shows also that λ increases as U decreases; at lower velocities the apparent viscosity η , Eq. (21), is higher and therefore thestresses relax more slowly. Thus, we expect the flow to evolve more slowly as U is reduced.The lid velocity is again increased gradually according to Eq. (68) with T = L/U . Figure 10a shows thatthe KE of the fluid builds up in an oscillatory manner, usually with an overshoot, due to elastic effects. Thelarger U is, the more prominent and persistent the oscillations. Figure 10b confirms that the stress evolutionis slower at lower U , as discussed above. The top row of Fig. 14 shows the corresponding near-steady-stateflow fields. Increasing the lid velocity leads to less unyielded material in the cavity, and the vortex has morefree space to move away from the lid (see the y c coordinate in Table 4). In each case there is a transitionzone between the bottom unyielded zone and the yielded material, with the former not having yet expandedthroughout the near-zero velocity region. The number and density of streamlines in Fig. 14 shows thatat higher Bn the flow is weaker and more confined to a thin layer below the lid while circulation is veryweak in the rest of the domain. This is reflected also in the normalised KE diagram 10a, and also in thenormalised vortex strengths listed in Table 4. That Table also shows that the vortex lies slightly to the leftof the centreline, as is typical of viscoelastic lid-driven cavity flows, although this shift is not as pronouncedas for the Oldroyd-B flows of Table 1. Interestingly, as U is increased the vortex centre moves towards thecentreline despite Wi increasing; this could be due to inertia effects that will be discussed in Sec. 6.3. Since the HB equation is used extensively to model viscoplastic flows, we compare its predictions againstthose of the SHB equation in order to get a feel of the error involved in neglecting elastic effects. The HB24 a) u , τ and τ d at x = L/ (b) p , τ and τ at x = L/ (c) u , τ and τ d at x = 0 . L . (d) p , τ and τ at x = 0 . L . Figure 13:
Profiles of some dependent variables along vertical lines, for the base case (Sec. 6.1): (a), (b) at thevertical centreline ( x = 0 . L ); (c), (d) close to the right wall ( x = 0 . L ). Dashed / continuous lines denote profiles at t = 30 s / 210 s, respectively. The vertical dash-dot lines in (a) and (c) mark the yield stress value τ y . The unyieldedand transition regions are shaded. The inset in (c) shows a close-up view of the τ d profile at y ∈ [0 . , . simulations were carried out by performing Papanastasiou regularisation [75], implemented as [76]: τ = η ˙ γ , η = τ ˙ γ ≈ ( τ y + k ˙ γ n )(1 − e − ˙ γ/(cid:15) )˙ γ (75)where the regularisation parameter (cid:15) determines the accuracy of the approximation. This method was usedfor simulating lid-driven cavity Bingham flows in [24, 34], where it was noticed that the equations becametoo stiff to solve for (cid:15) < / (cid:15) is progressively halved, solutions for (cid:15) = 1/128,000 were obtained. We performed steady-state simulations(since the steady-state of classic viscoplastic flows does not depend on the initial conditions) where the HBparameters τ y , k and n have the same values as for the SHB fluid (Table 2).The flowfields, depicted in the second row of Fig. 14, are much more symmetric than their SHB counter-parts. In HB lid-driven cavity flow, the only source of asymmetry is inertia. Inertial effects are imperceptible25 able 4: Dimensionless coordinates of the centre of the vortex (˜ x c , ˜ y c ) ≡ ( x c /L, y c /L ) and value of the streamfunctionthere ˜ ψ c ≡ ψ c / ( LU ), for various cases, at steady-state.SHB HB˜ x c ˜ y c ˜ ψ c ˜ x c ˜ y c ˜ ψ c U = 0.025 m / s 0 .
495 0 . − . .
500 0 . − . U = 0.100 m / s 0 .
495 0 . − . .
500 0 . − . U = 0.400 m / s 0 .
500 0 . − . .
505 0 . − . U = 0.100 m / s, slip 0 .
482 0 . − . .
500 0 . − . (a) SHB, U = 0 .
025 m / s, t = 90 s (b) SHB, U = 0 .
100 m / s, t = 210 s (c) SHB, U = 0 .
400 m / s, t = 60 s (d) HB, U = 0 .
025 m / s (e) HB, U = 0 .
100 m / s (f ) HB, U = 0 .
400 m / s Figure 14:
Top row: SHB flow snapshots at the times indicated, for different lid velocities. Bottom row: corre-sponding HB steady-state flowfields. The unyielded regions are shown shaded; the lines are instantaneous streamlines,plotted at streamfunction values ψ/ ( U L ) = 0.04, 0.08, 0.12 etc., plus a ψ = 5 × − m / s streamline just above thelower unyielded region. for U = 0.025 and 0.100 m / s (Figs. 14d and 14e) and only slightly noticeable for U = 0.400 m / s (Fig. 14d),which is in accord with the corresponding Re values listed in Table 3. So, for U = 0.4 m / s the vortex centreis shifted slightly towards the right (Table 4), whereas at the lower U ’s it is almost exactly on the centreline.Also, in Fig. 14f one can see that the upper unyielded (plug) region is somewhat stretched towards theleft. These features are opposite to those of the SHB case, where elasticity causes the vortex centre tomove towards the left and the plug region to stretch towards the right. The opposite effects of inertia andelasticity have been noted also in [36, 77, 29].Figure 14 and Table 4 show that the union of the lower unyielded and transition zones in SHB flow isslightly larger than the corresponding HB stationary regions, pushing the SHB vortex upwards and loweringits strength compared to the HB case. Figures 15a and 15b show that the SHB and HB velocity profiles26 a) U = 0 .
025 m / s (b) U = 0 .
400 m / s (c) τ /S Figure 15:
Profiles of u and τ d for U = 0.025 m / s (a) and U = 0.400 m / s (b), and profiles of τ /S for both U ’s (c).All profiles are along the vertical centreline ( x = L/ t = 90 s ( U = 0.025 m / s)and t = 60 s ( U = 0.400 m / s); dash-dot lines: SHB profiles at 30 s earlier; dashed lines: HB steady-state profiles. In(a) and (b) the vertical lines mark the yield stress. along the vertical centreline are very similar; the τ d profiles are somewhat more dissimilar but still not farapart, especially in the yielded regions. In the unyielded regions they cannot be expected to be similar dueto the stress indeterminacy of the HB model (the currently predicted HB stress field in the unyielded regionsis one of infinite possibilities, the one conforming to the Papanastasiou regularisation). Profiles of τ areplotted in Fig. 15c; for HB flow, due to the symmetric flowfield and τ being proportional to ∂u /∂x ,this stress component is nearly zero, except inside the plug region for U = 0.400 m / s, which is somewhatasymmetric (Fig. 14f). The SHB stresses, on the other hand, are significant due to elasticity, especially inthe higher U case which corresponds to higher Wi . Finally, we note that in Fig. 15 two SHB profiles areplotted for each lid velocity: at t = 60 and 90 s for U = 0.025 m / s, and at t = 30 and 60 s for U = 0.400m / s. These profiles are hardly distinguishable, indicating that the steady state for U = 0.025 m / s has beenpractically reached at t = 60 s, and for U = 0.400 m / s at t = 30 s.Next, we performed a couple of simulations similar to the base case, only that we artificially increasedthe elastic modulus to 10 and 100 times the original value listed in Table 2. This reduces the relaxation timesby factors of 10 and 100, respectively, so that the SHB material becomes more inelastic and its predictionsshould tend to match those of the HB model. Due to the smaller relaxation times, we also used smallertime steps: starting from an initial time step of 5 × − s in both cases, we allowed the time step toincrease up to ∆ t max = 1 . × − s and 10 − s in the G = 4,000 Pa and G = 40,000 Pa cases, respectively.The simulations were stopped at t = 50 s and 10 s respectively, and in the G = 40,000 Pa case the lidacceleration period (Eq. (68)) was reduced to 0.01 s. Figure 16 shows that the SHB and HB streamlinesconverge rapidly with increasing G , while Fig. 17 shows that convergence is not as rapid for the stresses.Of course, for τ d < τ y (top row of Fig. 17) we cannot expect the SHB and HB stresses to ever match,due to the aforementioned indeterminacy of the HB stress tensor in the unyielded regions. The yield lines(second row of Fig. 17) appear to converge, but at a slow pace. Despite the shorter duration of the G =40,000 Pa simulation, Fig. 17f shows that during this time the transition zone has evolved further than inthe lower G cases, due to the much shorter relaxation times. The same figure reveals some slight spuriousstress oscillations on the lower yield line (can also be seen in Fig. 17d). They could be due to the near-zerovelocities there, which complicates the application of upwinding within CUBISTA. Finally, the last row ofFig. 17 shows that within the yielded region the τ d predictions of the two models are quite close, even for G = 400 Pa. A persistent discrepancy between the two models at the sides of the plug region seen in Fig. 17imay be due to spatial or temporal discretisation errors, regularisation, etc.27 a) G = 400 Pa (b) G = 4 ,
000 Pa (c) G = 40 ,
000 Pa
Figure 16:
Continuous lines: SHB streamlines for various values of G , as indicated (the rest of the parameters areas in Table 2, and the lid velocity is U = 0.1 m / s). Dashed lines: corresponding HB streamlines, drawn at the samestreamfunction values as for the SHB model. Since the HB predictions are independent of G , the dashed lines do notchange between plots. The streamlines are drawn at fixed streamfunction intervals. Carbopol gels are quite slippery [11, 78, 74]; the previous, no-slip results can be considered to correspondto roughened walls, like those used in rheological measurements to avoid slip [33]. However, if the wallsare smooth then noticeable slip is expected. Viscoplastic and viscoelastic fluids can exhibit complex slipbehaviour, e.g. power-law slip [79, 80], pressure dependence [81], or “slip yield stress” [82, 83]. However,experiments with Carbopol in [78] (0.2% by weight) and [74] (0.08% by weight) showed nearly Navier slipbehaviour, u s = 5 . × − τ . w and u s = 2 . × − τ . w , respectively, where u s is the slip velocity (theleft-hand side of Eq. (11)) and τ w is the wall tangential stress (the term ( n · τ ) · s of the right-hand sideof Eq. (11)). In [74] it is conjectured that this may be due to the formation of a thin layer of Newtonianfluid (solvent) that separates the wall surface from the gel micro-particles. In the present section we imposea Navier slip condition, Eq. (11), on all walls, with a slip coefficient β = 5 × − m / Pa s. For this case,concerning stress extrapolation to the walls, we could not get SIMPLE to converge with the D coefficientsgiven by (46) – (47) and we set them equal to zero.Figure 18a shows the flowfield. A distinguishing feature is that the two unyielded zones touch the cavitywalls and yet contain moving material which is sliding over the wall. This feature is common also to theHB flowfield, shown in Fig. 18b, which is, again, much more symmetric. Figure 18c shows that SHB andHB contours of τ d are more similar for τ d > τ y . Due to slip, the flow induced by the lid is much weaker thanin the no-slip case, as can be seen from the lower kinetic energy in Fig. 10, the smaller vortex strength inTable 4, and the wider streamline spacing in Fig. 19b than in Fig. 19a. On the other hand, the streamlinesin Fig. 19b are more evenly spaced near the bottom of the cavity compared to Fig. 19a as slip allows thecirculation to extend all the way down to the bottom wall, where now there is no stationary unyielded zone.Consequently, there is also no transition zone. Indeed, comparing in Fig. 18a the yield lines at t = 30 s(dashed lines) and t = 60 s (the boundary of the shaded regions) one sees that there is very little change,and the bottom unyielded zone even slightly contracts with time. The absence of a transition zone can beexplained by examining the ˙ γ distributions of Fig. 19; whereas in the no-slip case (Fig. 19a) the velocitygradients are practically zero inside the bottom unyielded zone, in the slip case (Fig. 19b) they are non-zerothroughout the domain, thus excluding the transition zone situation where the constitutive equation reducesto Eq. (71). Figure 19 also shows that in the upper part of the cavity ˙ γ is much lower in the presence ofslip. Due to the shear-thinning nature of the fluid, this results in higher viscosities and associated relaxationtimes which may explain why the stress evolution is slower in the slip case than in the no-slip one as seen inFig. 10b (and also in Figs. 18a and 20, where there are slight changes between t = 30 and 60 s). Surprisingly,Fig. 10b shows that tr( τ ) avg is actually higher in the slip case, probably due to the absence of a bottomstationary zone. For benchmarking purposes we plot profiles of some dependent variables along the verticalcentreline in Fig. 20. 28 a) τ d = 40 Pa; G = 400 Pa (b) τ d = 40 Pa; G = 4 ,
000 Pa (c) τ d = 40 Pa; G = 40 ,
000 Pa (d) τ d = 70 Pa; G = 400 Pa (e) τ d = 70 Pa; G = 4 ,
000 Pa (f ) τ d = 70 Pa; G = 40 ,
000 Pa (g) τ d = 80 Pa; G = 400 Pa (h) τ d = 80 Pa; G = 4 ,
000 Pa (i) τ d = 80 Pa; G = 40 ,
000 Pa
Figure 17:
Continuous lines: contours of τ d predicted by the SHB model (top row: τ d = 40 Pa; middle row: τ d = τ y = 70 Pa; bottom row: τ d = 80 Pa) for various values of the elastic modulus G (left column: G = 400 Pa; middlecolumn: G = 4,000 Pa; right column: G = 40,000 Pa), the rest of the parameters having the values listed in Table 2,and the lid velocity being U = 0.1 m / s. Dashed lines: corresponding contours predicted by the HB model. It was mentioned in Sec. 2 that Cheddadi et al. [32] found that the steady state of a SHB flow (includingthe velocity field) can depend in a continuous manner on the initial conditions. That was observed insimple cylindrical Couette flow, where the shear stress is determined by the geometry and kinematics butthe circumferential normal stress depends also on its initial value, which affects the final extent of the yieldzone and therefore also the final velocity field. We performed analogous experiments to see if this can beobserved also in the present, more complex flow. So, we solved two more cases that are identical to the“base case” of Sec. 6.1, except for the initial stress conditions: τ = τ = + √ τ y (first case) or −√ τ y a) SHB flowfield (b)
HB flowfield (c) τ d = 75 Pa Figure 18: (a) SHB flowfield for the slip case; shaded region: unyielded material at t = 60 s; dashed lines: yieldlines at t = 30 s; continuous lines: streamlines (they are drawn at equal streamfunction intervals). (b) Correspondingsteady-state HB flowfield (the streamlines are drawn at the same values of streamfunction as in (a)). (c) The τ d = 75Pa contour of the SHB (continuous line) and HB (dashed line) flowfields. (a) No slip (b)
Slip
Figure 19:
Contours of ˙ γ [s − ] for U = 0.100 m / s without wall slip (a) (base case) or with wall slip (b). The thickcontinuous lines are yield lines and the dashed lines are streamlines. The latter are plotted at equal streamfunctionintervals (the same intervals in both figures). (second case) throughout, the remaining stress components being zero. These initial conditions correspond,respectively, to tension and compression states that are isotropic in the xy plane; from Eq. (9) it followsthat τ d ( t = 0) = τ y everywhere, i.e. the material is on the verge of yielding. The material is initially at rest( u = 0). The simulations were carried on until time t = 60 s.Interestingly, the velocity fields of both of these new cases arrive at practically the same steady state;this can be seen in Fig. 21, where the kinetic energies of both cases converge, and more clearly in Fig. 22awhere the respective streamlines can be seen to be identical. However, this steady-state is not the same asthat arrived at in the base case ( τ = τ = 0); Fig. 21 shows that unlike in the two new cases, the kineticenergy of the base case does not exhibit an overshoot and eventually increases to a value which is about2.25% larger than that of the other two cases. In Fig. 22a the base case streamlines are also not identicalto those of the other two cases. As far as the main vortex is concerned, the new cases have it located at(˜ x c , ˜ y c ) = (0 . , . ψ c = − . a) u and τ d at x = L/ (b) τ and τ at x = L/ Figure 20:
Profiles of some dependent variables along the vertical centreline, for the slip case. SHB results at t = 60and 30 s are plotted in continuous and dashes lines, respectively, and steady state HB results in dash-dot line. a large transition zone; this is shown in Fig. 22c where the transition zone is approximately the green region,throughout which τ d is slightly above τ y . Inside the transition region the τ d fields of the two new cases aresimilar (Fig. 22c), but the individual τ and τ stress components have opposite signs which they haveinherited from the initial conditions, as seen in Fig. 23; in fact, near the bottom wall these stress componentsretain values close to the initial ones, ±√ τ y ≈
121 Pa. On the other hand, the base case stresses are closeto zero in that region, again an inheritance of the initial conditions. This reflects on the much lower stresstrace values seen in Fig. 21. Outside of the transition region, where the rates of deformation are non-zero(including the plug zone), the two new cases have identical steady-state stress fields, as seen in Figs. 22cand 23, whereas those of the base case deviate slightly. This includes the yield line that forms the boundaryof the plug zone: Fig. 22b shows that the plug zones of the two new cases are identical, but that of the basecase is slightly smaller, which is likely the cause of the slightly different velocity field.
Our final investigation concerns the sudden stopping of the lid and the study of the subsequent flowdecay. Classic viscoplastic models are known to predict flow cessation in finite time after the driving agentis removed [84]. This can be proved rigorously using variational inequalities (e.g. upper bounds for thecessation times of some one-dimensional flows are derived in [85–87]), but roughly speaking it is due to the (a) t ∈ [0 ,
4] s (b) t ∈ [0 ,
30] s
Figure 21:
Histories of N.K.E. (Eq. (69)) and tr(˜ τ ) avg (Eq. (70)) for the cases with initial conditions τ = τ = −√ τ y (continuous lines), τ = τ = + √ τ y (dashed lines) and τ = τ = 0 (dash-dot lines, the base case). a) streamlines (b) yield lines (c) τ d Figure 22:
Steady-state streamlines (a) (streamfunction interval δφ/ ( ρU ) = 3 . × − ), yield lines (b), and τ d contours (c) for the cases with initial conditions τ = τ = −√ τ y (continuous lines), + √ τ y (dashed lines), and zero(dash-dot lines, only in (a) and (c)). In (c) the values of τ d , in Pa, are indicated next to each contour. The unyieldedregion of the τ = τ = −√ τ y case is shaded blue. The region where (cid:107) u (cid:107) < × − U is shaded green. (a) τ (b) τ (c) pressure Figure 23:
Normal stress components (a)-(b) and pressure (c) along the vertical centreline ( x = L/
2) for the caseswith initial conditions τ = τ = −√ τ y (continuous lines), + √ τ y (dashed lines), and zero (dash-dot lines). effect of the yield stress on the rate of energy dissipation: it keeps it high enough during the flow decay suchthat all of the fluid’s kinetic energy (KE) is dissipated in finite time – a rough explanation is sketched inAppendix C. The cessation of viscoplastic (Bingham) lid-driven cavity flow was studied in [40], accordingto the findings of which we would expect HB flow to cease completely at some time t c (cid:28) T c ≡ ρU L/S ≈ T c is the time needed for a force of magnitude SL to bring a mass of momentum ρU L to rest).The situation concerning SHB flow is expected to be somewhat different. Firstly, even though the energyconversion rate (cid:82) Ω τ : ∇ u dΩ is still expected to be large enough to convert all the KE of the material infinite time, this rate now includes not only energy dissipation, but also energy storage in the form of elastic(potential) energy, which can later be converted back into KE [88]. In fact, once all of the material becomesunyielded there is no mechanism for energy dissipation and it is expected that the remaining energy willperpetually change form from elastic to kinetic and vice versa, resulting in oscillatory motions. Furthermore,the formation of transition regions in previous simulations makes it uncertain whether all of the materialwill become unyielded in finite time. To investigate these issues, using the “steady-state” base flow (Sec.6.1) as initial condition, we suddenly stopped the lid motion and carried out a simulation to see how theflow evolves. We repeated the procedure also with the slip case of Sec. 6.4. Figure 24a shows the evolutionof the KE with time ( t = 0 is the instant the lid is halted). In neither case does the flow cease in finite time.In the no-slip case, the KE history is highly oscillatory, confirming the anticipated perpetual back-and-forth conversion between kinetic and elastic energies. The top diagram of Fig. 24b shows a clearer pictureover a narrower time window. The KE peaks 29 times within that 2 s window, i.e. a single oscillation lastsabout 69 ms; the time step is kept to about ∆ t = 2 × − s by the adjustable time step scheme, so each KE32 a) (b) Figure 24: (a): Fluid kinetic energy (Eq. (69)), normalised by its value at the instant the lid is halted, versus timeelapsed (the lid is halted at time t = 0) for the no-slip and slip cases. (b): The top diagram is a close-up view of aportion of the no-slip curve in (a). The bottom diagram plots the corresponding values of the force exerted by thefluid on the lid, normalised by τ y L . oscillation period is resolved into about 350 time steps. The diagram shows that the KE variation does notconsist of a single harmonic component, and Fig. 25 shows that the flow is indeed quite complex. Usually atany given instant there appears one main vortex rotating in either the clockwise or anticlockwise sense, butits location and shape are not fixed, while smaller vortices may also appear. The material is not completelyunyielded, and the yielded zones are transported along as the material oscillates while their size varies withtime. The existence of these yielded zones allows some energy dissipation, and hence the mean KE in Fig.24 very slowly drops. The bottom diagram of Fig. 24b shows the force F lid exerted by the fluid on the lid.This force is negative, pushing the lid towards the left, i.e. in the opposite direction than it was movingprior to its halting. The magnitude of the force oscillates slightly above the value τ y L , and an inspectionshows that this is because the magnitude of τ on the lid slowly drops towards τ y (the lid is in touch witha transition zone). Figure 25i shows that the magnitude of F lid drops when the lid touches a clockwisevortex (Figs. 25b, 25d, 25g) and it increases when it touches a counter-clockwise vortex (Figs. 25c, 25e). Acomparison between the KE and F lid diagrams in Fig. 24b shows that F lid oscillates at a frequency that isroughly half of that of the KE, which is expected since during a single F lid period both the clockwise andanticlockwise vortex velocities are maximised, which leads to two KE peaks.In the slip case, the KE peaks immediately after the lid halt (Fig. 24a) but then the oscillations dieout relatively quickly and the KE diminishes. The KE peak is due to the material recoiling right after thelid halts (Fig. 25j), which is possible as the walls provide limited resistance, due to slip. Then, the KEis dissipated through a mechanism that is absent in the no-slip case: as the material oscillates inside thecavity, it slides past the walls due to slip, and the friction at the wall / material interface dissipates the KE.Even in this case, however, transition zones exist long after the lid has halted (Fig. 25k).
7. Conclusions
This work presented a FVM for elastoviscoplastic flows described by the SHB model. The method isapplicable to collocated structured and unstructured meshes; it incorporates a new variant of momentuminterpolation, a variant of “two sides diffusion”, and the CUBISTA scheme, in order to suppress spuriouspressure, velocity, and stress oscillations, respectively. The method was shown to achieve this goal, and alsoto be consistent on both smooth and irregular meshes. Similar stabilisation techniques have recently beenapplied in a Finite Element context to allow use of equal-order polynomial basis functions for all variables[89]. An implicit temporal discretisation with adaptive time step is employed.33 a) t = 3 s (b) t = 3.05 s (c) t = 3.10 s (d) t = 3.15 s (e) t = 3.20 s (f ) t = 3.25 s (g) t = 3.30 s (h) colour map (i) History of F lid (j) (slip) t = 0.10 s (k) (slip) t = 60 s Figure 25: (a)-(g): Snapshots of the no-slip cessation flow field, with colour contours of τ d , of different base colourfor yielded (yellow) and unyielded (blue) regions (see the colour map (h)), and instantaneous streamlines plotted atstreamfunction intervals of δψ/LU = 0 . ψ = 0 as one of the contours. (i): Part of the history of the forceexerted by the fluid on the lid, normalised by τ y L , with the instants corresponding to snapshots (a)-(g) marked onthe plot. (j)-(k): Snapshots of the slip cessation flow field. The method was applied to the simulation of EVP flow in a lid-driven cavity, with the SHB parameterschosen so as to represent Carbopol. The results can serve as benchmark solutions, but furthermore thesimulations were designed so as to allow investigation of several aspects of the SHB model behaviour.Complementary simulations with the classic HB model were also performed for comparison. It was noticedthat the SHB model can predict “transition zones”, regions where the velocity is near-zero and it takesa very long, or infinite, time to transition from a formally yielded to an unyielded state. The differencesbetween the SHB and HB velocity fields are rather small for the lid velocities tested, although some elasticeffects are noticeable in the SHB case, such as a slight upstream displacement of the vortex, a downstreamstretching of the plug zone, and a small swelling of the bottom unyielded zone, compared to the HB results.These differences diminish by increasing the elastic modulus G , although the convergence of the SHB andHB yield lines is rather slow and requires very large values of G . The velocity fields converge much faster.We also applied slip at the walls, as EVP materials are slippery, in which case no transition zones developed.Motivated by the observation of Cheddadi et al. [32] that different initial conditions of stress can leadto different steady states even as concerns the velocity, we performed two additional simulations where theinitial state of the material was stationary but on the verge of yielding, with either compressive or tensile34esidual stresses, respectively. Indeed, although both of these simulations converged to the same steady-statevelocity field, that field was slightly but noticeably different than that obtained when the initial stresseswere zero. Due to this property of the SHB model, an accurate calculation of a steady state requires notonly sufficient spatial resolution (grid spacing) but also sufficient temporal resolution (time step size), whichis facilitated by the use of an adaptive time discretisation scheme such as the one proposed herein.Finally, simulations of the cessation of the flow after the lid is suddenly halted were performed, whichshowed that, unlike what is predicted by the HB model, the flow does not cease in finite time but there isa perpetual back-and-forth conversion between kinetic and elastic energies. In the no-slip case, a very slownet energy dissipation was observed, which is possible due to the persisting existence of yielded regions. Theenergy dissipation is much faster under wall slip, due to friction between the material and the walls.Incorporation of additional rheological phenomena such as thixotropy is planned for the future. Acknowledgements
This research was funded by the LIMMAT Foundation under the Project “MuSiComPS”.
Appendix A The SHB extra stress tensor in the limit of infinite elastic modulus
The SHB constitutive equation (7) distinguishes between τ and τ d , unlike the HB equation (4). For thetwo models to become equivalent as G → ∞ it is necessary that τ → τ d or, equivalently, tr( τ ) → G → ∞ the first term on the left-hand side of Eq. (7) diminishes. This makesyielded material (where the “max” term is non-zero) tend to behave as a generalised Newtonian fluid sothat tr( τ ) → ⇒ τ → τ d and the two models do become equivalent in yielded regions. On the other hand,in unyielded regions (where the “max” term is zero) Eq. (7) tends to reduce to ˙ γ = 0, which does not, atfirst glance, imply that τ → τ d . However, the following observation can be made. In an unyielded regionthe SHB model can be written as (substituting Eq. (8) into Eq. (7)) DτDt ≡ ∂τ∂t + u · ∇ τ = G ˙ γ + G (cid:0) ( ∇ u ) T · τ + τ · ∇ u (cid:1) (A.1)The material derivative Dτ /Dt is the rate of change of the stress tensor with respect to time in a fluidparticle moving with the flow. This derivative is applied component-wise to the stress tensor, and thereforethe rate of change of tr( τ ) is D (tr( τ )) Dt = tr (cid:18) DτDt (cid:19) (A.2)or, using Eq. (A.1), D (tr( τ )) Dt = G tr( ˙ γ ) + G tr (cid:0) ( ∇ u ) T · τ + τ · ∇ u (cid:1) (A.3)In the right-hand side, tr( ˙ γ ) → G → ∞ forces ˙ γ →
0. For the secondterm in the right-hand side, in index notation, we have( ∇ u ) T · τ + τ · ∇ u = (cid:88) i =1 3 (cid:88) j =1 (cid:32) (cid:88) k =1 ∂u i ∂x k τ kj + (cid:88) k =1 ∂u j ∂x k τ ik (cid:33) e i e j (A.4)where e i is the unit vector in the direction i . Taking the trace of the above expression we obtaintr (cid:0) ( ∇ u ) T · τ + τ · ∇ u (cid:1) = (cid:88) i =1 3 (cid:88) j =1 (cid:18) ∂u j ∂x i + ∂u i ∂x j (cid:19) τ ij (A.5)The terms in parentheses in the sum (A.5) are just the components of ˙ γ , and so they tend to zero as G → ∞ .Thus the whole right-hand side of Eq. (A.3) tends to zero in the unyielded regions as G → ∞ . Therefore,under these conditions, D (tr( τ )) /Dt →
0, i.e. tr( τ ) remains constant in any particle moving in an unyieldedregion. Given that tr( τ ) = 0 for yielded particles, as mentioned above, there appears to be no mechanismfor tr( τ ) to acquire any value different than zero. Thus we expect that τ → τ d also in the unyielded regions.But even if tr( τ ) (cid:54) = 0 ⇒ τ (cid:54) = τ d in such a region, for practical purposes the stress field would still beequivalent to a HB stress field where the isotropic part of τ is incorporated into the pressure.35 y c P c c N u u + u (cid:48) c h xy c P c c N u u + u (cid:48) c h (cid:48) = h/ Figure 26:
A local velocity perturbation u (cid:48) c with wavelength equal to the grid spacing h causes local perturbationsd u (cid:48) / d x to the velocity gradient that are proportional to u (cid:48) c /h = O ( h − ); it follows that the corresponding perturbationsto the second derivative of velocity d u (cid:48) / d x are proportional to O ( h − ) /h = O ( h − ). Appendix B Effect of grid refinement on the pressure stabilisation scheme
It follows from the definition of a mif (33) that as the grid is refined ( h f →
0) the viscous term in thedenominator dominates over the inertial one: at fine grids a mif ≈ h f / (2 κ + 2 η a ). This can be explained withthe aid of the simple one-dimensional example used for deriving the scheme. From Eq. (38), if we wantto locally perturb a velocity field at point c by u (cid:48) c (Fig. 26) by locally perturbing the pressure gradient by − d p (cid:48) / d x | c ≈ ( p (cid:48) P − p (cid:48) N ) /h then these perturbations are related by u (cid:48) c = 1 ρ (cid:18) d u d x (cid:12)(cid:12)(cid:12)(cid:12) c + 1∆ t (cid:19) + 2 η c h · p (cid:48) P − p (cid:48) N h = O ( h ) · d p (cid:48) d x (cid:12)(cid:12)(cid:12)(cid:12) c ⇒ d p (cid:48) d x (cid:12)(cid:12)(cid:12)(cid:12) c = O ( h − ) · u (cid:48) c (B.1)Thus, in order to effect a local velocity perturbation of wavelength equal to the grid spacing, as shownin Fig. 26, the pressure gradient must be adjusted locally by an amount that scales as ( h − ), i.e. it mustbecome larger and larger as the grid is refined. This is because it has to overcome the viscous resistancewhich scales as d u (cid:48) / d x | c = O ( h − ): grid refinement not only increases the velocity slopes (and associatedviscous stresses) but also changes these slopes over a shorter distance (Fig. 26). On the contrary, the partof the pressure gradient that is needed to balance the inertial contribution of u (cid:48) c to the momentum balanceis independent of the grid spacing. Conversely, Eq. (B.1) shows that fixed localised perturbations of thepressure gradient of wavelength equal to the grid spacing cause local velocity perturbations of magnitude O ( h ), which become smaller and smaller with grid refinement due to increased viscous resistance.The pseudo-velocities u p + f and u p − f defined by Eqs. (31) and (32), have magnitudes of O ( h ) (because a mif = O ( h ) and, in smooth pressure fields, both p P − p N and ∇ p · ( P − N ) are O ( h )), whereas the actualvelocity ¯ u c f in (30) does not diminish with refinement but tends to a finite value, the exact velocity atthe given point. Therefore, as the grid is refined the mass flux becomes dominated by the interpolatedvelocity ¯ u c f while the stabilisation pseudo-velocities u p + f and u p − f diminish. This may raise concerns thatthe stabilising effect of the scheme may also diminish. However, this is not the case. Suppose that spuriouspressure oscillations do arise, of amplitude ∆ p o (in the absence of preventive measures, this amplitude isunaffected by grid refinement [42]). The total pressure can be decomposed into a smooth part, which isclose to the exact solution, and a spurious oscillatory part: p = p s + p o . Then the sum of pseudo-velocitiesat point c of the grid shown in Fig. 4 is equal to u p + c − u p − c = a mic (cid:2) ( p sP + p oP ) − ( p sN + p oN ) − ¯ ∇ qh ( p s + p o ) c · ( P − N f ) (cid:3) = a mic (cid:104) ∆ p o + ( p sP − p sN ) − ¯ ∇ qh p sc · ( P − N f ) (cid:124) (cid:123)(cid:122) (cid:125) O ( h ) or O ( h ) (cid:105) ≈ a mic ∆ p o (B.2)because p oP − p oN = ∆ p o and ∇ qh p oP = ∇ qh p oN = 0 (the ∇ qh operator is insensitive to oscillations). Also,according to the discussion of Eq. (35), the underlined terms add up to O ( h ) (or O ( h ) on smooth structured36rids), which is small compared to ∆ p o . In other words, u p + f cancels out with the smooth part of u p − f butleaves out the oscillatory part. A similar consideration for point c P (Fig. 4) leads to u p + c P − u p − c P ≈ − a mic P ∆ p o (B.3)where the minus sign comes from the fact that, according to the nature of the spurious oscillation, if theoscillatory pressure component decreases by ∆ p o from P to N , then it increases by ∆ p o from P P to P .Then the combined contribution of faces c and c P to the image of the continuity operator for cell P (theleft-hand side of Eq. (34)) divided by the cell volume is1Ω P (cid:104) ˙ M c + ˙ M c P (cid:105) = 1 h ρ h (cid:2) (¯ u c + u p + c − u p − c ) − (¯ u c P + u p + c P − u p − c P ) (cid:3) which, using Eqs. (B.2) and (B.3), and also assuming that a mic ≈ a mic P ≈ Ch for some constant C ( C =2( κ + η a ) if the grid is fine enough, from Eq. (33)), becomes1Ω P (cid:104) ˙ M c + ˙ M c P (cid:105) = ρ ¯ u c − ¯ u c P h + 2 ρC ∆ p o (B.4)The first term on the right-hand side is the contribution of the actual velocities to the continuity image, andcan be seen to tend to the h -independent value ρ∂u/∂x | P with grid refinement. The second term on theright-hand side is the contribution of the pseudo-velocities, which can also be seen to be h -independent , i.e.it does not diminish with grid refinement, although the pseudo-velocities themselves do diminish comparedto the actual velocities, as mentioned above. If we take into account also the other two faces of cell P (thehorizontal ones) then the total contribution of the pseudo-velocities to the continuity image is 4 ρC ∆ p o . Ifwe repeat this analysis for the neighbouring cells N and P P , then it will turn out that the contributions ofthe pseudo-velocities to the continuity images of those cells are − ρC ∆ p o , i.e. they have opposite sign tothat of cell P . Overall, in the presence of pressure oscillations the contributions of the pseudo-velocities tothe continuity image over the entire grid have a checkerboard (oscillatory) pattern like the one shown in Fig.3, with amplitude proportional to the amplitude of the pressure oscillations, ∆ p o . As discussed in relationto Eq. (24), when solving the system of all continuity equations (34), because the right-hand side is smooth(it is zero), spurious pressure oscillations, which would produce an oscillatory right-hand side, are excluded.Thus, the effectiveness of momentum interpolation in suppressing spurious pressure oscillations does notdegrade with grid refinement. The amplitude of the oscillations produced by momentum interpolation tothe continuity image is proportional to the amplitude of the spurious pressure oscillations themselves, andindependent of the grid spacing h . Appendix C Energy dissipation and flow cessation
This appendix sketches an explanation for the finite cessation times of HB fluids, but for rigorous proofsthe reader is referred to the literature cited in Sec. 6.6. For inelastic fluids, once the lid motion ceases, therate of energy dissipation equals the rate of decrease of the kinetic energy of the fluid [88].dd t (cid:90) Ω 12 ρu · u dΩ = − (cid:90) Ω τ : ∇ u dΩ (C.1)For generalised Newtonian fluids ( τ = η ( ˙ γ ) ˙ γ ) we have τ : ∇ u = η ( ˙ γ ) ˙ γ : ∇ u = η ( ˙ γ ) ˙ γ (C.2)For a HB fluid, η ( ˙ γ ) = ( τ y + k ˙ γ n ) / ˙ γ and the energy dissipation rate (C.2) becomes ( τ y + k ˙ γ n ) ˙ γ ; in fact, thisexpression holds even in unyielded regions as it predicts zero energy dissipation there due to ˙ γ = 0. So, forHB flow the energy balance (C.1) becomesdd t (cid:90) Ω 12 ρu · u dΩ = − (cid:90) Ω ( τ y + k ˙ γ n ) ˙ γ dΩ (C.3) In HB unyielded regions the energy dissipation is zero due to ˙ γ = 0 since, by symmetry of the stress tensor, we have τ : ∇ u = τ : ∇ u + τ : ∇ u = τ : ∇ u + τ T : ∇ u T = τ : ( ∇ u + ∇ u T ) = τ : ˙ γ = 0. x at any instant t can be expressed as the productof a function of time, χ , and a function of space, u , as u ( x, t ) = χ ( t − t ) u ( x ) (C.4)where u ( x ) is the velocity at time t , if we set χ (0) = 1. We therefore assume that the velocity field retainsits shape in time, but just downscales by a factor of χ ( t − t ) at time t relative to time t . This assumptionis correct for the cessation of Newtonian flow if t is large enough [40], but obviously involves an error inthe viscoplastic case, where the unyielded regions expand in time; nevertheless, we will assume this errorto be acceptable, not invalidating the final conclusion. Then, let V ( t ) be a measure of the velocity in thedomain, e.g. the mean or maximum velocity magnitude. Whatever the precise choice of V , it follows from(C.4) that V ( t ) = χ ( t − t ) V (0). We can then normalise the velocity as u ( x, t ) V ( t ) = χ ( t − t ) u ( x ) χ ( t − t ) V (0) = u ( x ) V (0) ≡ ˜ u ( x ) ⇒ u ( x, t ) = ˜ u ( x ) V ( t ) (C.5)Substituting for u from (C.5) into the left and right side of Eq. (C.3) we get, respectively, (cid:90) Ω 12 ρu · u dΩ = V (cid:90) Ω 12 ρ ˜ u · ˜ u dΩ = C k V (C.6) (cid:90) Ω ( τ y + k ˙ γ n ) ˙ γ dΩ = τ y V (cid:90) Ω ˜˙ γ dΩ + k V n +1 (cid:90) Ω ˜˙ γ n +10 dΩ = C p τ y V + C v k V n +1 (C.7)where ˜˙ γ is the magnitude of ˜˙ γ ≡ ∇ ˜ u + ( ∇ ˜ u ) T (it is dimensional). Since ˜ u is not a function of t , neitherare the constants C k , C p and C v , and substituting (C.6) and (C.7) into (C.3) we are left with the followingordinary differential equation: d( V )d t = − C (cid:48) p τ y V − C (cid:48) v kV n +1 ≤ − C (cid:48) p τ y V (C.8)where C (cid:48) p = C p /C k and C (cid:48) v = C v /C k are positive constants. The above inequality means that the flow willdecay at least as fast as if the plastic term ( − C (cid:48) p τ y V ) was the only one present. In that case, solving theequation would give V ( t ) = V ( t ) − C (cid:48) p τ y t − t ) (C.9)so that the flow reaches complete cessation in finite time: V ( t c ) = 0 ⇒ t c = t + 2 V ( t ) / ( C (cid:48) p τ y ). Thus, theplastic dissipation (due to τ y ) alone suffices to bring the whole material to rest in finite time t c ; includingalso the viscous dissipation term, − C (cid:48) v kV n +1 , will only make this process faster.It is interesting to briefly look at the case τ y = 0, when the plastic term is absent from (C.8). If n = 1,the flow is Newtonian and the solution of the equation is V ( t ) = V ( t ) e − C (cid:48) vk ( t − t ) (C.10)This result is identical with what was found in [40], and means that the flow continuously decays but nevercompletely ceases. If n (cid:54) = 1 (power-law fluid) the solution is V ( t ) = (cid:20) V ( t ) − n − (1 − n ) C (cid:48) v k t − t ) (cid:21) − n (C.11)For n >
1, the exponent 1 / (1 − n ) is negative, the flow decays more slowly than in the Newtonian case asthe viscosity drops with flow decay, and complete cessation is never reached. However, in the n < / (1 − n ) is positive, andcessation is reached in finite time, V ( t c ) = 0 ⇒ t c = t + 2 V ( t ) − n / ((1 − n ) C (cid:48) v k ). Thus, the existence of ayield stress is not a prerequisite for finite cessation time (Fig. 27).38 igure 27: Decay of kinetic energy with time after lid motion cessation of two power law fluids (Eq. (1) with τ y =0, k = 1 Pa s n , and ρ = 1000 kg / m ) of exponents n = 1.0 and 0.5, respectively. The domain and grid are the sameas in Sec. 6, and the initial condition ( t = 0) is the steady state flow for a lid velocity of 1 m / s. The kinetic energy isnormalised by its value at t = 0. The lid is brought to a halt suddenly for n = 1, but gradually over a time intervalof 0.1 s for n = 0.5 (to avoid iterative solver convergence difficulties). Appendix D Second order accurate reconstruction of cell centre value from face values
Suppose that the exact (or approximated to at least second-order accuracy) values φ f of a quantity φ are known at the face centres of a cell P (Fig. 2), and we want to approximate from these values the valueof φ at the cell centre, φ P . One way to proceed is the following. The fact that the known values are locatedat the cell boundary provides an incentive to try to derive a scheme based on the divergence theorem. Let r ( x ) = x − P be the vector function that returns the relative position of location x relative to the centroid P . Then ∇ · r = ∇ · x = D where D = 2 or 3 in two- and three-dimensional space, respectively. We canthen use the product rule to the divergence of the product φr : ∇ · ( φr ) = ∇ φ · r + φ ∇ · r = ∇ φ · r + D φ
Integrating both sides over the whole cell P and applying the divergence theorem on the left-hand side weget F (cid:88) f =1 (cid:90) s f φ r · n f d s = (cid:90) Ω P ∇ φ · r dΩ + D (cid:90) Ω P φ dΩ (D.1)This is an exact equation, but now we will approximate each of the integrals by the midpoint rule: (cid:90) s f φ r · n f d s = φ ( c f ) ( r ( c f ) · n f ) s f + O ( h ) s f (D.2) (cid:90) Ω P ∇ φ · r dΩ = ∇ φ ( P ) · r ( P ) Ω P + O ( h ) Ω P = 0 + O ( h ) Ω P (D.3) (cid:90) Ω P φ dΩ = φ ( P ) Ω P + O ( h ) Ω P (D.4)where in (D.3) we have used that r ( P ) = P − P = 0. Substituting these into (D.1) we get φ ( P ) = 1 D Ω P F (cid:88) f =1 (cid:2) φ ( c f ) s f r ( c f ) · n f (cid:3)(cid:124) (cid:123)(cid:122) (cid:125) φ P + F (cid:88) f =1 O ( h ) s f Ω P (cid:124) (cid:123)(cid:122) (cid:125) = O ( h ) (D.5)The first term on the right hand side is the approximation φ P , as can be seen by comparing to Eq. (59).The second term on the right hand side is the difference between this approximation and the exact value φ ( P ), i.e. the error, which is O ( h ) because s f / Ω P = O ( h − ). Therefore, the approximation (59) appears tobe only first-order accurate, of the same order as simply setting φ P = φ ( c f ) for any arbitrary face f !39s it possible that this error estimation is too pessimistic and does not account for some error cancellationthat occurs when adding the truncation errors of the midpoint rule approximations (D.2) – (D.4)? In factthis happens to be the case and it turns out that the approximation (59) is second-order accurate, whichis easiest to prove by showing that it is exact for linear functions. Indeed, if it is exact for linear functionsthen expanding φ in a Taylor series about P we find that the error is φ ( x ) = φ ( P ) + ∇ φ ( P ) · r ( x ) + O ( h ) ⇒ φ P = [ φ ( P ) + ∇ φ ( P ) · r ( x )] P + O ( h ) P = φ ( P ) + ∇ φ ( P ) · r ( P ) + O ( h ) P = φ ( P ) + O ( h )where we have used the facts that our interpolation operator (59) is linear ([ φ + ψ ] P = φ P + ψ P ), that itis exact for linear functions (and φ ( P ) + ∇ φ ( P ) · r ( x ) is a linear function), that r ( P ) = 0, and finally that O ( h ) P = O ( h ) (which can be seen by replacing τ ij,f with O ( h ) in the formula (59) and noting that both s f ( c f − P ) · n f and Ω P are of the same order, O ( h ) in 2D and O ( h ) in 3D).It remains to show that the interpolation operator is exact for linear functions. Suppose a linear function φ ( x ) = φ ( P )+ g · r ( x ) where g is a constant vector (it is the gradient of φ ). We want to show that φ P = φ ( P ).Applying the approximation (59) to φ we get φ P = 1 D Ω P F (cid:88) f =1 (cid:2)(cid:0) φ ( P ) + g · r ( c f ) (cid:1) ( r ( c f ) · n f ) s f (cid:3) = φ ( P ) 1 D Ω P F (cid:88) f =1 (cid:2) ( r ( c f ) · n f ) s f (cid:3) + 1 D Ω P g · F (cid:88) f =1 (cid:2) r ( c f )( r ( c f ) · n f ) s f (cid:3) (D.6)To proceed further we need a couple of geometrical results. Figure 28 shows cell P divided into F trianglesby the dashed lines which connect the centroid P to the vertices v f . Consider the triangle ( P , v f , v f +1 ),which has face f as one of its sides. The product r ( c f ) · n f is the perpendicular distance from P to face f , and therefore the product ( r ( c f ) · n f ) s f is equal to the shaded area of Fig. 28, which is enclosed in arectangle with one side coinciding with face f , another parallel side through P , and two perpendicular sidespassing through the vertices v f and v f +1 , respectively. The area of the triangle ( P , v f , v f +1 ) is half that ofthis rectangle (in three dimensions, the volume of the cone-like shape obtained by joining P with the edgesof face f is one third of the volume of the shaded rectangular box). By adding the areas of all triangles weget the total area of cell P : 1 D F (cid:88) f =1 (cid:0) r ( c f ) · n f (cid:1) s f = Ω P (D.7)A second geometric result that will be needed is the following. The centroid of the triangle ( P , v f , v f +1 )is ( P + v f + v f +1 ) / P + 2 c f ) /
3. The centroid of the whole cell, P , is equal to the sum of the individualcentroids of the triangles, P f say, weighted by their areas, Ω f = (1 / r ( c f ) · n f ) s f : P = 1Ω P F (cid:88) f =1 P f Ω f ⇒ P Ω P = F (cid:88) f =1 (cid:0) P + 2 c f (cid:1) (cid:0) r ( c f ) · n f (cid:1) s f Substituting for Ω P from (D.7) (with D = 2), moving everything to the left hand side, and merging the twosums, we arrive at F (cid:88) f =1 r ( c f ) (cid:0) r ( c f ) · n f (cid:1) s f = 0 (D.8) A calculation involving Taylor expansions along the faces shows that, in the sum of the truncation errors of Eq. (D.2) overall cell faces, their leading order terms cancel out. Result (D.7) is also obtainable by setting φ = 1 in Eq. (D.1). f c f P v f v f +1 r ( c f ) Figure 28:
Notation concerning the geometry of cell P . where we have used also that c f − P = r ( c f ). In three dimensions, the exact same equation (D.8) holds, butis derived by noting that the volume of the cone- or pyramid-like shape whose base is face f and its apex isat P is Ω f = (1 / r ( c f ) · n f ) s f ( D = 3 in Eq. (D.7)) while its centroid is P f = ( P + 3 c f ) / φ P = φ ( P ): the interpolation scheme (59) is exact for linear functions, andis therefore second-order accurate. References [1] P. Coussot, “Bingham’s heritage,”
Rheologica Acta , vol. 56, pp. 163–176, 2017.[2] P. Coussot, “Slow flows of yield stress fluids: yielding liquids or flowing solids?,”
Rheologica Acta ,vol. 57, pp. 1–14, 2018.[3] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, “Yield stress materials in softcondensed matter,”
Reviews of Modern Physics , vol. 89, 2017.[4] E. C. Bingham,
Fluidity and Plasticity . McGraw-Hill, 1922.[5] W. H. Herschel and R. Bulkley, “Konsistenzmessungen von gummi-benzollsungen,”
Kolloid-Zeitschrift ,vol. 39, pp. 291–300, 1926.[6] K. v. Hohenemser and W. Prager, “ ¨Uber die ans¨atze der mechanik isotroper kontinua,”
ZAMM-Journalof Applied Mathematics and Mechanics/Zeitschrift f¨ur Angewandte Mathematik und Mechanik , vol. 12,pp. 216–226, 1932.[7] P. Saramito, “A damped newton algorithm for computing viscoplastic fluid flows,”
Journal of Non-Newtonian Fluid Mechanics , vol. 238, pp. 6–15, 2016.[8] E. Mitsoulis and J. Tsamopoulos, “Numerical simulations of complex yield-stress fluid flows,”
RheologicaActa , vol. 56, pp. 231–258, 2017.[9] P. Saramito and A. Wachs, “Progress in numerical simulation of yield stress fluid flows,”
RheologicaActa , vol. 56, pp. 211–230, 2017.[10] Y. Dimakopoulos, G. Makrigiorgos, G. Georgiou, and J. Tsamopoulos, “The PAL (Penalized Aug-mented Lagrangian) method for computing viscoplastic flows: A new fast converging scheme,”
Journalof Non-Newtonian Fluid Mechanics , vol. 256, pp. 23–41, 2018.[11] J. Piau, “Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges:Meso- and macroscopic properties, constitutive equations and scaling laws,”
Journal of Non-NewtonianFluid Mechanics , vol. 144, pp. 1–29, 2007. 4112] M. Dinkgreve, M. M. Denn, and D. Bonn, ““everything flows?”: elastic effects on startup flows ofyield-stress fluids,”
Rheologica Acta , vol. 56, pp. 189–194, 2017.[13] N. Dubash and I. Frigaard, “Propagation and stopping of air bubbles in Carbopol solutions,”
Journalof Non-Newtonian Fluid Mechanics , vol. 142, pp. 123–134, 2007.[14] W. F. Lopez, M. F. Naccache, and P. R. de Souza Mendes, “Rising bubbles in yield stress materials,”
Journal of Rheology , vol. 62, pp. 209–219, 2018.[15] J. Tsamopoulos, Y. Dimakopoulos, N. Chatzidai, G. Karapetsas, and M. Pavlidis, “Steady bubble riseand deformation in newtonian and viscoplastic fluids and conditions for bubble entrapment,”
Journalof Fluid Mechanics , vol. 601, 2008.[16] Y. Dimakopoulos, M. Pavlidis, and J. Tsamopoulos, “Steady bubble rise in Herschel–Bulkley fluids andcomparison of predictions via the Augmented Lagrangian Method with those via the Papanastasioumodel,”
Journal of Non-Newtonian Fluid Mechanics , vol. 200, pp. 34–51, 2013.[17] D. Fraggedakis, M. Pavlidis, Y. Dimakopoulos, and J. Tsamopoulos, “On the velocity discontinuityat a critical volume of a bubble rising in a viscoelastic fluid,”
Journal of Fluid Mechanics , vol. 789,pp. 310–346, 2016.[18] D. Fraggedakis, Y. Dimakopoulos, and J. Tsamopoulos, “Yielding the yield-stress analysis: a studyfocused on the effects of elasticity on the settling of a single spherical particle in simple yield-stressfluids,”
Soft Matter , vol. 12, pp. 5378–5401, 2016.[19] P. Saramito, “A new constitutive equation for elastoviscoplastic fluid flows,”
Journal of Non-NewtonianFluid Mechanics , vol. 145, pp. 1–14, 2007.[20] P. R. de Souza Mendes and R. L. Thompson, “A critical overview of elasto-viscoplastic thixotropicmodeling,”
Journal of Non-Newtonian Fluid Mechanics , vol. 187-188, pp. 8–15, 2012.[21] D. Fraggedakis, Y. Dimakopoulos, and J. Tsamopoulos, “Yielding the yield stress analysis: A thoroughcomparison of recently proposed elasto-visco-plastic (EVP) fluid models,”
Journal of Non-NewtonianFluid Mechanics , vol. 238, pp. 170–188, 2016.[22] I. Cheddadi and P. Saramito, “A new operator splitting algorithm for elastoviscoplastic flow problems,”
Journal of Non-Newtonian Fluid Mechanics , vol. 202, pp. 13–21, 2013.[23] S. L. Frey, M. F. Naccache, P. R. de Souza Mendes, R. L. Thompson, D. D. dos Santos, F. B. Link, andC. Fonseca, “Performance of an elasto-viscoplastic model in some benchmark problems,”
Mechanics ofTime-Dependent Materials , vol. 19, pp. 419–438, 2015.[24] A. Syrakos, G. Georgiou, and A. Alexandrou, “Solution of the square lid-driven cavity flow of a Binghamplastic using the finite volume method,”
Journal of Non-Newtonian Fluid Mechanics , vol. 195, pp. 19–31, 2013.[25] A. Syrakos, Y. Dimakopoulos, G. C. Georgiou, and J. Tsamopoulos, “Viscoplastic flow in an extrusiondamper,”
Journal of Non-Newtonian Fluid Mechanics , vol. 232, pp. 102–124, 2016.[26] P. J. Oliveira, F. T. Pinho, and G. A. Pinto, “Numerical simulation of non-linear elastic flows with ageneral collocated finite-volume method,”
Journal of Non-Newtonian Fluid Mechanics , vol. 79, pp. 1–43, 1998.[27] J. L. Favero, A. R. Secchi, N. S. M. Cardozo, and H. Jasak, “Viscoelastic flow analysis using the soft-ware OpenFOAM and differential constitutive equations,”
Journal of Non-Newtonian Fluid Mechanics ,vol. 165, pp. 1625–1636, 2010.[28] A. Afonso, M. Oliveira, P. Oliveira, M. Alves, and F. Pinho, “The finite-volume method in computa-tional rheology,” in
Finite-Volume Methods - Powerful Means of Engineering Design , ch. 7, pp. 141–170,In-Tech Open Publishers, 2012. 4229] A. Syrakos, Y. Dimakopoulos, and J. Tsamopoulos, “Theoretical study of the flow in a fluid dampercontaining high viscosity silicone oil: Effects of shear-thinning and viscoelasticity,”
Physics of Fluids ,vol. 30, p. 030708, 2018.[30] F. Belblidia, H. R. Tamaddon-Jahromi, M. F. Webster, and K. Walters, “Computations with viscoplas-tic and viscoelastoplastic fluids,”
Rheologica Acta , vol. 50, pp. 343–360, 2011.[31] P. Saramito, “A new elastoviscoplastic model based on the Herschel–Bulkley viscoplastic model,”
Jour-nal of Non-Newtonian Fluid Mechanics , vol. 158, pp. 154–161, 2009.[32] I. Cheddadi, P. Saramito, and F. Graner, “Steady Couette flows of elastoviscoplastic fluids arenonunique,”
Journal of Rheology , vol. 56, pp. 213–239, 2012.[33] L. Lacaze, A. Filella, and O. Thual, “Steady and unsteady shear flows of a viscoplastic fluid in acylindrical Couette cell,”
Journal of Non-Newtonian Fluid Mechanics , vol. 220, pp. 126–136, 2015.[34] A. Syrakos, G. C. Georgiou, and A. N. Alexandrou, “Performance of the finite volume method in solvingregularised Bingham flows: Inertia effects in the lid-driven cavity flow,”
Journal of Non-Newtonian FluidMechanics , vol. 208–209, pp. 88–107, 2014.[35] R. Sousa, R. Poole, A. Afonso, F. Pinho, P. Oliveira, A. Morozov, and M. Alves, “Lid-driven cavityflow of viscoelastic liquids,”
Journal of Non-Newtonian Fluid Mechanics , vol. 234, pp. 129–138, 2016.[36] R. da R. Martins, G. M. Furtado, D. D. dos Santos, S. Frey, M. F. Naccache, and P. R. de Souza Mendes,“Elastic and viscous effects on flow pattern of elasto-viscoplastic fluids in a cavity,”
Mechanics ResearchCommunications , vol. 53, pp. 36–42, 2013.[37] A. Syrakos, G. C. Georgiou, and A. N. Alexandrou, “Thixotropic flow past a cylinder,”
Journal ofNon-Newtonian Fluid Mechanics , vol. 220, pp. 44–56, 2015.[38] C. J. Dimitriou and G. H. McKinley, “A comprehensive constitutive law for waxy crude oil: a thixotropicyield stress fluid,”
Soft Matter , vol. 10, pp. 6619–6644, 2014.[39] I. Cheddadi, P. Saramito, C. Raufaste, P. Marmottant, and F. Graner, “Numerical modelling of foamCouette flows,”
The European Physical Journal E , vol. 27, 2008.[40] A. Syrakos, G. C. Georgiou, and A. N. Alexandrou, “Cessation of the lid-driven cavity flow of Newtonianand Bingham fluids,”
Rheologica Acta , vol. 55, pp. 51–66, 2016.[41] R. L. Thompson and E. J. Soares, “Viscoplastic dimensionless numbers,”
Journal of Non-NewtonianFluid Mechanics , vol. 238, pp. 57–64, 2016.[42] A. Syrakos and A. Goulas, “Estimate of the truncation error of finite volume discretization of theNavier-Stokes equations on colocated grids,”
International Journal for Numerical Methods in Fluids ,vol. 50, pp. 103–130, 2006.[43] J. Sun, M. Smith, R. Armstrong, and R. Brown, “Finite element method for viscoelastic flows basedon the discrete adaptive viscoelastic stress splitting and the discontinuous Galerkin method: DAVSS-G/DG,”
Journal of Non-Newtonian Fluid Mechanics , vol. 86, pp. 281–307, 1999.[44] A. Syrakos, S. Varchanis, Y. Dimakopoulos, A. Goulas, and J. Tsamopoulos, “A critical analysis ofsome popular methods for the discretisation of the gradient operator in finite volume methods,”
Physicsof Fluids , vol. 29, p. 127103, 2017.[45] C. M. Rhie and W. L. Chow, “Numerical study of the turbulent flow past an airfoil with trailing edgeseparation,”
AIAA Journal , vol. 21, pp. 1525–1532, 1983.[46] S. Zhang, X. Zhao, and S. Bayyuk, “Generalized formulations for the Rhie–Chow interpolation,”
Journalof Computational Physics , vol. 258, pp. 880–914, 2014.4347] G. B. Deng, J. Piquet, P. Queutey, and M. Visonneau, “Incompressible flow calculations with a con-sistent physical interpolation finite volume approach,”
Computers & fluids , vol. 23, pp. 1029–1047,1994.[48] C. Klaij, “On the stabilization of finite volume methods with co-located variables for incompressibleflow,”
Journal of Computational Physics , vol. 297, pp. 84–89, 2015.[49] H. M. Matos, M. A. Alves, and P. J. Oliveira, “New formulation for stress calculation: Application toviscoelastic flow in a T-junction,”
Numerical Heat Transfer, Part B: Fundamentals , vol. 56, pp. 351–371,2009.[50] M. Niethammer, H. Marschall, C. Kunkelmann, and D. Bothe, “A numerical stabilization frameworkfor viscoelastic fluid flow using the finite volume method on general unstructured meshes,”
InternationalJournal for Numerical Methods in Fluids , vol. 86, pp. 131–166.[51] F. Pimenta and M. Alves, “Stabilization of an open-source finite-volume solver for viscoelastic fluidflows,”
Journal of Non-Newtonian Fluid Mechanics , vol. 239, pp. 85–104, 2017.[52] C. Fernandes, M. S. B. Araujo, L. L. Ferras, and J. M. Nobrega, “Improved both sides diffusion(iBSD): A new and straightforward stabilization approach for viscoelastic fluid flows,”
Journal of Non-Newtonian Fluid Mechanics , vol. 249, pp. 63–78, 2017.[53] A. Syrakos, G. Efthimiou, J. G. Bartzis, and A. Goulas, “Numerical experiments on the efficiency oflocal grid refinement based on truncation error estimates,”
Journal of Computational Physics , vol. 231,pp. 6725–6753, 2012.[54] S. Patankar,
Numerical heat transfer and fluid flow . CRC press, 1980.[55] J. H. Ferziger and M. Peric,
Computational methods for fluid dynamics . Springer, 3rd ed., 2002.[56] F. Habla, A. Woitalka, S. Neuner, and O. Hinrichsen, “Development of a methodology for numericalsimulation of non-isothermal viscoelastic fluid flows with application to axisymmetric 4:1 contractionflows,”
Chemical Engineering Journal , vol. 207-208, pp. 772–784, 2012.[57] F. Moukalled, L. Mangani, and M. Darwish,
The Finite Volume Method in Computational Fluid Dy-namics . Springer, 2016.[58] M. A. Alves, P. J. Oliveira, and F. T. Pinho, “A convergent and universally bounded interpolationscheme for the treatment of advection,”
International Journal for Numerical Methods in Fluids , vol. 41,pp. 47–75, 2003.[59] X. Chen, H. Marschall, M. Schfer, and D. Bothe, “A comparison of stabilisation approaches for finite-volume simulation of viscoelastic fluid flow,”
International Journal of Computational Fluid Dynamics ,vol. 27, pp. 229–250, 2013.[60] F. Habla, M. W. Tan, J. Haßlberger, and O. Hinrichsen, “Numerical simulation of the viscoelastic flowin a three-dimensional lid-driven cavity using the log-conformation reformulation in OpenFOAM R (cid:13) ,” Journal of Non-Newtonian Fluid Mechanics , vol. 212, pp. 47–62, 2014.[61] R. Comminal, J. Spangenberg, and J. H. Hattel, “Robust simulations of viscoelastic flows at high weis-senberg numbers with the streamfunction/log-conformation formulation,”
Journal of Non-NewtonianFluid Mechanics , vol. 223, pp. 37–61, 2015.[62] P. H. Gaskell and A. K. C. Lau, “Curvature-compensated convective transport: SMART, a newboundedness- preserving transport algorithm,”
International Journal for Numerical Methods in Fluids ,vol. 8, pp. 617–641.[63] H. Jasak, H. Weller, and A. Gosman, “High resolution NVD differencing scheme for arbitrarily unstruc-tured meshes,”
International Journal for Numerical Methods in Fluids , vol. 31, pp. 431–449.4464] B. P. Leonard, “Bounded hihigh-order upwind multidimensional finite-volume convection-diffusion al-gorithms,” in
Advances in Numerical Heat Transfer (W. J. Minkowycz and E. M. Sparrow, eds.), vol. 1,pp. 1–57, Taylor and Francis, 1996.[65] W. Zhou, J. Ouyang, X. Wang, J. Su, and B. Yang, “Numerical simulation of viscoelastic fluid flows us-ing a robust FVM framework on triangular grid,”
Journal of Non-Newtonian Fluid Mechanics , vol. 236,pp. 18–34, 2016.[66] R. Fattal and R. Kupferman, “Constitutive laws for the matrix-logarithm of the conformation tensor,”
Journal of Non-Newtonian Fluid Mechanics , vol. 123, pp. 281–285, 2004.[67] A. Afonso, P. J. Oliveira, F. T. Pinho, and M. A. Alves, “The log-conformation tensor approach inthe finite-volume method framework,”
Journal of Non-Newtonian Fluid Mechanics , vol. 157, no. 1-2,pp. 55–65, 2009.[68] F. D. Vita, M. Rosti, D. Izbassarov, L. Duffo, O. Tammisola, S. Hormozi, and L. Brandt, “Elasto-viscoplastic flows in porous media,”
Journal of Non-Newtonian Fluid Mechanics , vol. 258, pp. 10–21,2018.[69] Y. Dimakopoulos and J. Tsamopoulos, “A quasi-elliptic transformation for moving boundary problemswith large anisotropic deformations,”
Journal of Computational Physics , vol. 192, pp. 494–522, 2003.[70] A. Syrakos and A. Goulas, “Finite volume adaptive solutions using SIMPLE as smoother,”
InternationalJournal for Numerical Methods in Fluids , vol. 52, pp. 1215–1245, 2006.[71] T.-W. Pan, J. Hao, and R. Glowinski, “On the simulation of a time-dependent cavity flow of anOldroyd-B fluid,”
International Journal for Numerical Methods in Fluids , vol. 60, pp. 791–808, 2009.[72] P. Saramito, “On a modified non-singular log-conformation formulation for Johnson–Segalman vis-coelastic fluids,”
Journal of Non-Newtonian Fluid Mechanics , vol. 211, pp. 16–30, 2014.[73] M. A. Hulsen, R. Fattal, and R. Kupferman, “Flow of viscoelastic fluids past a cylinder at high weis-senberg number: Stabilized simulations using matrix logarithms,”
Journal of Non-Newtonian FluidMechanics , vol. 127, pp. 27–39, 2005.[74] A. Poumaere, M. Moyers-Gonzlez, C. Castelain, and T. Burghelea, “Unsteady laminar flows of a Car-bopol gel in the presence of wall slip,”
Journal of Non-Newtonian Fluid Mechanics , vol. 205, pp. 28–40,2014.[75] T. C. Papanastasiou, “Flows of materials with yield,”
Journal of Rheology , vol. 31, pp. 385–404, 1987.[76] K. Sverdrup, N. Nikiforakis, and A. Almgren, “Highly parallelisable simulations of time-dependentviscoplastic fluid flow with structured adaptive mesh refinement,”
Physics of Fluids , vol. 30, p. 093102,2018.[77] M. R. Hashemi, M. T. Manzari, and R. Fatehi, “Non-linear stress response of non-gap-spanning mag-netic chains suspended in a newtonian fluid under oscillatory shear test: A direct numerical simulation,”
Physics of Fluids , vol. 29, p. 107106, 2017.[78] J. P´erez-Gonz´alez, J. J. L´opez-Dur´an, B. M. Mar´ın-Santib´a˜nez, and F. Rodr´ıguez-Gonz´alez, “Rheo-PIVof a yield-stress fluid in a capillary with slip at the wall,”
Rheologica Acta , vol. 51, pp. 937–946, 2012.[79] D. M. Kalyon, “Apparent slip and viscoplasticity of concentrated suspensions,”
Journal of Rheology ,vol. 49, pp. 621–640, 2005.[80] P. Panaseti, A.-L. Vayssade, G. C. Georgiou, and M. Cloitre, “Confined viscoplastic flows with hetero-geneous wall slip,”
Rheologica Acta , vol. 56, pp. 539–553, 2017.[81] P. Panaseti, K. D. Housiadas, and G. C. Georgiou, “Newtonian Poiseuille flows with pressure-dependentwall slip,”
Journal of Rheology , vol. 57, pp. 315–332, 2013.4582] S. G. Hatzikiriakos, “Wall slip of molten polymers,”
Progress in Polymer Science , vol. 37, pp. 624–643,2012.[83] M. Philippou, Y. Damianou, X. Miscouridou, and G. C. Georgiou, “Cessation of newtonian circularand plane couette flows with wall slip and non-zero slip yield stress,”
Meccanica , vol. 52, pp. 2081–2099,2017.[84] I. A. Frigaard,
Background Lectures on Ideal Visco-Plastic Fluid Flows , pp. 1–40. Springer InternationalPublishing, 2019.[85] R. Glowinski,
Numerical methods for nonlinear variational problems . Springer-Verlag, 1984.[86] R. Huilgol, B. Mena, and J. Piau, “Finite stopping time problems and rheometry of Bingham fluids,”
Journal of non-Newtonian Fluid Mechanics , vol. 102, pp. 97–107, 2002.[87] L. Muravleva, E. Muravleva, G. C. Georgiou, and E. Mitsoulis, “Unsteady circular Couette flow of aBingham plastic with the Augmented Lagrangian Method,”
Rheologica Acta , vol. 49, pp. 1197–1206,2010.[88] H. H. Winter, “Viscous dissipation term in energy equations,” in
Modular Instruction Series C: Cal-culation and Measurement Techniques for Momentum, Energy and Mass Transfer, Vol. 7 , pp. 27–34,American Institute of Chemical Engineers: New York, 1987.[89] S. Varchanis, A. Syrakos, Y. Dimakopoulos, and J. Tsamopoulos, “A new Finite Element formulation forviscoelastic flows: circumventing simultaneously the LBB condition and the high-Weissenberg numberproblem,”