Particle-laden thin-film flow in helical channels with arbitrary shallow cross-sectional shape

Particle-laden flows in helical channels are of interest for their applications in spiral particle separators used in the mining and mineral processing industries. In this paper, we extend the previous work of Lee, Stokes, Bertozzi (2013) by studying thin-film flows of mono-disperse particle-laden fluid in helically-wound channels of arbitrary centreline curvature and torsion, and arbitrary cross-sectional shape. In the case where the particles are uniformly distributed through the depth of the film, significant analytic progress can be made: the governing equations reduce to a single nonlinear ordinary differential equation, which is readily integrated numerically to obtain the solution subject to appropriate boundary conditions. Motivated by possible application to the design of spiral separators, we consider the effects of changing the channel centreline geometry, the cross-sectional shape and the particle density on the resulting flows and the radial distribution of particles. Our results support the findings in Arnold, Stokes, Green (2016) regarding the effect of channel centreline geometry and cross-sectional shape on flows in particle-free regions. In particle-rich regions, similar effects are seen, although the primary velocity is lower due to increased effective mixture viscosity. \ys{Of key interest, is the effect of channel geometry on the focusing of the particles, for given fluxes of fluid and particles. We find that introducing a trench into the channel cross-section, a feature often used in commercial spiral particle separators, leads to smaller radial width of the particle-rich region, i.e. sharper focusing of the particles, which is consistent with experiments showing that channel geometry influences particle separation in a spiral separator.


I. INTRODUCTION
The work described in this paper is motivated by the aim of improving our understanding of the operation of spiral particle separators, devices used in the mining and mineral-processing industries to separate ores and clean coal. 25,26 They consist of a helically-wound channel (as shown in Figure 1), down which a slurry of crushed rock/ore and water flows, driven by gravity. The complex secondary motion of the fluid, along with the particles' inertia and tendency to settle under gravity, acts to sort the particles across the channel width according to their density, assuming milling has produced particles of uniform size. Thus, particles of different densities can be separated from each other by using vertical splitter plates at the end of the channel to split the flow at appropriate positions across its width. It is known that the effectiveness of spiral particle separators is quite sensitive to the composition of the input slurry, and small changes to the size, density, or concentration of the particles can lead to dramatically impaired separation efficiency. Conversely, small changes to the separator geometry can greatly affect separation efficiency for a given feed slurry. Often, spiral particle separators operate in a regime very close to beaching, where particles are deposited on the channel bottom at the edge of the region occupied by fluid, and thus fail to reach the end of the channel 15 . Spiral particle separators are an established technology; an 1899 US patent 30 describes such a device. Yet they did not see widespread adoption until the late 1940's and 1950's, with the Humphrey's Spiral Concentrator 19,32 . Despite their widespread use, and experimental work undertaken by Holland-Batt 12,13 in the 1970s, 1980s and 1990s, mathematical modelling of the flow in spiral separators has been limited, and the design process remains largely experimental; Holland-Batt 14 gives a summary of the process. However, experimentation with spiral separators has, in general, proven difficult. Estimates of the errors in the measurement of flow velocities can be as high as 30%. Nevertheless, a series of experiments [16][17][18] showed that a complex secondary flow exists, that a fully-developed steady state flow profile appears relatively quickly (within 2-3 turns down the spiral), and that the fluid depth is typically very shallow (1-12mm) relative to the channel width (250-350mm). These last two observations in particular, suggest that, in mathematical models, the steady flow and thin film assumptions can be exploited, resulting in significant simplification of the governing equations. Recently, Boucher et al. 6,7 performed experiments using positron emission particle tracking to obtain more precise information about the motion of particles in spiral separators. Their results allow particles to be tracked in the radial and axial directions, but they do not track the vertical velocity of particles and/or fluid, nor give the fluid depth profile and the vertical particle distribution. Improving mathematical modelling of flow in spiral separators has the potential to aid the design process by providing quantitative predictions of the separation performance of different channel configurations.
Single-phase flows of Newtonian fluid (i.e. excluding particles) in helically wound pipes have been extensively studied by Berger et al. 5 , Germano 10 , Wang 40 and many more. In the case of an open helical channel however, the flow has a free-surface, which makes the problem more complex. Early progress came in the form of empirical formulae, notably from Holland-Batt and Holtham 16 . They performed single-phase flow experiments on spiral particle separators and tried to explain the results by using empirical equations based on a Manning equation near the inner wall and a free-vortex in the main section of the flow. Stokes 35,36 mathematically investigated flow in a helical channel with small centreline torsion and curvature, and with a semi-circular cross-section filled with fluid. Following this, efforts were made to find solutions valid under the more general conditions of non-small curvature and torsion, shallow fluid depth, and with different channel cross-sectional shapes. In particular, exploiting the thin-film approximation, which pertains to flows in spiral separators, has allowed significant progress for rectangular channels 37,38 , although all of these studies required some limits on the channel centreline, be it small torsion, small curvature, or small centreline slope. Subsequent work 1 showed that significant analytical progress was possible even when these assumptions were relaxed in the case of rectangular channels. Results for the more general case of an arbitrary channel cross-sectional shape have recently been obtained 2 .
The mechanics of particle-laden fluids has been a topic of research for many years. Einstein 9 developed a leading order approximation for the increase in effective viscosity of a dilute mixture of particles in a viscous fluid. The second order term was found in the 1970s by Batchelor and Green 4 and Batchelor 3 . For dilute suspensions, the particles have a relatively minor effect on the bulk flow, but as the volume fraction of particles increases, the particles have a larger effect. To describe more concentrated suspensions, models that allow for two-way coupling between the fluid and particle motion are required. Such models can be split into two main categories: discrete and continuum. Discrete models track individual particles, and are very computationally intensive 24 . Continuum models treat the fluid and particles as two separate continuous phases, with coupled equations for each.
Two extensively-studied continuum models for suspension flows are the diffusive-flux model introduced by Leighton and Acrivos 21 and the suspension-balance model of Nott and Brady 29 , which are discussed and compared in some detail by Murisic et al. 28 The diffusive-flux model gives a simple framework to study various phenomena acting on the particle phase, such as Brownian diffusion, settling under the effect of gravity, and shear-induced migration. The suspension-balance model is based on rigorously averaging the mass, momentum, and energy balance equations, is more general, and actually contains the diffusive-flux model as a special case 29 . As shown by Nott and Brady 29 , for 2D channel flow and pipe-flow the suspensionbalance model predicts a plug-like region in the centre of the channel/pipe, in keeping with experimental observation, and overcomes the non-physical infinite viscosity and cusp in the particle volume fraction at the centreline predicted by the diffusive-flux model (which demands that the particle fraction reach the maximum packing fraction there due to zero shear). For other flows, in particular wide-gap Couette flow, the suspension-balance model reduces to the diffusive-flux model. Murisic et al. 27,28 reported excellent agreement between experiments and a formulation of the diffusive-flux model similar to the well-known Phillips et al. 31 model, for free-surface monodisperse particulate flows on inclined planes (neglecting side boundaries), for which zero shear at the free surface does not necessarily correspond to the maximum particle volume fraction and there can be no cusp in the particle volume fraction. More recently Wong and Bertozzi 41 extended the diffusive-flux model to bidensity suspensions on inclined planes, finding good agreement with experimental results. The flows of interest in this present paper differ from these studies of flow on inclined planes in just two respects: the bottom of the channel is not necessarily flat and the channel is curved rather than straight. Therefore, we expect that the diffusive-flux model should provide reasonable predictions for our work and, as a simpler particle-flux model for coupling with a more complex flow, is worth investigation.
The diffusive-flux model has been used in several papers 8,27,28 to study particle laden thin films flowing down straight inclined planes. Murisic et al. 27,28 performed experiments, and observed three distinct behaviours depending on the slope of the channel and the particle volume fraction. The particles can settle to the bottom of the fluid, rise to the surface, or disperse evenly throughout the channel depth. For shallow angles of inclination and low particle concentrations, the particles settle to the bed and the suspending liquid flows over the top of the particles with higher speed. For steeper angles of inclination and higher particle volume fractions, the particles collect near the surface of the film, flowing at higher velocity over the top of the suspending liquid. Between the settled and ridged regimes is the well-mixed regime, where the particles disperse uniformly throughout the film thickness, and the fluid and particle velocities are similar. The diffusive-flux model performed extremely well in correctly predicting the flow behaviour; even where parameter values led to the maximum volume fraction being reached at the free surface (ridged regime) or at the channel bed (some, but not all, instances of the settled regime) the model predictions agreed closely with experiments.
Lee et al. 20 used the diffusive-flux model in a study of monodisperse particle-laden flows in helical channels of small slope and found solutions in the well-mixed regime. In this paper we extend their work to channels of arbitrary centreline curvature and torsion, and arbitrary (but shallow) cross-sectional shape, and demonstrate a method for coupling particle transport with fluid flow in helically-wound channels. Although, unlike in the work of Lee et al. 20 , our model does not dictate the well-mixed regime, we will require this as a simplifying assumption (to be discussed later), which precludes the particle volume fraction from reaching its maximum value and reinforces our expectation that the diffusive-flux particle-transport model is a reasonable choice for the present work. Then, motivated by spiral particle separators in which, assuming all particles in the flow are of the same size, the best separation of particles of different densities is achieved by the focusing of different particles into distinct non-overlapping horizontal bands, we use our model to investigate the effect of channel geometry on the horizontal distribution of particles and, in particular, the focusing of the particles for a given monodisperse feed, where sharper focusing is interpreted to mean a particulate region of smaller horizontal extent for the same total particle and fluid fluxes along the channel. Thus, we aim to use mathematical modelling to gain insight into the critical geometrical parameters that control particle focusing and so increase the understanding needed for improved separation.
Studies on the effects of varying channel cross-sectional shape on resulting flows have been performed for unsteady viscous gravity currents in straight channels, for example by Takagi and Huppert 39 and Golding and Huppert 11 . Likewise, we aim to undertake a systematic study of the impact of the channel geometry on particle focusing. Our model involves a number of governing parameters additional to cross-sectional shape, including slurry composition and channel centreline shape, leading to a large configuration space. Hence, we restrict ourselves to studying three basic types of channel cross-section which, nevertheless, provide a good understanding of how cross-sectional geometry may be used to improve particle focusing, namely rectangular, parabolic, and rectangular with a trench. This paper is structured as follows. We begin by describing the geometry of the spiral channel and define a coordinate system in Section II A. We then introduce the fluid/particle model, and derive two conservation equations in Section II B. We nondimensionalise and scale the system of governing equations in Section III, and present the solution in a special case in Section IV. We present and discuss results in Section V and make conclusions in Section VI. A table summarising the notation is provided in Appendix B.

II. MATHEMATICAL FORMULATION
A. Geometry and coordinate system We consider a helically wound channel with arbitrary shallow cross-section, such as shown in Figure 1. The choice of coordinate system is very important in studying flows in helical geometries. We are interested in the final configuration of fluid and particles near the channel exit and hence study the fully-developed steady-state, helically-symmetric flow. Helical symmetry necessitates the use of a non-orthogonal coordinate system 42 , and we choose to use a body-fitted helicoidal coordinate system as used in Arnold et al. 2 , which is similar to that used by Lee et al. 20 , Manoussaki and Chadwick 22 and Arnold et al. 1 A point x in the fluid domain is represented by its radius, r, from the vertical axis, angular position, β, measured from a reference angle, and a vertical displacement, z, from the channel bottom, as where 2πP is the pitch of the channel centreline, and H(r) is the channel cross-sectional profile. The final parameter required to specify the channel geometry is A, the radius of the channel centreline. Equation (II.1) gives a non-orthogonal coordinate system, and tensor calculus must be used in order to derive component forms of differential operators. Simmonds 34 is a good introduction into the techniques used. The velocity components in this coordinate system are denoted v r , v β , and v z , in the r, β, and z directions, respectively.

B. Particle-laden flow model
We consider a monodisperse suspension of uniformly-sized spherical particles with diameter d. The channel half-width is denoted a, and aδ is a characteristic fluid depth where δ ≪ 1, as the fluid depth is very small relative to the channel width. We require d ≪ aδ so that the particles are small relative to the fluid depth, and the slurry can be treated as a continuum. We will later give a further restriction on d to ensure that the particles remain large enough to avoid the colloidal limit where Brownian diffusion dominates gravitational settling.
Following Phillips et al. 31 , Murisic et al. 27 , and Lee et al. 20 , we assume that the slurry can be treated as a viscous fluid, and denote its velocity, v(x, t). We introduce the local particle volume fraction, φ(x, t), which can vary from zero (no particles) to a maximum value φ m (0 < φ m < 1) which corresponds to the closepacking limit for randomly-arranged spheres, where the mixture can no longer be modelled as a Newtonian fluid. We assume the densities of the particles and fluid are constant, and so the local effective mixture density is given by, where ρ p and ρ f are the densities of the particles and fluid, respectively. Note that the effective mixture viscosity, µ, will vary with φ (i.e., µ = µ(φ)) in a manner which will be discussed later.
We assume the motion of the particles involves two contributions: advection with the mixture velocity v, and motion relative to the mixture due to other physical effects, denoted by the particle-flux vector J. Conservation of mass for the particles and fluid, respectively, then gives where D/Dt = ∂/∂t + v · ∇ is the usual material derivative. We use equations (II.6) and (II.7) to impose conservation of mass for the two phases. Thus, in place of the usual continuity equation for incompressible flow, ∇ · v = 0, we have that v + J is solenoidal.
We obtain an equation for v by considering the momentum balance for the mixture. Assuming the flow is steady, and treating the slurry as a viscous fluid (but allowing its viscosity and density to vary in space) yields the steady Navier-Stokes equations, which can be written as, where p is the pressure, ρ = ρ(φ) is the effective density, µ = µ(φ) is the effective viscosity, S = ∇v + ∇v T is the rate-of-strain tensor, g is the acceleration due to gravity and e z is a unit vector pointing in the vertical direction.
The local effective viscosity, µ, of the fluid is dependent on the particle volume fraction and we use an empirical model to describe this dependence. Note that many such models exist (see Murisic et al. 28 for a discussion), however, where a specific functional form is required, we follow Murisic et al. 27 and adopt the Maron-Pierce equation 23 which is simply the Krieger-Dougherty equation with the power rounded from −1.8 to −2; µ f is the viscosity of the fluid without particles. Various values of φ m , the maximum packing volume fraction, are reported in the literature, and we will use φ m = 0.61, again as in Murisic et al. 27,28 As φ → φ m we have µ → ∞ so that this equation is strictly not valid for φ very near φ m . The Krieger-Dougherty equation was obtained by fitting data for 0.01 < φ < 0.5 31 , while the models of Murisic et al. 27,28 for suspension flows on inclined planes, which utilise the Maron-Pierce equation, have been validated against experiments for feeds with 0.2 < φ < 0.5 and local solids volume fraction as high as 0.6. Should φ anywhere increase above 0.6, the viscosity there will increase rapidly and the flow will go to zero, which can be intepreted as beaching and, consequently, a channel geometry that is not fit for purpose. We note that in the feed to spiral separators, which motivate our study, the total solids volume fraction of particles of all densities is typically 0.25-0.35 25 .
It remains to specify the particle-flux vector, J, which models the movement of particles relative to the suspending fluid. This motion is assumed to occur due to the effects of gravitational sedimentation and shear-induced migration, and hence we use the same form as Murisic et al. 27 and Lee et al. 20 : where K c and K v are diffusion coefficients andγ = 1 2 S ij S kl g ik g jl is the local shear rate (S is the rate-ofstrain tensor and g is the metric tensor). This particle-flux vector is valid when the diffusion due to shear is much larger than the Brownian diffusion of particles in the solvent. It is valid for unidirectional shear 31 , which, although not true in general for the problem of interest here, is the case at leading order. The first term in the expression for J corresponds to sedimentation, the settling of particles on the channel bottom due to the influence of gravity, and the second and third terms correspond to shear-induced migration, the tendency of particles to move away from areas of high shear and high particle volume-fraction (such as near the channel bottom) towards areas of low shear and particle volume-fraction (such as the free surface near the centre of the channel). The function ω(φ) is a hindrance function that models the resistance to motion of a single particle based on the local concentration of particles. There are many different models for the hindrance function, however here we use the simple function ω(φ) = 1 − φ, valid in the presence of shear 33 and used in Murisic et al. 28 and Lee et al. 20 .
We note that in our model, the velocity v is not solenoidal. This changes some aspects of the derivation of the Navier-Stokes equations in our body-fitted coordinate system from that undertaken in Arnold et al. 1,2 , for example, where the usual continuity equation was used to help simplify the equations. The derivation is lengthy and, because of the similarity to Lee et al. 20 and Arnold et al. 2 , is not presented here, however, the full helically-symmetric, steady-state Navier-Stokes and conservation equations are given in Appendix A. We further remark that in other studies using the shear-induced migration model, the velocity has been assumed to be solenoidal. Phillips et al. 31 and Murisic et al. 28 required that the velocity v be divergence-free, and an equation equivalent to Equation (II.7) was also given. Lee et al. 20 gave two conservation equations equivalent to those we have derived above, and also required v to be divergence-free. Thus two equations were added to the clear-fluid model, but only one variable (the particle volume fraction φ), implying that the system of governing equations obtained was overdetermined. The model was only self-consistent if ∇ · J = 0, which is not necessarily true. This is an important difference between these previous models and that presented here.

III. THIN-FILM EQUATIONS
We here nondimenionalise the steady-state, helically-symmetric governing equations given in Appendix A. As defined earlier, a is a representative channel half-width, aδ is a representative fluid depth, and A is the radius of the channel centreline. Since we are concerned with thin-film flow, we assume 0 < δ ≪ 1. Using carets to denote dimensionless quantities, we write where U is a characteristic velocity scale and ρ s = (ρ p − ρ f )/ρ f is the density difference between particles and fluid, relative to the fluid density. Note the new radial direction variable y, measured from the channel centreline. The dimensionless slope of the channel centreline is λ = P/A (the pitch of the centreline divided by its radius), and we define a curvature-like parameter ǫ = a/A. Then gives the local slope of the channel bottom at any point across its width, and Υ = 1+Λ 2 is used for notational convenience. The Reynolds and Froude numbers are defined as, and, as in Arnold et al. 2 , we set, The first of these relationships implies a significant gravity-driven axial flow at the channel centreline y = 0 and, along with the definition of the Reynolds and Froude numbers, gives an equation for the axial velocity scale U , The second relation in Equation (III.4) ensures sufficient cross-sectional flow for an interesting (non-flat) free-surface profile. The scaling of the particle-flux vector J includes the particle diameter d, which must be small, but we need to compare it to the fluid depth scaling, δ, in order to compare terms in the perturbation expansions of the governing equations. Previously we stated the requirement d ≪ aδ so that the particles are small relative to the fluid depth, justifying the continuum modelling approach. If the particles are too small, however, Brownian diffusion becomes a dominant effect and settling due to gravity is negligible. To ensure this does not occur, and as in Lee et al. 20 and Murisic et al. 28 , we require the settling distance along the channel, L settle , of the particles to be small relative to the axial length scale of the channel. We assume that the channel is much longer than it is wide, and set the axial length scale as a/δ. Considering the film thickness, settling velocity, and axial velocity scales in the spirit of equation (2.3) in Murisic et al. 28 , we obtain, For this settling length to be small relative to the axial length scale, we require or (d/a) 2 ≫ δ 4 . Combining this with the condition d ≪ aδ, we obtain, asymptotically as δ → 0, and so we set (d/a) 2 = O(δ 3 ) when comparing terms in the asymptotic expansions of the governing equations. Using Maple we obtain the following system of equations at leading order in δ, where we have, from here on, omitted the carets on all dimensionless variables. The Navier-Stokes equations (A.3)-(A.5) yield the conservation equations (A.1) and (A.2) give and the boundary conditions (A.6)-(A.10) give Although these equations are significantly simplified compared to the full Navier-Stokes equations, they cannot be solved analytically as, in general, the particle volume fraction will vary in both the y and z directions, and hence µ and ρ are functions of both y and z.
Experimental work by Zhou et al. 43 , Murisic et al. 27 and others on particle-laden thin film flows on inclined planes found three qualitatively different behaviours depending on the angle of inclination and initial particle volume fraction. The three regimes are termed settled, well-mixed, and ridged. In the settled regime, particles tend to settle quickly and clear fluid flows more quickly over the top of the particles. In the ridged regime, particles move towards the surface of the film, move more quickly than the surrounding fluid, and amass at the fluid front. Finally, in the well-mixed regime, the fluid and particle velocities are roughly the same, and the particles do not tend to congregate at the bottom or top of the fluid, but disperse uniformly throughout the film thickness. As is readily seen from the figures in Murisic et al. 27 , for a plane at a given angle of inclination, the settled regime applies for smaller initial particle volume fraction and the ridged regime for larger initial particle volume fraction, with the well-mixed regime applying at intermediate values of the initial particle volume fraction, between the settled and ridged regimes. The range of particle volume fraction associated with the settled regime decreases as the inclination angle increases.
In Lee et al. 20 , it was found that the only solutions that their thin-film governing equations could admit were of the well-mixed type. However, as discussed in Section II B, we believe their system of equations was overdetermined, and unnecessarily constrained the particle volume fraction φ to be independent of depth. For our model, this restriction need not necessarily apply. Nevertheless, motivated by the existence of well-mixed flows down inclined planes, we here seek solutions in helically-wound channels in the wellmixed regime, by assuming φ = φ(y) so that the particle volume-fraction is independent of depth. This has the benefit of making the governing equations significantly more tractable, enabling further analytic progress. This modelling approach is also reasonable given the present lack of experimental data on particle distribution in the vertical direction with which to compare, and the need for a model to aid separator design that gives understanding of the radial distribution of particles. Spiral separators work because of a radial separation of particles of different densities, while vertical distribution appears not to be important. We thus consider whether there is a region in the parameter space in which such solutions exist and, if so, whether they differ significantly from those found by Lee et al. 20 for channels of small slope. The more general case where the particle volume-fraction varies with both depth and radial position is left for future work.
The steps outlined below are broadly similar to Lee et al. 20 , although the equations are different. With φ = φ(y), both the density and viscosity become functions of y only, and this allows the thin-film governing equations to be solved by sequential integration. Equation (III.11) can be solved directly to give and Equation (III.10) gives Note this equation is used to set the first of the scalings in Equation (III.4) to ensure sufficiently high axial velocity at the channel centreline. Now Equation (III.12) can be integrated with respect to z to give J z = C(y) for some C(y), and the boundary condition J z | z=0 = 0 gives J z = 0, so gravitational settling is balanced by shear-induced migration. We first note that at leading order in the thin-film limit, the shear rate isγ Then, writing the z-component of the particle-flux vector explicitly gives an equation involving φ, µ, and u, which can be rearranged as Now ∂φ/∂z = 0 by assumption, we use (III.10) to replace µ∂u/∂z, and substitute the hindrance function This is a quadratic equation for φ(y) with solution (IV.7) The particle volume fraction given in Equation (IV.7) is a function of Λ(y) (the local slope of the channel bottom at position y) and must satisfy 0 ≤ φ < φ m = 0.61 in the particle-rich region in order for our model to be valid. Figure 2 shows a typical profile for φ in terms of Λ(y). Note that φ(y) → 1 as Λ → 0 (corresponding to y → ∞) and as Λ → ∞ (corresponding to y → −1/ǫ + ), and takes a minimum value at Λ = 1 (where y = (λ − 1)/ǫ). In the limit Λ 2 → 0, when Υ = 1, (IV.7) is equivalent to Equation (11) in Lee et al. 20 . This is to be expected; their model was derived in the small-slope limit λ 2 → 0, when Λ 2 is also small (so long as ǫ is not too close to 1). Their φ profiles were monotonic increasing, but now allowing a region of higher shear at the inside of the channel (with non-small λ), a higher concentration of particles can be maintained there. Increasing the local fluid velocity therefore allows for a higher concentration of particles to be sustained in the well-mixed regime.
To ensure φ(y) < φ m wherever particles occur, the particle-rich region must lie fully within the open interval (y b , y t ), the end-points of which satisfy φ(y b ) = φ m and φ(y t ) = φ m . We find where and, for example, χ = 1.8421 with K c = 0.41, ρ s = 1.5 and φ m = 0.61. This restricts our results somewhat, which will be discussed further in Section V. Finally we note that there is a maximum ρ s above which φ(y) ≥ φ m for all y. To have any particles in the flow, we require points y b and y t to exist and hence require χ > 1. This is satisfied as long as ρ s < 12.04 (with the given K c and φ m ). By way of comparison, glass particles with ρ s ≈ 2.5 and ceramic particles with ρ s ≈ 3.8 have been used in experiments 28 .
Having determined φ, we next find the radial velocity, v, by integrating Equation (III.9) twice, yielding, v = − Re 3 ǫρ 3 Λ 2 Fr 4 (1 + ǫy)µ 3 Υ 5 z 120 (IV.10) Equation (III.13) may be written as and integrating with respect to z, rearranging using Liebniz's integral rule, and applying boundary conditions, gives where C is a constant which effectively controls the net radial flow at any position y. We set C = 0 so that there is no net flow into or out of the fluid domain, to obtain, We now introduce a streamfunction ψ, defined by Using (IV.15) and requiring ψ = 0 on the channel bottom z = 0, we find, The second equation in (IV.16) can be used to find the vertical velocity w (not given due to its complexity).
As we are dealing with particle-laden flow, two fluxes are relevant: the total mixture flux and the flux of particles. The total (dimensionless) flux down the channel Q (scaled by δa 2 U ) is given by where y * l and y * r represent the boundaries of the particle-rich region of the flow, which do not necessarily correspond to the whole fluid domain. The particle-rich region must be contained in the fluid region, so we must have y l ≤ y * l < y * r ≤ y r . The negative signs in the two fluxes given above are due to the axial coordinate direction pointing up the channel, meaning u is always negative.

V. RESULTS
Having determined that there is a region of the parameter space in which there is a solution to the thin-film system of equations for particles uniformly dispersed in the vertical direction, we now present and discuss some particular cases. For a given channel geometry, ǫ, λ, and H(r), we have a set of 8 (not all independent) parameters: {y l , h l , y r , h r , Q, y * l , y * r , Q p }. To specify a solution, we choose 5 of these parameters: y l or h l , y r or h r , Q, y * l , and Q p . From these, we find the remaining parameters: y l or h l , y r or h r , and y * r . For a channel with a vertical wall, we specify the value(s) y l and/or y r , and for a channel with a curved bottom, we set h l = 0 and/or h r = 0. Finding the unknown parameters involves a numerical root-finding problem to find the extents of the partical-laden and clear fluid regions for the chosen total and particle fluxes.
Different choices of parameters could be used to specify a solution. For example, we could specify y * l , y * r and three out of the following group of four parameters: {y l , h l , y r , h r }. The parameter not picked, as well as the fluxes Q and Q p could be calculated by solving the governing equations. This choice would be simpler than specifying fluxes because it can be directly calculated by solving the free-surface differential equation in the particle-rich and particle-free regions. Motivated by spiral particle separators, which are often described in terms of the quantity of material that they can process per hour, we choose to specify our solutions in terms of the fluid and particle fluxes. In addition, this choice enables an examination of how channel geometry affects the radial focusing of particles for given fluid and particle fluxes.
By considering spiral particle separators, we expect that, in steady state conditions, the particle-rich region of the flow should be near the inside wall of the channel, and extend part or all of the way to the outside channel wall. The volume fraction of particles φ is dependent on radial position y, and the flux of particles, together with φ(y), determines how wide the particle-rich zone will be. Lee et al. 20 provided an intuitive explanation as to why this configuration of particles is expected to be stable. For completeness, we briefly recapitulate their argument here.
Let us first assume that the particles collect near the inside channel wall, and there is a clear fluid region in the outer region of the channel cross-section, as in Figure 3a (top image). The outward flow near the free surface will carry particles radially outwards into the particle-free region (second image), where, because the local slope of the channel bottom will be smaller and the particle volume fraction will be low, the settled regime will apply and the particles will rapidly equilibrate to the channel bottom (third image). The equilibration will be fast due to the assumption of short settling length. The secondary flow, which is radially inwards near the channel bottom, then carries the particles back into the particle-rich region (bottom image). So in practice we would expect that if such a configuration of particles was reached, it would remain stable. Now assume that there is a configuration with a particle-free region near the inside wall and a particlerich region further outwards, as in Figure 3b (top image). On the boundary between the particle-rich and particle-free regions, the flow near the channel bottom will carry particles into the particle-free region (second image). Again, the settled regime will apply due to the low volume fraction of particles so these particles will remain near the channel bottom and continue to be carried to the inside of the spiral by the secondary flow. At the same time, in the outer particle-rich region, particles will settle into the vacated area near the channel bottom (second figure). The process will continue over time (third and bottom images) and the particle-rich outer region will first enter the settled regime, when the particle volume fraction reduces sufficiently, and finally disappear. This configuration of particles is therefore not stable.
From this we deduce that, over time, the particles will migrate to a region next to the inside wall of the channel. Hence, we search for steady solutions with the particle-rich region near the inside channel wall and the particle-free region near the outside wall. This means we set y * l = y l . We note that this configuration is consistent with the known fact that heavy particles migrate to the inside wall of a spiral separator.
We solve the governing equations separately in the particle-rich and particle-free regions, and couple them by requiring that the free-surface be continuous at the boundary between the two regions. At the interface between the particle-rich and particle-free regions of the flow, the fluid properties, the streamlines (with the exception of the free surface), and the axial flow contours are discontinuous. In practice there will be a thin boundary layer between the two regions, which is not captured by our leading order model, across which quantities change continuously and which would remove the non-physical cusp in the free surface that is seen in our results at the interface.
We will be particularly interested in considering how well a channel can focus the particle-rich region of the flow into a narrow band of the channel cross-section. Certainly, for polydisperse slurries, separation of particles of different densities requires the focussing of each particle species into non-overlapping bands in the cross-section. This suggests that minimising the band width is advantageous. Below we will consider two basic channel cross-sections, rectanglar and parabolic, the first as the simplest case that has vertical side walls at which the fluid depth must be determined, and the second as the simplest case that requires determination of the horizontal extent of the fluid with a zero fluid depth at each side. We then examine the effect of a trench in the channel bottom. As will be seen, this choice of channel cross-sectional geometries provides a good understanding of how cross-sectional geometry may be used to improve particle focusing.
Throughout this section, we use dashed streamlines for the particle-laden section of the flow, and solid streamlines for the particle-free fluid. A solid vertical line indicates the boundary between particle-free and particle-rich regions at y = y * r . In all plots we set K c = 0.41, which has been empirically determined from experiments 31 , and φ m = 0.61 as previously discussed. Unless stated otherwise, we will use ρ s = (ρ p − ρ f )/ρ f = 1.5.

A. The effect of centreline geometry
Changing ǫ and λ in particle-free flows has been studied in detail in Arnold et al. 1,2 , and we now seek to understand the effect of ǫ and λ on particle-rich flows. Figure 4 shows the effect of changing ǫ on the resulting flow in a channel of rectangular cross-section. The total flux and particle flux are the same in each plot, Q = 0.15 and Q p = 0.003. The addition of the particlerich zone causes significant differences from clear-fluid solutions, with a sharp change in the free-surface slope at the interface between the particle-rich and particle-free regions, separate secondary flow profiles in the two regions, and much lower axial velocity in the particle rich region. The streamlines in the particle-rich and particle-free regions are not plotted at the same contour levels, as the magnitude of the streamfunction in the particle-free region is larger than in the particle-rich region. There are some similarities with the row shows the free-surface profile and stream-function of the secondary flow in the particle-rich (dashed lines) and particle-free (solid lines) regions. The middle row shows the particle volume fraction φ(y). The bottom row shows contours of the axial velocity u in both the particle-rich (dashed lines) and particle-free regions (solid lines). The scaling of the streamlines and axial-flow contours differs between the regions, with the magnitude of the flow velocity being significantly smaller in the particle-rich region.
clear-fluid solutions from Arnold et al. 1,2 ; increasing ǫ tends to push the mixture towards the inside of the channel, increasing the depth there and decreasing it at the outside of the channel. The particle volume fraction φ(y) is also shown in Figure 4. In all three cases, the particle-rich region of the flow has a high volume fraction of particles, near to φ m = 0.61, but everywhere below 0.6, which, as discussed earlier, we interpret to mean no beaching. The effect of ǫ on the location of the interface between the particle-rich and particle-free regions of the flow (y * r ) is complex. The value of ǫ modifies the local slope Λ(y) = λ/(1 + ǫy) at each position y and, therefore, the volume fraction of particles at each y changes in a nonlinear manner according to the relation (IV.6) plotted in Figure 2. As ǫ increases, y * r initially decreases, but as ǫ increases further, y * r starts increasing again. Depending on the slope of the channel centreline λ, we can observe three behaviours: y * r decreasing as ǫ increases (which occurs for large λ's), y * r increasing as ǫ increases (which occurs for small λ's), and y * r decreasing then increasing as ǫ increases (which occurs for intermediate λ's). In Figure 4, y * r decreases as ǫ increases, but if ǫ were increased enough, y * r would eventually start to increase. Figure 5 shows the dependence of y * r on the channel centreline geometry (ǫ and λ), within the valid portion of this parameter space, for a channel of rectangular cross-section with the same total and particle fluxes as used for Figure 4 (Q = 0.15, Q p = 0.003). From (IV.8) and noting that y * ℓ = y ℓ (where y ℓ = −1 for the rectangular channel), the valid part of the (ǫ, λ) parameter space is given by which is the coloured region; outside of this coloured region there is no solution to the problem that maintains φ(y) < φ m in the particle-rich region. This plot helps to make sense of the three behaviours noted above. We see that, for fixed ǫ, there can be a decrease in y * r as λ increases, followed by an increase in y * r as λ increases further. Further useful data that can be obtained from this plot are the values of ǫ and λ for which the particle-rich region is smallest, i.e. y * r is minimised. As mentioned earlier, reducing the radial extent of the particle-rich region could potentially be beneficial in spiral particle separators. For ρ s = 1.5, the case considered in Figure 5, the radial extent of the particle-rich region is minimised when ǫ is maximised. This agrees with our physical intuition, the larger ǫ, the more gravity dominates inertia and, hence, particles (and fluid) tend to aggregate near the inside channel wall.
Also shown in Figure 4 are contours of the axial velocity. Axial velocity is always zero at the channel bottom, but is non-zero on the side walls. This is a result of our thin-film scaling, with does not allow us to impose boundary conditions on the side walls. The only condition we can impose is that there is no net flux through the side walls. In practice we would expect thin boundary layers to form at the side walls, as well as at the interface between the particle-rich and particle-free regions where the fluid properties and velocity are not continuous. The axial velocity contours show that the velocity increases with height above the channel bottom, but there are other geometric effects that move the point of maximum axial velocity. For example, in the particle-free region, the point of maximum axial velocity is at the outer wall for ǫ = 0.2 and ǫ = 0.4, but moves to near the middle of the channel for ǫ = 0.6. Figure 6 shows the effect of changing λ on the resulting flow in a parabolic channel with H(y) = y 2 and fixed ǫ = 0.25. The particle flux and total mixture fluxes are the same for all three plots, Q = 0.08 and Q p = 0.001. There is little qualitative difference in the solutions except for the change in the location of the interface between the particle-rich and particle-free regions. The axial velocity contours do change slightly, with the location of maximum axial velocity moving outwards as λ increases up to λ = 1. Figure 7 shows The top row shows the free-surface profile and stream-function of the secondary flow in the particle-rich (dashed lines) and particle-free (solid lines) regions. The middle row shows the particle volume fraction φ(y). The bottom row shows contours of the axial velocity u in both the particle-rich (dashed lines) and particle-free regions (solid lines). The scaling of the streamlines and axial-flow contours differs between the regions, with the magnitude of the flow velocity being significantly smaller in the particle-rich region.
values of y * r in the valid region of the (ǫ, λ) parameter space as given by (V.1). Note that, the solution process for this channel geometry differs from that for the rectangular channel, in that we set h ℓ = h r = 0 and determine the free-surface contact points y ℓ and y r , and the location of the fluid-particle interface y * r . Hence y * r = y r varies with ǫ and λ. Not only is the effect of λ relatively minor, but the effect of ǫ is even smaller, with the contour levels nearly horizontal in much of the parameter space. For channels with moderate steepness (λ less than roughly 2 or 3), the extent of the particle-rich region is relatively insensitive to both ǫ and λ. Compared to the rectangular channel case Figure 5, y * r is less sensitive to λ for parabolic channels. The centreline geometry that minimises the radial extent of the particle rich region is ǫ ≈ 0.7 and λ ≈ 0.7.
A comparison of Figures 4 and 6 suggests that the shape of the channel cross-section has a more significant impact on the secondary flow and particle volume fraction distribution than do the parameters ǫ and λ. The parabolic channel allows less freedom for the secondary flow and free-surface profile to change than in the rectangular channel case. Note too that the no-slip boundary conditions are satisfied on the entire surface of the parabolic channel which, unlike the rectangular channel, has no vertical walls on which the boundary conditions must be relaxed. In the next section, we will further explore how the channel shape can be changed in order to reduce the radial extent of the particle-rich region.
B. The effect of a trench in the channel bottom Motivated by spiral particle separators, which often feature a trench near the inside wall, we consider a rectangular channel with a trench. The trench is a quartic polynomial that joins smoothly to the flat channel bottom, with maximum depth k and half-width b. The equation of the channel bottom is, and the area of the trench is (16/15)kb. Figure 8 shows the results for such channels with ǫ = 0.3, λ = 0.5, b = 0.3, Q = 0.184, Q p = 0.0035 and k = 0, 0.15, 0.3. The particle-rich region contracts as the trench deepens, and the effect is significant. For the no-trench case, the particle-rich region extends over half the width of the channel (up to y = 0), but with a trench of depth 0.3 and width 0.6, the interface between the particle-rich and particle-free regions is at y * r ≈ −0.35, i.e. particles occupy roughly one third of the width of the channel. The free-surface profile does not change significantly as the trench depth increases, (as seen for clear-fluid flows in Arnold et al. 2 ), but the trench means there is more cross-sectional area near the inside wall, and so more particles can fit within a certain distance from the inside wall of the channel. Since the depth is increased over the trench, the axial velocity (also shown in Figure 8) is increased in this region and, hence, the flux of particles there is magnified. This effect is important in decreasing y * r and, hence, the width of the particle-rich region. The trench has little effect on the secondary flow in the particle-free region, but does entrain a rotating cell of particle-rich fluid. Whilst we expect boundary layers near the vertical channel walls which would cause all the streamlines to connect up and form closed curves (discussed in detail in Arnold et al. 1 ), the trench moves the centre of rotation radially outwards from the inside wall. Figure 9 shows that changing the width of the trench while holding the depth constant has a similar effect to increasing the depth while holding the width constant, namely reducing the extent of particles y * r compared to the rectangular channel with no trench. The trench in Figure 9c has a larger area than any in Figure 8 and is large enough that the particle-rich region does not fill the trench. It is possible to create a very wide trench which is mainly filled with particle-free fluid. This is a poor configuration as the benefit of the trench in focussing the particles is reduced once the deepest part of the trench falls outside the particle-rich region. For a given fluid and particle flux, the optimal trench width appears to be the width such that the edge of the trench roughly coincides with y * r , the edge of the particle-rich region. The channels in Figure 8 and 9 feature trenches with different areas; comparing trenches with the same area is also interesting. Figure 10 shows the effect of changing the width of the trench on y * r , whilst keeping the total area of the trench constant. It is clear that decreasing the width of the trench (and hence increasing the depth) decreases y * r and, hence, focusses the particles more effectively. The specific shape of the trench is not of critical importance for particle focussing. A similar trench using the cosine function was tested, with very similar results. The results obtained when adding a trench to the parabolic channel were also qualitatively similar, suggesting that the effect of adding a trench to a channel doesn't depend significantly on the shape of the channel, and the mechanism of increased depth allowing increased axial velocity adequately explains the results.
C. The effect of particle density Figure 11 shows a rectangular channel, with ǫ = 0.2, λ = 0.7, Q = 0.192 and Q p = 0.0018, carrying particles of increasing density. We first note that ρ s is not the actual particle density, but is defined in terms of the particle and fluid densities as ρ s = (ρ p − ρ f )/ρ f . The effective nondimensional mixture density is ρ = 1 + φρ s , so increasing ρ s increases the density of the mixture for a given concentration φ.  Figures 8b and 8c, respectively, while that in (c) is larger. The scaling of the streamlines differs between the regions, with the magnitude of the flow velocity being significantly smaller in the particle-rich region. As ρ s increases, the volume fraction φ(y) that leads to a well-mixed regime, given by Equation (IV.7), increases. Both ρ and µ (see (II.9)) increase with φ, with µ increasing very rapidly as φ → φ m . The axial flow velocity u (Equation (IV.2)) includes a factor ρ/µ, and if this factor decreases, the axial velocity will also decrease. For the full range of parameters considered in this section (ǫ and λ restricted so that φ(y) < φ m everywhere in the channel), µ grows more rapidly that ρ and ρ/µ decreases as ρ s increases. If the axial velocity decreases the particle-rich region, hence y * r , must grow to accommodate the same particle flux down the channel. This is shown in Figure 11.
The solutions shown in Figure 11 are quite sensitive to ρ s . As ρ s increases, the region in the ǫ-λ plane where φ(y) < φ m everywhere in the channel gets smaller, so that the choice of centreline geometry becomes increasingly restricted. In these cases, φ(y) is relatively large throughout the particle-rich region of the channel, so we are close to the limit of validity of the model equations where µ → ∞ as φ → φ m and the Maron-Pierce equation (II.9) ceases to be a valid representation of the viscosity. As discussed earlier we interpret this to imply beaching and a channel geometry that is not fit for purpose, from which we conclude that our model is only valid for particles with a density ρ p that is not too much larger than the fluid density ρ f .

VI. CONCLUSIONS
We have derived a system of equations (A.3)-(A.5) governing particle laden flow in helically-wound channels with arbitrary centreline pitch and radius, and arbitrary shallow cross-sectional shape. The thin-film scaling was used to give a much simpler system of equations (III.9)-(III.11), although these were still too complicated to be solved analytically. Motivated by available experimental results, the specific spiral particle separator application which relies on a radial separation of particles of different densities, and a previous study 20 , we restricted our attention to finding solutions with uniform particle volume fraction in the vertical direction. This assumption allowed significant analytic progress to be made, with only the free-surface differential equation requiring numerical solution. The results we have obtained are entirely consistent with current knowledge on spiral particle separators.
The solution to the simplified system of equations provides insights into the effect of particles on the fluid flow. The basic features of flows in helical channels that were observed and discussed in Arnold et al. 1,2 still apply to particle-laden flows, but the change in the density and viscosity of the fluid (due to the presence of particles) changes the balance between gravitational and inertial effects. Increasing ǫ or λ tends to push the fluid closer to the inside wall of the channel, and vice versa.
We have considered two basic channel cross-sections, rectangular and parabolic. The first is the simplest case that has vertical side walls at which the fluid depth must be determined for given total and particle fluxes, and the second is the simplest case that requires determination of the horizontal extent of the fluid by finding the contact points of the free surface at which the fluid depth is zero. At a vertical side wall the asymptotic model requires zero net flux through the wall but the no-slip boundary conditions cannot be enforced and, in practice, there is a boundary layer near to the wall that is not resolved by the model. For channels with no vertical walls, so that the fluid depth is zero at points of contact of the free surface with the channel bottom, the no-slip boundary conditions are enforced along the entire channel surface. Similar qualitative behaviour was found for rectangular and parabolic cross-sections. We have also shown the effect of a trench in the bottom of an otherwise rectangular channel and examined the effect of trench depth and width. Although not shown, the effect of a trench in an otherwise parabolic cross-section is similar to that seen in the rectangular case.
Our model provides insight into the radial distribution of particles, in terms of a depth-independent particle volume fraction, and the focussing of the particles within the channel cross-section. For a given flux of particles and total mixture flux, we have shown that the particle-rich region can be focussed into a smaller region of the channel cross-section by changing the channel centreline pitch and radius. The particle rich region can also be focussed by introducing a trench near the inside wall of the channel, essentially increasing the cross-sectional area there. Commercially available spiral particle separators often feature deeper trenches near their inside wall, and our modelling is able to show the benefit in reducing y * r that arises from them. A reported benefit of the trenches in commercial separators is that particles settle into them and are trapped there, however our assumption that particles are evenly distributed in the vertical direction means we cannot investigate this phenomenon. Our results were relatively sensitive to particle density, something that is also seen in practice, as spiral particle separators are known to be sensitive to the composition of the input slurry.
There is significant scope for further research into flows of particle laden fluids in helical channels. It is not certain that the particle volume fraction in spiral particle separators is independent of depth, and relaxing this assumption in our model would give more widely-applicable results. Numerical solutions of the thin-film system of equations (III.9)-(III.13), with φ = φ(y, z) would allow a wider range of channel centreline geometries and flow regimes to be studied. Non-helically-symmetric simulations of the thin-film Navier-Stokes equations in the spirit of Murisic et al. 28 could provide significant insight into particle and fluid motion in spiral separators, for example by studying the transition to equilibrium configurations.
As mentioned in Section V C, the particle volume fraction can become very large, and the Maron-Pierce equation for the viscosity may not be realistic. The very large viscosity it predicts could potentially affect the asymptotic expansion of the governing equations. The model used in this chapter treated the slurry as a Newtonian fluid, and at very high particle volume fractions, this is not necessarily the case, and the particle-rich fluid may be better modelled as non-Newtonian or even as a granular material. Nevertheless, very large particle volume fraction is likely to be associated with beaching and consequently indicates that the spiral separator design is not appropriate to the given feed. Thus our model may be a useful design tool for determining channel geometries that minimise or avoid occurrences of beaching.
The shear-induced migration model used in this paper is just one particle-laden flow model that can be considered. Different models (such as the suspension balance model of Nott and Brady 29 ) could be investigated to see whether they would predict different behaviour. Our model is for monodisperse slurries, with a single particle species, however spiral particle separators are intended to separate different species of particles. Models for flows of slurries consisting of two or more species of particles are less well understood than monodisperse slurries, and there is significant scope to improve modelling of such complex mixtures. Combining a recent model for bidensity suspensions 41 with our equations for the helical geometry is seen as a useful next step towards understanding particle separation in curved channels. simplify the divergence of the rate of strain tensor, and secondly, there are new terms arising from ∇µ · S, The kinematic condition is v · n = 0 at the free surface F (r, z) = 0, and gives v r F ,r + v z F ,z = 0. (A.10)