A Wall Boundary Condition for the Simulation of a Turbulent Non-Newtonian Domestic Slurry in Pipes

The concentration (using a lesser amount of water) of domestic slurry promotes resource recovery (nutrients and biomass) while saving water. This article is aimed at developing numerical methods to support engineering processes such as the design and implementation of sewerage for concentrated domestic slurry. The current industrial standard for computational fluid dynamics-based analyses of turbulent flows is Reynolds-averaged Navier–Stokes (RANS) modelling. This is assisted by the wall function approach proposed by Launder and Spalding, which permits the use of under-refined grids near wall boundaries while simulating a wall-bounded flow. Most RANS models combined with wall functions have been successfully validated for turbulent flows of Newtonian fluids. However, our experiments suggest that concentrated domestic slurry shows a Herschel–Bulkley-type non-Newtonian behaviour. Attempts have been made to derive wall functions and turbulence closures for non-Newtonian fluids; however, the resulting laws or equations are either inconsistent across experiments or lack relevant experimental support. Pertinent to this study, laws or equations reported in literature are restricted to a class of non-Newtonian fluids called power law fluids, which, as compared to Herschel–Bulkley fluids, yield at any amount of applied stress. An equivalent law for Herschel–Bulkley fluids that require a minimum-yield stress to flow is yet to be reported in literature. This article presents a theoretically derived (with necessary approximations) law of the wall for Herschel–Bulkley fluids and implements it in a RANS solver using a specified shear approach. This results in a more accurate prediction of the wall shear stress experienced by a circular pipe with a turbulent Herschel–Bulkley fluid flowing through it. The numerical results are compared against data from our experiments and those reported in literature for a range of Reynolds numbers and rheological parameters that are relevant to the prediction of pressure losses in a sewerage transporting non-Newtonian domestic slurry. Nonetheless, the application of this boundary condition could be extended to areas such as chemical and food engineering, wherein turbulent non-Newtonian flows can be found.


Introduction
This research is aimed at experimentally and numerically supporting the design of a novel sanitation system for urban areas with emphasis on the reduction of water usage, the improved recovery of nutrients and biomass and the efficient transport of domestic slurry.The experiments provide information on the rheology of the slurry and the pressure losses it brings about (due to wall friction) while flowing (turbulent) through circular pipes and bends.Using our experiments, we wish to develop a suitable simulation methodology to enable industrial studies and applications such as the design of sewerage carrying concentrated domestic slurry.The scope of this article is thus restricted to exploring the potential of the current state-of-the-art in industrial computational fluid dynamics for modelling and not the theoretical study of the turbulent properties of domestic slurry.
We used a clay-based slurry for the experiments, which was confirmed to be rheologically similar to concentrated domestic slurry.Similarly to many industrial slurries (clay, coal, iron oxide, etc.), the experimental slurry also departs from Newtonian behaviour [1].A Newtonian fluid is a fluid for which the shear stress and the strain rate are related linearly through In Equation (1), µ is the molecular viscosity that depends on temperature and pressure but not the flow; τ is the stress tensor, whereas is the strain rate tensor, defined as Variables representing tensors are printed in bold, whereas vectors are underlined.Above, u is the velocity vector and ∇u is the velocity gradient tensor or ∂u j ∂x i .In the fields of rheology and non-Newtonian fluid mechanics, it is common practice to replace 2 as follows: We use γ from now on to conform with rheological practices.The stress due to viscosity in any fluid that deviates from Newtonian behaviour can be defined as where the function η(| γ|) is not directly related to γ but to its second invariant | γ|, the shear rate; γ and | γ| are related through γ : γ = tr( γT γ) (6) γ : γ is the Frobenius product of γ and tr is the trace of a square matrix, which is the sum of its diagonal elements.T is the transpose of a matrix.
In Equation ( 4), η is the apparent viscosity that, unlike the molecular viscosity, depends not only on the fluid's temperature and pressure but also on the flow conditions such as the shear rate and even the duration of the shear [1].

Herschel-Bulkley Fluids
For succinctness, this article does not provide a detailed description of non-Newtonian fluids and their properties, for which the reader is referred to Chabbra and Richardson [1].Nonetheless, a rheogram that plots the relation between the shear stress τ and the shear rate γ for various non-Newtonian fluids is shown in Figure 1.
One notices that a Herschel-Bulkley fluid displays shear-thinning (pseudoplastic behaviour), which is the reduction in the apparent viscosity with increasing shear rate.Further, Herschel-Bulkley fluids require a minimum shear stress before they flow as a fluid does.This minimum stress is called the yield stress, and hence such fluids are also called yield pseudoplastic fluids.Herschel and Bulkley [2] introduced a mathematical definition for Herschel-Bulkley fluids: τ = τ y + m γn (7) where τ y is the yield stress, m is the fluid-consistency index and n is the fluid-behaviour index (all are standard terms as used in Chabbra and Richardson [1]).Further, all terms in Equation ( 7) are scalars, with γ being the shear rate and τ being the second invariant of the stress tensor τ.Equation ( 7) is only valid when |τ| = τ ≥ τ y .If τ < τ y , γ = 0 (see Figure 1).In three dimensions with full tensor notation [3], Equation (7) reads Figure 1.Rheograms of time-independent flow behaviour (adapted from [1]).
If a Herschel-Bulkley fluid undergoes a laminar flow through a circular pipe, one can simply calculate the velocity profile and shear stress through the balance between the shear stress and the pressure gradient across the pipe.From Figure 2, one can write where p is the pressure; τ rz is the shear stress against the direction of the flow; R and D are the radius and the diameter, respectively; L is the length of the section under consideration; r is the radial distance from the pipe's centreline; and V z (r) is the axial velocity at a radial distance r from the centreline.Further, the average flow velocity V is defined as where the volumetric flow rate Q is Using Equations ( 7), ( 9) and (11), one obtains an implicit relation between V z (r), r and the wall shear stress τ W [4]: where φ = τ y τ W .
(The readers are referred to Chabbra and Richardson [1], Govier and Aziz [5] and Bird et al. [6,7] for details and relevant literature.Heywood and Cheng [8] provide a review of known expressions for the laminar flow of Herschel-Bulkley fluids through circular horizontal pipes.)Akin to Equation (12), expressions that relate the wall stress to the flow velocity for the turbulent flow of Herschel-Bulkley fluids have also been proposed.Heywood and Cheng [8] mention that such expressions are correlations of experimental data and are not fully theoretical.Two such semi-theoretical expressions were proposed by Torrance [9] and Hanks [10].Further, Wilson and Thomas [11] proposed a correlation based on the relation between the dissipative small-scale turbulent eddies and the thickness of the viscous boundary layer.This was extended to Herschel-Bulkley fluids [12] and compared against experiments by Slatter [13].Slatter [13] reported that the correlation by Thomas and Wilson [12] is accurate in the early turbulent regime but diverges when the shear stress (or the flow velocity) increases, leading to more turbulence.It is also noted that semi-theoretical models are sensitive to the rheological properties of the fluids.

Computational Fluid Dynamics (CFD) and Non-Newtonian Fluids
With the development of computational resources, analyses of many industrial flows are being increasingly made as numerical computations of the Navier-Stokes (NS) equations.However, turbulent flows display a range of length and time scales; resolving all of these is not feasible at high Reynolds numbers [14].Therefore, these scales are mathematically altered (a process called decomposition) to reduce the range of length and time scales to a manageable quantity while ensuring that the decomposed scales represent the macroscopic properties of the flow under consideration.
One such approach is called Reynolds-averaged Navier-Stokes (RANS) modelling, in which the turbulent scales are ensemble-averaged (an average over many instances of the flow) to obtain a time-independent representation of the flow [15].It is assumed (with most RANS models) that the flow under consideration is fully turbulent.The following equations in indicial notation result from the Reynolds-decomposition and time-averaging of the NS equations for an incompressible Newtonian fluid in the absence of body forces.
In the above equations, the velocity u i has been decomposed as u i = u i + u i : the sum of an ensemble-average, u i , and a quantity that fluctuates across instances, u i .One now solves the above equations to obtain an ensemble-averaged velocity and pressure field, u and p.If the molecular viscosity is not constant but is a function of the shear rate, the Reynolds decomposition of the diffusive term will result in an ensemble-averaged and a fluctuating value of the effective viscosity [16].Upon including this term, Equation ( 14) reads where η is the ensemble-averaged effective viscosity of the non-Newtonian fluid under study.
To solve Equation ( 14), one needs a closure for −ρu i u j , known as the Reynolds stress.
RANS methods model −ρu i u j in terms of the known velocity field (the ensemble-averaged u in this case) to close Equations ( 13) and (14).Wilcox [17] presents a detailed description of RANS models, of which the κ- [18] and the Reynolds stress model (RSM; see Launder et al. [19] and references therein) are used in this article.For brevity, the readers are requested to refer to Wilcox [17] and the references therein for details.
On the other hand, to solve Equation (15), one also needs a closure for the term η γ ij mentioned as the fluctuating-viscosity stress tensor in literature [20,21].Known closures include those proposed by Pinho et al. [21] and Leighton et al. [22]; the latter was based on a priori direct numerical simulation.However, the validation of such closures is restricted to very low Reynolds numbers, as direct numerical simulation is computationally intensive [23].In fact, higher-Reynolds-number-based experimental (or numerical) studies that could be used to ascertain the contribution of the fluctuating viscosity tensor to RANS dynamics and validate appropriate closures are hitherto unavailable, thus limiting the application of such closures to theoretical studies.Following an appropriate definition of the Reynolds number used by Rudman and Blackburn [23] for their numerical studies, only the lowest Reynolds numbers observed in our experiments are perhaps comparable with studies used to validate the proposed closure.
Further, the closure proposed by Gavrilov and Rudyak [16] assumes that the fluctuations in effective viscosity (η ) are smaller than the ensemble-averaged viscosity.In a turbulent flow, the turbulent viscosity brings about a higher amount of dissipation as compared to the molecular (or effective) viscosity [24].Under such circumstances, the fluctuations themselves, if considered small, may not affect the dissipation of energy to a great extent and could possibly be ignored.
On the contrary, there is relatively more literature on the usage of RANS models (without a closure for the fluctuating viscosity) for the simulation of various types of non-Newtonian fluids, of which Malin [25][26][27] is amongst the first; κ-and κ-ω [24] were used.Both models were noted to deliver similar results despite different formulations and simulation methodologies.It is fair to mention that the model was not specifically calibrated for non-Newtonian fluids but was used with model constants determined for Newtonian fluids; however, the turbulent viscosity was numerically dampened near the walls with a suitable function.For Herschel-Bulkley fluids, at high Reynolds numbers (∼10 5 ), the simulated predictions converged onto values determined by a semi-empirical equation proposed by Dodge and Metzner [28] (see Section 1.3).
Similarly, Bartosik [29] numerically modified the eddy viscosity to account for non-Newtonian behaviour (viscous effects) near the walls, where it is more dominant than turbulence.Bartosik [30] used the modifications to simulate experiments conducted by Slatter [13] on Herschel-Bulkley fluids.It was noted that the concerned modification improves the accuracy of the κ-model even when the yield stresses are high (τ y ∼ 5 Pa or more).The above-mentioned numerical studies have experimental counterparts that support a mathematical modification of the different behaviour of non-Newtonian fluids near walls [28,31,32] (see Skelland [4] for a detailed overview).However, the modifications these studies suggested are applicable to power-law fluids (see Figure 1).
Other computational fluid dynamics (CFD)-based analyses include those of Singh et al. [33] (and references therein) and Gnambode et al. [34], which are based on more computationally expensive techniques, namely, direct numerical simulation and large eddy simulation, respectively.Apart from these, there is literature relevant to fields such as biofluid mechanics, anaerobic digestors, and so forth, concerning the use of standard RANS models for non-Newtonian fluid flows.

Wall-Bounded Flows
This work concerns the wall-bounded turbulent flow of Herschel-Bulkley fluids.For any fluid, as the Reynolds number of the flow increases, the effects of viscosity diminish [35].Nonetheless, in the presence of walls, all flows slow down to zero velocity at the wall (the no-slip condition), and in doing so, even a turbulent flow passes through a laminar region near the wall [36].Therefore, RANS models such as the κ-that are meant for the analysis of fully turbulent flows require a correction to model the transition into a laminar regime and ultimately zero velocity.
Such empirical corrections have permitted the use of RANS models for a range of industrial applications.This article examines one possible correction that may enable the application of RANS models for the simulation of Herschel-Bulkley fluids.The effects of the fluctuating viscosity tensor have been ignored in light of the fact that the Reynolds number in this study was relatively high and little experimental support is available to verify this effect.Another alternative would also be a modification for the model constants of various RANS models.However, this would also require detailed experimental data such as velocity fields and their fluctuations in time, for fully turbulent flows of non-Newtonian fluids (at Reynolds numbers higher than those reported in existing literature).
We therefore work on a simple modification of the behaviour near walls, where turbulence diminishes (as done by Malin [27], Bartosik [30] for high Reynolds numbers) as an alternative numerical model.We would like to mention again that our purpose is to explore suitable modelling techniques and not to formulate a theoretical basis for the turbulence of non-Newtonian fluids.
To begin with, Equation ( 16) defines the wall function developed by Launder and Spalding [18] on the basis of the velocity profile of a turbulent Newtonian fluid in a smooth tube proposed by Prandtl [37,38], which was later extended to the universal law of the wall [39].For a Newtonian fluid, where K ∼ 0.41 is the von Kármán constant and C ∼ 5.5 for a smooth wall.Therefore, Equation ( 16) can also be expressed as by merging the constant in Equation ( 16) with ln(y + ).E is another constant that equals 9.973.Equation ( 16) is not valid in the immediate vicinity of the wall but in the regions close to the wall where the flow is still turbulent.Immediately next to wall (see Figure 3), The terms y + and u + are used to non-dimensionalise the distance from a wall (y) and the velocity at that distance parallel to the wall (u).
In Equations ( 20)-( 22), u τ is called the friction velocity.With the κ-model, the first computational node is placed at a distance of 30 ≤ y + ≤ 300 from a wall, and the above equations (procedure defined in Launder and Spalding [18]) are used to estimate the wall shear stress from the local velocity field [40].

Non-Newtonian Wall Functions
Equivalents of Equation ( 16) for non-Newtonian fluids have also been proposed.Dodge and Metzner [28] used dimensional analysis and experimental data to propose an equation similar to Equation ( 16) for a power-law fluid (τ = m γn ) but with different coefficients.
where A n = 5.66 n 0.75 ln10 (24) and Clapp [31] (see Skelland [4] for details and more references) used Prandtl's concept of mixing length [37,38] to derive an expression for the turbulent region of a wall-bounded flow of a power-law fluid.Clapp's equation has a form similar to that of Equation ( 23) but with the following experimentally determined constants: A n = 2.78 n (26) The experimental data reported by Clapp [31] deviated by no more than 5% from Equations ( 23), ( 26) and ( 27) for 0.698 ≤ n ≤ 0.813.It is important to note that the y + in the above equations is not the same as in Equation ( 20), but is Other expressions for the velocity profile of wall-bounded non-Newtonian flows have also been proposed.Skelland [4] provides a detailed summary of such expressions.Additionally, Sawko [32] compares the above-mentioned wall models and a revised wall modelling approach (developed by Sawko [32]) based on analytical and experimental analyses of turbulent boundary layers.

Approach
The preceding section mentions the few strategies that exist for the modelling of wall-bounded non-Newtonian flows, although only those applicable to power-law fluids that yield to a non-zero stress.In other words, expressions such as Equations ( 23)- (25) have not yet been proposed or tested for yield viscoplastic and yield pseudoplastic fluids, of which the latter concerns us here.
This section elucidates the derivation and subsequent validation of an expression that could be used for modelling domestic slurries with Herschel-Bulkley rheology in a turbulent flow regime within a circular horizontal pipe.The approach is similar to that followed by Clapp [31], which is based on Prandtl's concept of mixing length [37,38] but includes the influence of yield stress as well.It is fair to mention that the proposed approach is based on certain approximations and that the intention is not to propose an exact and theoretically correct expression for the turbulent behaviour of Herschel-Bulkley fluids.
This article is divided into the following sections.Section 2 summarises the derivation and describes the NS solver and the related numerical challenges.Section 3 presents an overview of the experimental test cases considered in this article.Section 4 contains the results obtained from numerical analyses on the test cases mentioned in Section 3, and Section 5 reflects the conclusions and recommendations.

Methodology
This section provides details on the solver, the RANS models and numerical methods used for this article.This is followed by a derivation for the wall modification for Herschel-Bulkley fluids.

Solver and Numerics
ANSYS FLUENT [40], which is based on a finite-volume method, has been used to solve the RANS Equations ( 13) and ( 14) mentioned in Section 1.2.The spatial discretisation was done with a second-order upwind scheme to ensure numerical stability.The pressure was resolved with second-order accuracy but was decoupled from the velocity field with the semi-implicit method for pressure-linked equations (SIMPLE).
Amongst the RANS models, κ-was used with the standard model constants, whereas the RSM was used in two different configurations.(Please refer to ANSYS [40] for details.We did not use the κ-ω model because of the associated computational costs.)

•
The Reynolds stresses in the wall-adjacent cells are calculated explicitly in terms of the wall shear stress.This setting shall be referred to as RSM1.

•
A transport equation for κ is solved to obtain the Reynolds stresses in the cells adjacent to the wall, which shall be quoted as RSM2.
To implement the constitutive relation Equation ( 7), FLUENT uses a biviscosity model similar to that proposed by Tanner and Milthorpe [41] (see Mitsoulis [42] for a summary of similar approaches) [40].This approach prevents numerical instability arising from the yield stress at zero strain.

Wall Modelling: Specified Shear Approach
The following equations describe the wall model implemented in ANSYS FLUENT (from ANSYS [40]), which is based on Launder and Spalding [18]: where the subscript 1 refers to the first grid point next to a wall; κ is the turbulence kinetic energy (TKE) and E is an empirical constant equal to 9.793, as mentioned in Section 1.3; u * and y * are similar to u + and y + (see Equations ( 20)-( 22)) but with the wall shear expressed in terms of turbulence parameters.
As per ANSYS [40], the quantities y + and y * are nearly equal in equilibrium boundary layers (verified through post-processing).
Once a flow field is initialised, Equation ( 31) is used to obtain y * , after which u * is calculated from Equation ( 29) and then used to obtain τ W from Equation (30).This τ W is used in the RANS equations to obtain a new velocity field, following which the process is repeated until a converged solution is obtained.We consider a solution as converged once the iterative residuals for continuity, velocity, κ and are below 10 −6 .

An Appropriate Reynolds Number
To categorise the flow analysed in this article, one requires a suitable Reynolds number.Because the viscosity of a non-Newtonian fluid changes with the local flow conditions, one cannot use the standard Reynolds number.Instead, Rudman et al. [43] have proposed a Reynolds number based on the wall effective viscosity defined below.This viscosity can only be determined once the wall shear stress is known, and hence it serves only as a post-experimental (numerical) parameter.
We name this apparent viscosity η W . Then using Equation ( 7), replacing τ with the wall stress τ W and expressing γ as τ W η W (see Figure 2 for details), one obtains Rudman et al. [43] further used this η W to define a Reynolds number as where V is the average flow velocity and D is the pipe diameter.Rudman et al. [43] mentioned that this definition of the Reynolds number was more pertinent to the analysis of wall-bounded flows as it reflected the behaviour of Herschel-Bulkley fluids in the near-wall region, as compared to the more commonly used Metzner-Reed Reynolds numbers [44].Re R is used while presenting the results of the simulations.

Specified Wall Shear Based on Prandtl's Mixing Length
Following the approach put forth by Launder and Spalding [18], we derive a relation between the wall shear stress and the velocity near a wall for a Herschel-Bulkley fluid.The derivation is based on Prandtl's theory [37,38].The comprehensive book by Skelland [4] presents a clear derivation for the original law of the wall for a Newtonian fluid.
Prandtl began by describing the overall shear stress (see Figure 2) at any distance y from the wall as the sum of laminar and turbulent stresses: Prandtl's concept of mixing length (proposed originally for a Newtonian fluid) expresses τ rzturb as The total stress τ rz can be written in terms of τ W using Equation ( 9): Replacing τ rzlam with Equation ( 7), Near the wall, one may make certain assumptions regarding the laminar and turbulent stresses, as is usually done for Newtonian fluids and is done for non-Newtonian fluids by Dodge and Metzner [28] and Clapp [31].Firstly, y << R. Secondly, immediately near the wall within the laminar sublayer and up to the transition zone, the laminar stress is more than the turbulent stress (as also done by Prandtl).Therefore, From Equation (39), Integrating with respect to y and using the boundary condition u = 0 at y = 0, one obtains If τ y = 0, Equation ( 41) reduces to that proposed by Clapp [31] for a power-law fluid.A further n = 1 reduces the equation to that for a Newtonian fluid (law of the wall, Prandtl).
In contrast to the laminar region, within the turbulent core but still near the wall, the turbulence-related stresses outweigh the laminar stresses.
Although the laminar stress is low, τ y in Equation ( 42) is retained to include τ W − τ y in this analysis, which is also present in Equation (33).
Upon integration with respect to y, where C is the constant of integration.Further mathematical modifications to the logarithmic term lead to u In the previous step, the term comprising τ W and τ y is absorbed inside C. As illustrated in Skelland [4] for a Newtonian fluid, one may express the constant C in terms of other known terms (and the unknown τ W ).
Rewriting Equation (45) using Equation ( 46), If τ y = 0 and n = 1, Equation ( 47) reduces to the law of the wall or Equation ( 18) with u + and y + having their Newtonian definitions.This also implies that one must set C in Equation ( 47) as E .Doing so, one finally obtains Equation ( 48) is referred to as ψ 1 .
While deriving ψ 1 , it was assumed in the step corresponding to Equation ( 42) that retaining τ y does not alter the assumption that the laminar stresses are much lower than the turbulent stresses in the core turbulent region.However, this also implies that τ y must itself be low as compared to τ W .
Before discussing the effect of neglecting the magnitude of the yield stress, one must understand the flow regime of a Herschel-Bulkley fluid inside a smooth pipe.Figure 4 depicts the velocity profile of a Herschel-Bulkley fluid in a laminar flow.
At the centre of the pipe, the shear stress must reduce to zero, which in the case of a Herschel-Bulkley fluid implies a region that does not yield but flows as a plug, as shown in Figure 4.In the case of a turbulent flow, this plug is not as smooth as is shown in the Figure 4 for a laminar flow.However, direct numerical studies by Rudman and Blackburn [23] (at Re R = 7000) show the existence of a region characterised by a high effective viscosity and a smoother, nearly uniform velocity profile.Further, this region also displays reduction in turbulence and hence turbulent mixing, which possibly could lead to laminar pockets coexisting with unsteady turbulent flow until the wall, followed by a transition to a laminar flow consistent with a wall boundary.As a result, the effective mixing region that could sustain turbulence is reduced to that shown in Figure 4.The extent of this region is controlled by the value of the yield stress defined by Equation (37).Assuming this region of low (or zero) turbulence is circular is shape (not the case, as shown by direct numerical simulation, but simplifies the explanation provided below), the region will have a stress equal to the yield stress.Therefore, While deriving Equation ( 44), we neglected the fact that the turbulent stress term involves a mixing length or Ky that is directly proportional to the distance from the wall [37,38].This mixing length is an estimate of the distance over which a parcel of fluid retains its properties before mixing with the surrounding fluid.For a Newtonian flow, this implies that the maximum mixing length is R.
However, for a Herschel-Bulkley fluid, at a distance y = R − r from the pipe's wall, a parcel of fluid must merge into the unyielding plug; one may hypothesise that the mixing length could be thought as being reduced by a factor that depends on the size (radius) of the plug.One may rewrite y as y = R(1 − r R ) or y = R(1 − τ y τ w ).For clarity, ζ is defined as follows: Upon introducing this into Equation ( 42), thus redefining the mixing length for a wall-bounded Herschel-Bulkley fluid as Ky(1 − ζ).
Following the procedure used to derive ψ 1 with Equation (51) as the basis for the turbulent regime, one obtains another equation, ψ 2 , that reads Further modifications given by [28] to define a suitable y + for a power-law fluid could also be extended to Herschel-Bulkley fluids as Equation ( 53) and its relation with u + would require experimental verification, which we are unable to perform.Therefore, we shall avoid commenting on the velocity profile of the slurries that this article concerns.
One notices that both ψ 1 and ψ 2 are implicit in terms of τ W for a given value of u in the first cell near a wall boundary.Thus unlike the original wall function described in Section 1.3, ψ 1 and ψ 2 must be implemented as specified shear boundary conditions.
An example of such an approach is the use of the Monin-Obukhov similarity theory [45] to relate the shear stress on the ground to the velocity field at the first computational node, while simulating the atmospheric boundary layer using large-eddy simulation [46,47].As the simulation progresses in time, the implementation of a theoretical (or semi-empirical) relation between the instantaneous wall shear stress and the instantaneous velocity at the first grid point results in an accurate velocity profile (not instantaneous but time-averaged) of the atmospheric flow near the earth's surface.
This study aims to do the same but with RANS.As the simulation progresses towards a time-independent ensemble-averaged solution, the implementation of ψ 1 or ψ 2 will impose a theoretical constraint on the velocity field (the iterated values) and the wall shear.Thus, upon beginning the simulation with a uniform velocity field equalling the average flow velocity in the pipe, one expects to converge onto a solution that satisfies ψ 1 or ψ 2 near the wall for a Herschel-Bulkley fluid, which is equivalent to using the original wall function for a Newtonian fluid.
Although the original wall function provided an input to generate an accurate velocity profile of the flow with the distance from the wall, our aim is to only accurately predict the wall shear stress that could exist for a given average flow velocity in the pipe.This is purely because we are unable to verify the velocity (or turbulence) profile as a result of the lack of relevant experimental data.

Mesh
The computational grid for simulating fluids in horizontal pipes is shown in Figure 5. Instead of simulating the entire circular cross-section, two symmetry boundaries are used to simulate a quarter-pipe.The boundaries mathematically reduce the gradient across them to zero and allow no flux to leave or enter the quarter-pipe.The inlet has a specified velocity and turbulence intensity, whereas the outlet has a simple outflow condition with a diffusion gradient and a mass flow balance.The pipe walls are modelled with the specified spear approach mentioned in Section 2.2, that is, As one can see in Figure 5, the cells along the pipe increase in height from the pipe's wall towards the centre.This pattern is a part of the standard treatment used with wall functions while simulating turbulent boundary layers [17,40].In every case, the grid follows the same pattern from the wall until half the radius, following which the cells are more uniformly arranged into an O-grid.
The height of the first cell, 2y 1 (y 1 was previously used as the first grid point, which is located at half the height of the first cell) must be set as per guidelines similar to those for Newtonian fluids (see Section 1.3).
The η W and Re R values were calculated for each test case, following which, the distance for the first grid point or y 1 was obtained by setting Equation ( 53) to values from 20 to 120.Once the first grid cell's height was known, cells were added with an expansion ratio of between 1.05 and 1.1 until the cells reached the central O-grid, following which the O-grid was created.In either case, the cells along the pipe were sized to maintain an aspect ratio below 40.There is no specific literature on aspect ratios.However, on the basis of multiple guidelines available freely, the maximum aspect ratio was restricted to 40.
To ascertain the grid convergence, the height of the first grid point was varied between y + = 20 and y + = 200 (higher values of y + were possible only with higher Reynolds numbers, Re R ).It is noted that for Re R < 2 × 10 4 , y + for a converged value of the wall shear stress was closer to 20, whereas for Re R > 10 5 , y + was at least 60 (as high as 200).Thus, the first grid point was placed progressively further from the wall (in terms of y + ) as the Re R value increased.
For the standard Newtonian wall function, one should be aware that the first cell must be placed between 30 < y + < 300 (Newtonian definition of y + ) to avoid being completely within the fully laminar or turbulent regimes [14,35].This has been corroborated through various studies.As a result of a lack of experimental or numerical evidence, it is impossible to conclude the same for a non-Newtonian fluid, and hence we refrain from discussing any trends in the relation between y + , u + and τ W .Therefore, the data reported here is that obtained with values of y + (Herschel-Bulkley definition), beyond which the τ W did not change (as mentioned previously, this y + increased with Re R ).

Experiments
Table 1 presents an overview of the various experimental test cases considered in this article.All cases deal with Herschel-Bulkley domestic slurries in circular horizontal pipes.As per Equation ( 7), the unit of m for a given behaviour index n is Pas n .Table 1.The test cases considered for this article.The first two are from Slatter [13], PARK1 is from Park et al. [48], and the last two belong to experimental data supplied by the co-authors [49].In all simulations, the length of the pipe was set to L = 30 to 45 m (longer for a thicker pipe, L D > 300) to ensure that the flow reached a fully developed state and to avoid interference from the outflow boundary.The centreline velocity was monitored as it developed across the length of the pipe to a value that stopped varying with the downstream distance (after a certain distance from the inlet, thus reaching a fully developed state) before collecting any data.Slatter [13] mentions a maximum experimental error of about 7% in measuring the wall shear stress when it is about 1 Pa.As the flow velocity and hence the wall shear increase, the maximum error reduces to under 1% for a wall shear of 10 Pa.The error in the rheological measurements is a maximum of 1.46% [50].Park et al. [48] mention that their rheological parameters were obtained by linear regression of viscometric data, with a coefficient of determination of 0.996.

Case
The rheological parameters of the clay slurries S17 and S21 were determined using a Couette-flow co-axial rheometer of the Searle type, with a rotating inner cylinder and a stationary outer cylinder (Anton Paar MCR302; cup diameter of 29.29 mm, bob diameter of 27 mm and bob length of 40.5 mm; see www.anton-paar.comfordetails).(Clay slurries that were determined to possess a rheology similar to domestic slurry were used for the actual experiments.)The temperature was controlled to an accuracy of 0.1 • C using a Peltier system.After a rheogram was obtained using the viscometer, the rheological parameters were estimated for a Herschel-Bulkley behaviour with an absolute normalised residual of 2% and a root-mean-square normalised residual of 3%.Further, the pressure sensors had a maximum error of 2%, and the flow meter that measured the volumetric flow rate to assess the flow velocity had an maximum error of 5%.The experimental set-up is shown in Figure 6.The set-up involved a 100 mm pipe with a long horizontal section to permit the flow to reach a fully developed state, before the flow passed through two pressures sensors separated by a horizontal section of 15 m, two 90 • bends and another horizontal section of 12.85 m.
The numerical solver was first validated against KERS2408, KERS0608 and PARK1, following which, it was used to simulate the set-up shown in Figure 6.Upon simulating the fluid S17 (see Table 1) for a range of velocities, the pressures calculated numerically matched their experimental counterparts at the locations shown in Figure 6, while accounting for the presence of two vertical pipe bends.This provided additional validation of the numerical method proposed in this article.More details on our studies on pipe bends will comprise another article in the future.Figure 7 shows the residuals obtained for one test case.A simulation was considered as converged only after the residuals (not normalised) dropped below 10 −6 (in most cases below 10 −7 ).The sudden increase in the middle of the plot is the result of changing the wall function from the standard Newtonian function [18] to ψ 1 .In either case, the simulation reached convergence but the values of the wall shear predicted were very different (as described in Section 4).

Observations
We present our observations using Figures 8-15.Each figure has two y-axes.The left y-axis represents τ W , while the right y-axis represents the Reynolds number Re R .Therefore, each figure simultaneously shows the relationship of the velocity V on the x-axis with τ W and Re R .For simplicity, Re R is shown with a dashed line with crosses at the relevant data points.
The experimental values of τ W are indicated with grey horizontal lines.We have chosen to use a line to mark the experimental τ W to avoid cluttering the plots with too many markers.This does not imply that a given experimental τ W value holds for all velocities V on the x-axis.Rather, the velocity at which a given experimental τ W is known is made clear with a green error bar on the relevant grey line, representing a 5% deviation from the experimental value, and with an additional magenta bar that indicates a 15% deviation.Finally, if a number of data points for a certain legend entry are missing, this means that the simulation did not converge or the numerically obtained data was very inaccurate and could not "fit within a concise plot".
Figure 8 concerns the case PARK1, which had a very high τ y value.Although not shown on the plot (for the sake of resolution), the use ψ 1 in this case led to an over-prediction in the shear stress.Upon using ψ 2 , the τ W value was predicted more accurately.One notices that κ-remained accurate, whereas the accuracy of RSM1 increased with Re R .
Figures 9 and 10 represent the case S17.Although the τ y value was low, n was also low, indicating a stronger non-Newtonian behaviour [11,12].One notices that as Re R increased, the accuracy of RSM2 with ψ 1 increased, whereas that of κ-with ψ 1 was reduced.At a very high Reynolds number (Re R ∼ 7.5 × 10 5 ; see Figure 10), either RANS model was inaccurate.It is interesting to note that using ψ 2 in such a case (low τ y ) led to an under-prediction in τ W .
At a lower Reynolds number (Re R < 8 × 10 4 ), the accuracy of either RANS model could not be deduced with certainty, although the accuracy of κ-(see Figure 9) seemed to be acceptable.However, when the flow was more turbulent (the regime in which RANS models are generally valid), RSM was more accurate, possibly as a result of its anisotropic formulation that includes more physics than the isotropic κ-.A similar trend is noted in Figures 11 and 12, which pertain to the case S21.The accuracy of κ-was good at low Reynolds numbers (Re R < 7 × 10 4 ; see Figure 11) but diminished at higher Reynolds numbers.RSM2 showed a trend similar to that noted with S17.Further, RSM1 has also been added to elucidate the difference between RSM1 and RSM2 (see Section 2.1).RSM1 remained mostly inaccurate and could, in effect, be replaced by RSM2 or κ-depending on Re R .It is important to note that as Re R increased, in both S17 and S21, the accuracy of both RSM2 and κ-was low, around 6 × 10 4 < Re R < 8 × 10 4 .
KERS0608 represents a case that was mildly non-Newtonian with a high n (for KERS0608, although τ y was low, it was comparable to τ W when Re R was low).As shown in Figure 13, RSM2 was generally inaccurate.Although one sees at Re R ∼ 1.2 × 10 5 that RSM2 with ψ 2 was accurate, the general increasing trend of RSM2 (both ψ 1 and ψ 2 ) suggests that RSM2 was inaccurate at higher values of Re R .Further, RSM1 was much more accurate with ψ 1 , whereas with ψ 2 , it under-predicted τ W .With regard to κ-, ψ 2 is the logical choice at lower velocities, as τ y = 1.88 Pa was comparable to τ W . However at higher velocities, ψ 1 becomes more accurate.Finally, as opposed to cases with lower values of n, κ-is noted to have been more accurate than RSM for KERS0608.
KERS2408 in Figures 14 and 15 is similar to KERS0608 with a high value of n = 0.80 and a τ y value that is comparable to τ W at low velocities (low Re R ).One notices that RSM2 was inaccurate with both ψ 1 and ψ 2 ; κ-seems to be the better choice with ψ 2 at lower velocities and ψ 1 at higher velocities, in keeping with what was observed for KERS0608 (both had a high value of n).
As noted for previously described cases, when Re R increased to values above 1.2 × 10 5 , the accuracy of all the methods was reduced, and the prediction deviated by more than 5% from the experimental value.However, in this case, κ-still remained between 5% and 10% accurate.
We also simulated more test cases from Slatter [13] and observed the trends mentioned above.We have chosen to include only the cases mentioned in Table 1 for succinctness.

Conclusions and Outlook
On the basis of the mentioned observations, one may conclude that a combination of RANS modelling and a specified shear-based wall boundary condition suitable for non-Newtonian flows can provide accurate estimates (within 5% in most cases) of the wall shear stress and the pressure drop through circular horizontal pipes carrying a Herschel-Bulkley fluid.
When used with the proposed specified shear boundary condition (ψ 1 ), the RSM is observed to deliver more accurate results with fluids that have a smaller behaviour index or n ∼ 0.6, alluding to a stronger nonlinear behaviour.On the other hand, κ-seems to work better with fluids that have a value of n greater than 0.8.
Further, the wall boundary condition must be modified (ψ 2 ) to account for the influence of a fluid's yield stress or τ y when it is comparable in magnitude to the wall shear stress or τ W .Therefore, when τ y is large, both RSM and κ-are observed to be accurate when used with ψ 2 , depending on the value of n.It is fair to mention that the accuracy of these numerical methods depends to a great extent on the rheological input.As observed in Slatter [13], models based on rheological inputs measured with laminar flows tend to deviate from experiments for turbulent regimes.It is therefore important to perform a sensitivity analysis of ψ 1 and ψ 2 in terms of the rheological inputs, which will also help to define the envelope within which these equations can be used.Further, the effect of the pipe diameter must be taken into account should these equations be extended to thicker pipes, as it is well known that flow regimes could change with increasing pipe diameter.We aim to address these issues in future publications.
Future work will additionally concern the extension of the specified shear functions for the study of the flow through bends in circular pipes carrying a Herschel-Bulkley slurry and possibly in an unsteady setting to study pulsatile effects.Lastly, using a commercial code adds restrictions to the numerical schemes that a user could use.For instance, upwind schemes used with RANS in ANSYS FLUENT are meant to guarantee numerical stability but may influence the final solution through numerical diffusion.Therefore, we aim to implement the specified shear conditions (ψ 1 and ψ 2 ) in OpenFOAM, which could provide more flexibility in terms of using higher-order (third or higher) upwind schemes or a blend of central and upwind schemes.

Figure 2 .
Figure 2. A schematic of a circular horizontal pipe.

Figure 4 .
Figure 4.The effective mixing region in a pipe carrying a Herschel-Bulkley fluid with yield stress τ y .The unyielding plug represents the region wherein τ < τ y .

Figure 6 .
Figure 6.A schematic of the experimental set-up showing the location of the pressure sensors and the dimensions of the various sections [49] (courtesy: Stichting Deltares).

Figure 7 .
Figure 7. Normalised residuals for one of the test cases.Applying ψ 1 improved the numerical prediction of the wall shear stress.