An Eulerian-Eulerian formulation of suspension rheology using the finite volume method

Accepted Manuscript

An Eulerian-Eulerian Formulation of Suspension Rheology Using the Finite Volume Method N.J. Inkson, D. Papoulias, M. Tandon, V. Reddy, S. Lo PII: DOI: Reference:

S0377-0257(17)30010-1 10.1016/j.jnnfm.2017.05.002 JNNFM 3896

To appear in:

Journal of Non-Newtonian Fluid Mechanics

Received date: Revised date: Accepted date:

9 January 2017 4 May 2017 8 May 2017

Please cite this article as: N.J. Inkson, D. Papoulias, M. Tandon, V. Reddy, S. Lo, An Eulerian-Eulerian Formulation of Suspension Rheology Using the Finite Volume Method, Journal of Non-Newtonian Fluid Mechanics (2017), doi: 10.1016/j.jnnfm.2017.05.002

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • An Euler-Euler two phase model for particles in a Newtonian liquid • Mixture constitutive model of Morris and Boulay is reformulated by splitting the stress

CR IP T

• Successfully predicts correct particle volume fractions and velocity profiles of NMR experiments

AC

CE

PT

ED

M

AN US

• Discussion of anisotropic normal stress difference effects

1

ACCEPTED MANUSCRIPT

An Eulerian-Eulerian Formulation of Suspension Rheology Using the Finite Volume Method

CR IP T

N. J. Inkson, D. Papoulias, M. Tandon, V. Reddy and S. Lo. Siemens PLM Software, Trident House, Trident Park, Basil Hill Road, Didcot, Oxfordshire, OX11 7HJ, United Kingdom [email protected]

1

AN US

May 8, 2017

Abstract

CE

PT

ED

M

We have adapted a well-known constitutive model formulated by Morris and Boulay[32] that describes the stress tensor for a mixture of particles in a Newtonian liquid into a two-phase finite-volume solver. The two-phase model treats each phase with a separate continuity and momentum equation that splits the stress of the constitutive model between the two phases, with the particle pressure applied only to the particle phase in addition to source terms related to the drag of particles in the fluid. We compare the resulting model using a variety of NMR experiments, considering the flow of neutrally buoyant monodisperse particles at low Reynolds number in the following geometries: 2D simple shear flow in a Couette device, axisymmetric and 3D pressure driven flow down a pipe and in a 2D and 3D asymmetric channel bifurcation. We find excellent agreement with velocity and volume fraction profiles in all of the comparisons with experimental data without much further numerical fitting beyond the original suspension stress model.

2

Introduction

AC

Suspensions of particles in liquids have been of great industrial importance to transport large quantities of solids through pipes; for example sand, polymer and waste slurries. Suspensions are also involved in many industries such as paints, pharmaceuticals and foods, they are also ubiquitous through nature in the form of blood. Colloidal suspensions (with particles of size . 1µm) have been studied since the early days of alchemy, in particular gold colloids which were investigated in the nineteenth century by Faraday who made note of their properties to disperse light and their great stability over time. Research by Thomas Graham (1861)[13] established two distinct classifications of a variety of materials that could diffuse through a membrane as either crystalliods or colloids. Research gathered pace in

2

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

the twentieth century where Robert Brown, discovered Brownian motion which was quantified by Albert Einstein [10] and Paul Langevin[24] into stochastic differential equations determining the motion of the particle. Further pioneering research on suspensions followed by Batchelor and Green[4, 3, 17] where they considered the micro-hydrodynamics occurring when spherical particles meet in flow and made corrections to Einstein’s viscosity dependence on particle volume fraction. For particles of larger size than colloidal particles, of which we are concerned with in this work, we can ignore the role of Brownian diffusion as it is negligible. Relative (mixture) viscosity models were developed for suspensions of particles in a liquid by Krieger and Dougherty (1959)[22] where the relative viscosity depends only upon the volume fraction of the particles. As the concentration of particles increases, an exponential increase in viscosity occurs as the particles move towards contact causing moments of hydrodynamic lubrication forces to become large. Eventually jamming occurs at the critical volume fraction (maximum packing). Shear-migration of particles was first observed in pipe flow by Karnis et. al. in 1966[21] by noting that the profile of particles was not parabolic as expected for a Newtonian fluid with matched material properties. Eckstein et. al. (1977)[9] were able to measure the shear-migration phenomenon by using radioactive tracer particles in Couette flow. The shear-induced migration effect arises from the spatial hindrance of particle paths when the concentration of particles is large. As the concentration increases then stronger interactions occur from hydrodynamic forces as two particles approach. In the 1980s Acrivos’ team experimentally studied shear-induced migration extensively and found that a settled bed of particles can re-suspend through this mechanism in shear flow, and noticed that particles tended to migrate from regions of high shear rate towards regions of lower shear rate, for instance in the Couette device the particles move away from the rotating wall[12, 25]. Phillips (1992)[37] built upon this concept of diffusion of particles in shear and by considering the balancing of diffusive fluxes developed the “Diffusive flux model”. They confirmed results using nuclear magnetic resonance (NMR) imaging to measure the particle profiles across the Couette cell over time. The diffusive flux model possesses a simple form, however it predicted a volume fraction profile that showed a high peak near the centreline in flow of a suspension in a pipe that did not match experimental data due to a discontinuity in volume fraction where the shear-rate is zero. Subsequent modifications to the theory have been applied to solve this problem [1]. Nott and Brady (1994)[33] further developed the idea of particle diffusion with the “Suspension balance model” which successfully reformulates the problem in terms of the stress tensor and has the concept of suspension temperature in common with other granular theories. It has been shown that the suspension balance model reduces to the diffusive flux model in simple pipe flows. Later theory by Morris and Boulay[32] builds upon this work using the stress tensor as the sole mechanism for the flux.

3

Theory

The particle stress model of Morris and Boulay[32] derived from the “suspension balance model” of Nott and Brady[33] has been adapted into an Euler-Euler

3

ACCEPTED MANUSCRIPT

3.0.1

Relative Viscosity

AN US

CR IP T

two-fluid[19] formulation for multiphase flows and applied to Siemens PLM R Software’s STAR-CCM+ finite volume software. The model was adapted to split the continuity and momentum equations following the Euler-Euler theory of Ishii[19] and resulted in successfully predicting particle flow in suspensions through pipes in addition to modeling the high pressure drops in stabilized crude-oil and water emulsions due to the relative viscosity with the Krieger and Dougherty[22] and Morris and Boulay models[32]. The application to emulsions in laminar pipe flow was published recently by Inkson. et. al.[18]. Further work has completed the formulation of suspension rheology and been applied to simulate many experiments in order to comprehensively validate the method. The work of Morris and Boulay[32] goes beyond considering purely shear viscosity and takes into account normal stress contributions that causes shear induced migration of particles in shear flow and was shown to give equivalent results to the suspension balance approach. This was performed entirely through the constitutive model and for this reason we feel this approach is superior and easier to implement than previous diffusion techniques. Recent research has inferred a cross-over between the regime of particle kinetic energy dominated granular models, and viscosity dominated rheology models by Boyer et. al. (2011)[5].

M

In suspension rheology, the dimensionless relative (effective) viscosity is required to describe the numerics of the mixture viscosity which tends to infinity as the dispersed phase reaches the limit of maximum packing. The relative viscosity is defined as η , (1) ηr = µc

AC

CE

PT

ED

where η is the mixture viscosity and µc is the continuous phase (Newtonian) viscosity. The relative viscosity prescribes how the mixture viscosity increases as more particles are interacting with each other, along with hydrodynamic forces present in the fluid between the particles. As the maximum packing of the particles is approached networks of contact form between the particles and the mixture jams [6], the relative viscosity tends to infinity essentially encapsulating a yield stress behaviour, although in these relatively simple models the actual value of the yield stress does not directly enter the equations. There are many models for relative viscosity with a similar form, one of the earliest being the Krieger-Dougherty model [22] which can be written as ηr (φ) =

 −[η]φm φ , 1− φm

(2)

where [η] is the intrinsic viscosity, [η] = 2.5 for spherical particles and φm is the maximum critical packing fraction (φm = 0.645) for the random close packing of mono-disperse hard spheres. In the Morris and Boulay model[32] the relative viscosity is defined  −1  2  −2 φ φ φ ηr (φ) = 1 + 2.5φ 1 − + Ks 1− , φm φm φm 4

(3)

ACCEPTED MANUSCRIPT

this form ensures the Einstein limit is retained as φ tends to zero[10]. The “normal relative viscosity” is defined  2  −2 φ φ ηn (φ) = Kn 1− . (4) φm φm

3.0.2

CR IP T

which vanishes in the limit of φ → 0. Ks and Kn are the shear and normal contact strengths and were found to be approximately 0.1 and 0.75 respectively by comparing with experiment. In our numerical model we apply maximum values to limit the shear and normal relative viscosities. The decision of having a maximum value is critical to the simulations likelihood of converging as the finite volume solver will struggle as the suspension starts to jam, this is similar to defining a yield viscosity for a Bingham fluid with the Generalized Newtonian Herschel-Bulkley model typically utilized in finite volume solvers[38]. Suspension Stress

Dhui = ∇ · Σp , Dt

(5)

M

hρi

AN US

From studying the theoretical developments for suspensions in the literature we discovered that many Brownian dynamics simulations and experiments had been performed for colloids and suspensions of larger particles [29]. Following the “Suspension balance model” approach of Nott and Brady (1984) [33] and the formulation of the stress tensor by Morris and Boulay[32] that reproduces the effects of the suspension balance model through the normal stress. The momentum balance of the mixture is

ED

where hρi is the average mixture mass density and hui is the average mixture velocity. The total stress tensor for the mixture is given by Σp = −µc ηn (φ)γQ ˙ + 2µc ηr (φ)D,

(6)

AC

CE

PT

where φ is the dispersed particle volume fraction, γ˙ is the shear strain-rate, ηn (φ) is the normal relative viscosity, D is the rate-of-deformation tensor of the mixture 1 D = (∇u + ∇uT ), (7) 2 and Q is the anisotropy tensor   λ1 0 0 Q =  0 λ2 0  , (8) 0 0 λ3

λ1 = 1.0, λ2 = 0.8 and λ3 = 0.5 are the anisotropy parameters. Note that this was derived for simple shear flow in the x−direction. An improvement to this approach would be to rotate the anisotropy tensor in the direction of the local flow, which has been previously studied for 2D flow by Miller et. al. (2009)[30]. The first normal stress difference is defined as N1 = Σ11 − Σ22 ,

(9)

and the second normal stress difference as N2 = Σ22 − Σ33 . 5

(10)

ACCEPTED MANUSCRIPT

CR IP T

Anisotropy of normal stresses implies that the first and second normal stress differences are non-zero. This anisotropy is required to cause the Weissenburg effect where the liquid suspension will dip down where a rod is rapidly rotating in the liquid. This dip would be due to a linear combination of N1 and N2 being negative as explained by Zarraga et. al. 2000[47] and is opposite to the rod climbing effect displayed by polymer solutions due to a positive N1 and N2 combination. Simulations in axisymmetric pipe flow (see below) revealed that the anisotropy in Q is subtle compared to an isotropic form of Q where it is set equal to the identity tensor Q = I. (11)

AN US

We find most axisymmetric simulations performed with both equation (8) and (11) give a similar solution to the migrated volume fraction profile (see Figure 3) because the term is dominated by the comparatively large ηn (φ), however for certain flows where normal stress differences are important we might expect some deviation from our results (for example a Weissenburg effect in rotating flow), this an aspect in need of further investigation. For the isotropic form the scalar particle pressure is simply Π = µc ηn (φ)γ˙ d .

(12)

3.0.3

M

A similar shear-induced particle pressure has been experimentally measured by Deboeuf et. al. [7] using a fine screen gap between concentric cylinders. Note that this particle pressure is applied in addition to the hydrodynamic pressure of the rest of the system. Euler-Euler Formulation

AC

CE

PT

ED

The models for suspension rheology presented so far are usually derived from a “single phase” perspective solving the momentum equation of the mixture, with the modeling of the volume fraction of particles as an extra transport equation. We decided to generalize the model into a fully multiphase model by following Lhuiller [26] who formulated the most complete two-phase model version for suspensions which emphasises how the particle stresses are linked with hydrodynamic stresses. The Euler-Euler methodology defines phases as inter-penetrating continua coexisting in the domain[19]. Each phase has its own distinct velocity, temperature and physical properties. Conservation equations are solved for each phase with additional closure laws to model the interactions between the phases. This approach is well established for bubbly flow in air and water, granular flow in air-solid mixtures, and oil and water for pipelines [27]. We solve the continuity and momentum equation for each phase, k. For the rest of the paper we use the suffix c for the continuous liquid phase and d for the dispersed particle phase. The conservation equations for mass and momentum take the following form ∂ (αk ρk ) + ∇.(αk ρk uk ) = 0 ∂t ∂ (αk ρk uk ) + ∇.(αk ρk uk uk ) = ∂t −αk ∇p + αk ρk g + ∇.τ k ± F k 6

(13)

(14)

ACCEPTED MANUSCRIPT

where F k are the interaction forces between the phases and are comprised of contributions from drag force, lift force, virtual mass force, and wall lubrication force. In this work we only considered the effect of drag force, F k = F d . For a Newtonian fluid the stress tensor is   1 τ N = 2αk µk D k − (∇.v k )I . (15) 3

Σshear = 2µc ηr (αd )D.

CR IP T

For suspensions and emulsions we have an additional stress tensor similar to the Morris and Boulay equation (6) which has separate equations for the dispersed particle and continuous fluid phases. The mixture shear stress is (16)

AN US

Note that traditionally the dispersed particle phase volume fraction is φ which in multiphase notation is now αd . In order to split this stress between the phases, we make the assumption that the averaged sum of the stress of the continuous phase and the stress of the dispersed phase is equal to the mixture stress, which is justified in that the relative viscosity accounting for the mixture appears equally in both terms (weighted by the respective volume fraction). Σp = αc τ c + αd τ d .

(17)

M

Comparing equation (17) with equation (6) we rearrange to find an expression for the continuous liquid (τ k = τ c ) that contains a shear stress contribution from the relative viscosity and using the Newtonian stress equation (15) we derive   1 τ c = 2αc µc ηr (αd )D c − (∇.v k )I . (18) 3

ED

The dispersed particle phase stress (τ k = τ d ) has a split contribution from shear along with the entire normal stress contribution (note that the Newtonian contribution for the particle phase is zero as µd = 0)

PT

τ d = 2αd µc ηr (αd )D d − µc ηn (αd )γ˙ d Q.

(19)

AC

CE

This completes the formulation of the stress tensors. The only term remaining in the momentum equations is the inter-phase drag. Whilst the flow for the suspensions considered here are in the Stokes regime with the particle Reynolds number, the ratio of inertial forces to viscous forces being small, Re << 1. The drag laws are implemented using a generic structure that covers a full range of flow velocities so is formulated in terms of Re. In our case the suspension rheology is restricted to Stokes or laminar flows so Re will never be large. We choose a drag force defined F D = AD |∆v r |.

(20)

where |∆v r | is the relative slip between phases |∆v r | ≡ |uc − ud |.

(21)

The direction of the drag force added to the momentum equation (14) is opposite to the drag of the other phase. 7

ACCEPTED MANUSCRIPT

The linearized drag coefficient, AD is given by AD =

a000 ρc |∆v r |CD . 8

(22)

where a000 is the symmetric interaction area density and is defined as 6αc αd d

(23)

CR IP T

a000 =

Rem =

AN US

and d is the particle diameter. The drag coefficient, CD is a the drag coefficient for particles. The Schiller-Naumann drag (for particles in a Newtonian fluid)[42] is given as 24 CD = (1 + Re0.681 ). (24) Re We take a novel approach of modifying the drag to account for the increase in viscosity of the mixture by including the relative viscosity so that the mixture single particle Reynolds number is ρc |∆v r |d Re = , µc ηr ηr

(25)

Finite Volume Method

CE

3.0.4

PT

ED

M

where Re is the Reynolds number for a particle in the Newtonian continuous liquid phase and d is the particle diameter. The modified drag for a suspension is now 24ηr CD = (1 + (Re/ηr )0.681 ). (26) Re This modification of drag has been shown to effectively slow the settling of the particles and is equivalent to the sedimentation function in Stokes flow. We limit the minimum value of Re to avoid the singularity in the function. In multiphase flow this has been traditionally handled with a correction to the drag to account for the presence of the other particles[40][41]. However in the context of rheological models it is the relative viscosity that accounts for the presence of the network of particles that increases the drag. This observation that the mixture viscosity replaces the solvent viscosity of the continuous suspending liquid has been observed previously by Lee et. al. [20] when considering the increase in Peclet number around the onset of shear-thickening of silica suspensions.

AC

The code is built using discrete equations using the finite volume method and the Phase-Coupled SIMPLE (PC-SIMPLE) algorithm[43, 45] which is in turn based on the well known SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm first proposed by Patankar and Spalding[36, 35]. All the phases share a common pressure field p (which is independent of the particle pressure discussed previously). The SIMPLE algorithm, that we iterate through to step n, reads as follows: 1. Set the boundary conditions. 2. Compute the gradients of phase velocities and guess the pressure, p∗ . 3. Solve the discretized momentum equations to compute the intermediate velocity fields, u∗ . 8

ACCEPTED MANUSCRIPT

4. Compute the uncorrected mass fluxes at faces (face mass flux is defined m ˙ f = ρuA, where A is face area). 5. Solve the pressure correction equation to produce cell values of the pressure correction, p0 .

7. Update the boundary pressure corrections p0b . 8. Correct the face mass fluxes: m ˙ n+1 =m ˙ ∗f + m ˙ 0f f

CR IP T

6. Update the pressure field: pn+1 = pn + urf · p0 (where urf is the underrelaxation factor for pressure).

9. Correct the cell phase velocities: We solve avP (un+1 − u∗ ) = −V ∇p0 for un+1 where V is the cell volume and avP is the diagonal entry for the linearized momentum equation coefficients (see Patankar[35] for full details).

Experiments

M

4

AN US

Diffusion terms are discretized for unstructured polyhedral meshes following Mathur and Murphy[28], as discussed in depth by Demirdˇzi´c[8]. Rhie and Chow[39, 11] interpolation is necessary for the pressure fields in the finite volume approach to avoid checker-boarding artefacts in pressure fields as discussed in greater depth by Tandon[44]. We found that particle pressure (equation (12)) was also in need of this numerical treatment.

AC

CE

PT

ED

A short overview of the experiments shall be followed by our detailed simulations of those experiments in the Results section. The Morris and Boulay formulation has been previously used within the context of a single phase mixture model applied with commercial finite element code by P. Mirbod[31] for the Couette experiment performed by Phillips et. al.[37], where Mirbod also extrapolated the model to simulate eccentric cylinders. We chose to model the Phillips Couette experiment to validate our code for simple shear in a rotating flow. We analyzed a detailed experimental data set of pressure driven pipe flow in several concentrations combined with nuclear magnetic resonance (NMR) data for the particle concentration by Hampton et. al. (1997)[14]. Finally the most complex geometry studied was the asymmetric forked bifurcation experiment of Xi and Shapely[46] which has been previously simulated by Yazaz Ahmed et. al.[2] as a single phase mixture using a finite volume method with a scalar volume fraction field and the diffusive flux model. To our knowledge this is the first time anyone has modeled these experiments with a Eulerian-Eulerian (multiphase) approach. Whilst this is computationally more expensive, as we have to solve a momentum equation for all phases, it allows slipping of the phases which occurs in many situations, for instance for the flow of a suspension of heavy particles around a bend, also we chose this formulation as it readily suits the modeling of surfactant stabilized emulsions by the mixing of two Newtonian liquids.

9

ACCEPTED MANUSCRIPT

0.65

0.65 0.6

0.6

0.55

0.55

0.5 0.5

φ

φ

0.45

0.45

0.4

0.35

0.3 0.25

0.3

0.2 0.4

0.6

0.8

1

0.4

r/R

CR IP T

0.4

φ = 0.55 φ = 0.50 φ = 0.45

0.35

n=50 100 200 800 12000

0.6

r/R

0.8

1

Figure 1: (a) Transient evolution of initial particle volume fraction, for n rotations of the inner cylinder, φ = 0.55 compared to experiment. (b) Steady state particle volume fractions for three initial volume fractions.

Results

AN US

5

5.1

Couette Flow

M

Most simulations used a maximum packing of the suspension of φm = 0.68. This value is higher than the theoretical monodisperse value of random close packing which is around 0.645 in order to reflect that some polydispersity that exists in the particles used in the experiments, it is also the same value used in previous studies[32].

AC

CE

PT

ED

We simulated a circular 2D cross section of the Couette device with the same dimensions as the original Phillips et. al. experiment[37] consisting of a cell of inner radius 0.64cm and outer radius as 2.38cm. We performed the simulation on a coarse mesh of 1600 cells and fine mesh of 21600 cells to ensure no dependence on the mesh and found we got the same results, we present the results from the finer mesh. Three simulations were performed at initial volume fractions of φ = 0.45, φ = 0.50 and φ = 0.55 using a time step of dt = 0.05s. The material data was matched to the experiment where Poly(methylmethacrylate) (PMMA) spherical particles of density 1182kgm−3 were suspended in a Newtonian tetrabromoethane/oil mixture with the same density as the particles and a dynamic viscosity of µc = 4.95P a.s. We found the steady state results for this geometry were largely independent of Kn , so we present the results using the value of Kn = 0.75 from the original formulation[32]. The simulations were initially run with maximum packing volume fraction of φm = 0.68, however near the wall at the highest volume fraction near steady state when the particle phase tries to exceed a volume fraction of 0.6 we found that the normal stress repels the particles away effectively blocking further migration. Increasing φm to 0.72 allowed a better result. This increase in φm is not unphysical as several experimental studies have shown that as the particles experience a strong shear flow then they align into rows exhibiting a hexagonal close packing[34] that can increase φm in the rheology[16]. The Phillips et. al. paper does not clearly state the angular frequency for each run, only that it varied between 17 and 117rpm. We ran all of the

10

ACCEPTED MANUSCRIPT

5.2

AN US

CR IP T

simulations at 17rpm, as higher speeds resulted in an over-prediction in volume fraction of particles above the experimental data points, which is something we need to investigate further in later work. We ran three simulations at initial volume fractions of φ = 0.45, φ = 0.50 and φ = 0.55. Fig.1(a) shows the transient evolution of the particle concentration. When the simulations starts a large wave of particles pushes away from the wall and travels towards the outer edge. For the case of φ = 0.55 transient data is given in terms of number of revolutions of the inner wall. This case was simulated at a frequency of f = 27rpm = 27/60Hz. Therefore the period of rotation is T = 60/27 = 2.222s. 50 revolutions corresponds to a time of 111.11s. Fig.1(b) shows the steady state values for all three initial volume fractions studied. We see a good agreement with the experimental data with no fitting or adjustment of the Morris and Boulay model parameters. It is not clear from the experimental data what happens to the particles at the static wall, as the last data point stops somewhat short, however most of the simulations predict a smooth continuation in volume fraction up to the wall. In general the results of the simulations give a good agreement with the data, we find a slight under-prediction of particle concentration nearer the inner cylinder, and a slight over-prediction at the outer edge. It may be possible to lower Kn to get a more accurate fit here.

Pressure Driven Pipe Flow

AC

CE

PT

ED

M

Hampton et. al.[14] performed a series of experiments coupled with NMR data of a suspension of spherical PMMA particles flowing in a pipe where the particles are neutrally buoyant with the suspending oil. The experiments where performed with particles of sizes 650µm in a narrow pipe of internal diameter 25.4mm which has a ratio of diameter of particle to the diameter of the pipe, a/R = 0.265 (here labeled pipe-1) and a wider pipe of diameter 50.8mm with larger particles on average 3175µm which gives a ratio a/R = 0.0625 (labeled pipe-2). The Newtonian oil has a density 1180.7kgm−3 and molecular viscosity 2.1P a.s. The simulations were run using the Morris and Boulay model outlined above until steady state was achieved. The cross section of the volume fractions were recorded at a length of 120 diameters of the pipes downstream, giving lengths L = 3.048m for the pipe-1 and L = 6.096m for pipe-2. Each pipe had data recorded for suspensions of initial particle volume fractions of 0.2, 0.3 and 0.45. The simulations were transient with time-steps between dt = 0.001 − 0.01s and were run until steady state was achieved. The suspension was modeled in an axisymmetric pipe of radius 1.27cm and length 3.048m with 15000 cells. The fine mesh was necessary to correctly capture the particle migration along the radial direction. We limited the maximum value of the normal relative viscosity ηnmax = 2000 as the simulation could diverge if the value becomes too high (due to a high volume fraction). All simulations were initially run with default values from Morris and Boulay (1999) of φm = 0.68, Kn = 0.75. Some of the results did not give enough migration of the particles towards the centre of the pipes. In these cases the value of the normal stress contact parameter Kn was increased. The Hampton paper did not specify the real velocity of the experiments so we chose a default input plug flow velocity of 0.05ms−1 unless otherwise stated, we also normalize the results by this inlet velocity. The results were generally independent of velocity, however some cases reached 11

ACCEPTED MANUSCRIPT

0.6

2 Particle volume fraction Experimental data

0.5

1.5

0.3

1

u

φ

0.4

0.2

CR IP T

0.5 0.1

Normalized particle velocity Experimental data

0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

r/R

0.4

0.6

0.8

1

r/R

Figure 2: Pipe-2 with initial particle volume fraction of φ = 0.2. Outlet steady state axisymmetric (a) particle volume fraction, (b) Normalized particle velocity profiles. 2

AN US

0.7

Particle volume fraction Kn = 0.75 Particle volume fraction Kn =2 Particle volume fraction Kn =2, anisotropic Experimental data

0.6 0.5

φ

u

0.4

1.5

0.3 0.2

1

0.5

0.1 0 0

0.2

0.4

0.6 r/R

M

Normalized particle velocity Experimental data

0.8

1

0

0

0.2

0.4

0.6

0.8

1

r/R

ED

Figure 3: Pipe-2 with initial particle volume fraction of φ = 0.3. Outlet steady state axisymmetric (a) particle volume fraction, (b) Normalized particle velocity profiles.

AC

CE

PT

steady state earlier with higher speeds. The results for the Pipe-2 are fairly consistent with all runs using the default values from the Morris and Boulay paper. Fig. 2 shows the result for initial particle volume fraction of φ = 0.2 where we do not predict such a large depletion of the particles away from the wall. In the experimental data φ approaches zero at the wall r/R = 1. Clearly we need additional physics to account for this effect, the explanation of which is likely to be attributed to an inertial lift force as observed by Han et. al.[15, 23]. There is an over-prediction of the center line particle volume fraction with all the simulations in general. The velocity profiles demonstrate a remarkable agreement with the experimental data. Fig. 3 shows the result for initial particle volume fraction of φ = 0.3, we indicate how increasing the value of Kn results is a broadening of the distribution of particles away from the pipe axis, along with a reduction in the peak particle volume fraction, note this also raises the profile of the volume fraction above the data elsewhere across the radius. Another way to reduce the peak of the volume fraction at the center line is to reduce the maximum packing volume fraction, although it is hard to justify why you would have to do this with

12

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

M

Figure 4: Pipe-2 3D directed mesh of 1,250,000 cells.

0.7

2

PT

0.6

φ

0.4

CE

0.3

1.5

u

0.5

0.2

0.5

Axisymmetric particle volume fraction 3D particle volume fraction Experimental data

0.1

Axisymmetric normalized velocity 3D normalized velocity Experimental data

0

AC

0

0.2

1

0 0.4

0.6

0.8

1

0

r/R

0.2

0.4

0.6 r/R

0.8

1

Figure 5: Pipe-2 with initial particle volume fraction of φ = 0.45. Outlet steady state axisymmetric and 3D (a) particle volume fraction, (b) Normalized particle velocity profiles.

13

ACCEPTED MANUSCRIPT

0.6

2 Particle volume fraction Kn = 3.5 Particle volume fraction Kn = 0.75 Experimental data

0.5

1.5

0.3

u

φ

0.4

1

0.2

CR IP T

0.5 0.1

Normalized particle velocity Experimental data

0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

r/R

0.4

0.6

0.8

1

r/R

Figure 6: Pipe-1 with initial value of φ = 0.2. Outlet steady state (a) particle volume fraction, (b) Normalized particle velocity profile. 2

0.6

AN US

0.7

1.5

0.5

u

φ

0.4 0.3 0.2

1

0.5

Particle volume fraction Kn = 3.5 Particle volume fraction Kn = 0.75 Experimental data

0.1

Normalized particle velocity Experimental data

0

0

0.2

0.4

0.6 r/R

0.8

1

M

0

0

0.2

0.4

0.6

0.8

1

r/R

ED

Figure 7: Pipe-1 with initial value of φ = 0.3. Outlet steady state (a) particle volume fraction, (b) Normalized particle velocity profile.

AC

CE

PT

the spherical particles studied. We also compare a run using the anisotropic form of the anisotropy tensor Q where the central peak is slightly higher and is equivalent to a run with a slightly smaller value of Kn , nevertheless the effect is subtle for pipe flow. For the initial particle volume fraction of φ = 0.45 we ran an additional 3D simulation to compare with the axisymmetric case. The 3D mesh consisted of 1, 250, 000 cells (Fig. 4). The comparison (Fig. 5) shows a very good agreement and acts to validate the axisymmetric result against the 3D case (which is more accurate due to the finer mesh). The simulation using pipe-1 required greater adjustment in normal contact Kn to get an agreement with experimental data. Fig. 6 shows the simulation of the case with initial particle volume fraction of φ = 0.2 with initial plug flow velocity of 0.5ms−1 . To achieve a match with the experimental volume fraction profile at the end of the pipe it was necessary to increase the value of the normal contact contribution Kn = 0.75 to Kn = 3.5, otherwise the particles migrated into a sharp peak at the centreline. Increasing Kn might compensate for the lack of wall depletion of particles with our model or it might be that our profile section is not being taken at a great enough distance down the pipe, as the inner build up of particles increases downstream with distance. Fig. 7 shows the simulation of the case with initial particle volume fraction 14

ACCEPTED MANUSCRIPT

0.7

2

0.6 1.5

0.5

u

φ

0.4 1

0.3 0.2

Normalized particle velocity Experimental data

0

0 0

0.2

0.4

0.6

CR IP T

0.5 Particle volume fraction Kn = 3.5 Experimental data

0.1

0.8

1

r/R

0

0.2

0.4

0.6

0.8

1

r/R

Figure 8: Pipe-1 with initial value of φ = 0.45. Outlet steady state (a) particle volume fraction, (b) Normalized particle velocity profile.

AC

CE

PT

ED

M

AN US

of φ = 0.3 with initial plug flow velocity of 0.3ms−1 . Again we find an underprediction of the width of the peak volume of fraction of particles near the centreline of the pipe. We compensate for the lack of migration of particles from the wall by increasing Kn to 3.5, although this causes the volume fraction to drop below the experimental value on the centreline which would be remedied by the availability of more particles being pushed away from the wall. Fig. 8 shows the simulation with initial particle volume fraction of φ = 0.45 with initial plug flow velocity of 0.2ms−1 . We find relatively good agreement with the data, the volume fraction is slightly under-predicted between r/R = 0.2 to 0.4 but if we correctly captured the wall effects the extra particles would have brought this line up to the data. In start-up of flow there is a push of particles away from the wall which drives the migration of particles towards the centre of the pipe. The normal stress increases linearly with the strain-rate, however from observation of the data it is clear this force alone is not enough to account for the presence of the wall. Perhaps the wall force varies quadratically with strain rate as observed in some viscoelastic fluids? Unfortunately adding an additional inertial lift force is beyond the scope of this paper and something we would like to address in future work. Regardless, it is impressive that the pipe data matches the predictions of the Morris and Boulay constitutive model with very minor fitting. The necessary increase in Kn observed may imply that the value of the normal relative viscosity is too small at low volume fractions, or that there is variation with the size or roughness of the particles or it could be that the plane section has not been taken far enough downstream from the pipe inlet to allow the full profile to develop. It is reassuring that the adjustment to Kn was consistent for the whole set of wide pipe/large particle (pipe-2) or thin pipe/small particle (pipe-1) sets of data at all three initial volume fractions.

5.3

Asymmetric Flow Bifurcation in a Channel

The asymmetric forked bifurcation experiment of Xi and Shapely[46] has been previously simulated as a single phase mixture with the finite element method by Yazaz Ahmed et. al.[2] where they programmed the original diffusive flux model of Phillips[37] considering the volume fraction of particles as a scalar

15

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 9: Velocity field of the mixture showing the center-plane of the geometry with “inlet”, “main” and “side” planes shown along with planes on the outlets.

AC

CE

PT

ED

M

quantity. Whilst this method was fairly successful at predicting the velocity profiles, the particle volume fractions did not show the accuracy that we have presented below, as we improve upon those results with the implementation of the Morris and Boulay constitutive model within a full Euler-Euler multiphase framework. The geometry consists of 30mm long sections of rectangular channel of 3mm width and 11.5mm height with a gap of 4mm between the two branches. We simulated the experimental geometry in both in 2D and 3D. The 2D mesh consisted of 13,320 square cells and the 3D mesh consisted of 2.639 million regular hexahedral cells. The inlet boundary conditions for velocity and volume fraction area were scanned from the 2D experimental data and interpreted as a data map using MATLAB, this saved us from modeling the larger upstream pipe which would have resulted in a much larger simulation. Note that the pipe upstream of the branched flow section is offset resulting in an unusual non-uniform particle volume fraction distribution and velocity at the inlet. The initial volume fraction was set to the average value of 0.4 throughout the domain, whilst the fluid and particles are initially at rest. The inputs parameters were matched to the experiment, with the liquid composed as a glycerin in water solution matching the density of the particles, ρ = 1180kgm−3 with a Newtonian viscosity of µc = 0.023P a.s. The interaction length scale was set to 85µm which was the average diameter of the PMMA particles. The initial average particle volume fraction is 0.4. Two simulations were performed, one run as steady state with Kn = 0.75 run for 11000 iterations and one transient with Kn = 2 and a time-step of dt = 0.1s. The results between the two were almost identical in terms of velocity and particle volume fraction profiles. We present the results for Kn = 0.75 as this is consistent with the other experiments studied. Fig. 9 shows the velocity distribution of the mixture scanned at a position 7mm upstream of the branch section which is the inlet to our mesh, a corresponding particle volume fraction was also scanned and used as the boundary condition. 7mm downstream of the branch section we take a plane and compare our velocity distribution to the experimental field. Fig. 10 shows that the predicted flow is in good agreement with the experiment showing the fluid being slowed where the particles are present due to the increase in relative viscosity. The experiment shows the flow is slightly weaker in the side branch but we do

16

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 10: Axial velocity map (cms−1 ) downstream of the fork (y = 7mm) showing the experimental NMR scan on the left[46] and our Computational Fluid Dynamics (CFD) simulation on the right, lines marking lateral and spanwise directions are shown on the main and side planes. 3D CCM+ inlet A1 3D CCM+ main A2 3D CCM+ side A3 Experiment inlet A1 Experiment main A2 Experiment side A3

8

6 5 4 3 2 1 0.1

0.2

M

Velocity (cm/s)

7

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

ED

-z/D+1

PT

Figure 11: Velocity profiles on the center lines A1 inlet branch (y = −7mm), A2, main branch (y = 7mm) and A3, side branch(y = 7mm). Lines are from the 3D simulation and points are extracted from Xi and Shapely’s experimental data.

AC

CE

not see this as clearly in the simulation. Fig. 11 illustrates the velocity of the mixture on the line probes defined in Fig. 10. The A1 line has been mapped from the 2D inlet boundary data and agrees with the data points from the original paper. The curves at A2 and A3 are the velocities spanning from the bottom to the top of the channels in the main and side branches predicted by simulation where we can see that we slightly underestimate the peak and trough in velocity around the particles. Note that there is little slip between the phases with the velocity of the liquid equal to the velocity of the particles in all the data presented (although we plot the particle velocity). Fig. 12 compares the mixture velocity on a plane across the channel, with B2 across the main branch and B3 across the side branch. Here we show the 2D simulation on the left and the 3D on the right we can see we follow the same trend in asymmetry as the experiment. The inflow profile of particle volume fraction demonstrates elliptically shaped

17

ACCEPTED MANUSCRIPT

Velocity (cm/s)

6 5 4 3 2

6

1 0

3D CCM+ Inlet B1 3D CCM+ Main B2 3D CCM+ Side B3 Inlet branch data B1 Main branch data B2 Side branch data B3

7

Velocity (cm/s)

Inlet branch data B1 Main branch data B2 Side branch data B3 2D CCM+ Inlet B1 2D CCM+ main B2 2D CCM+ side B3

7

5 4 3 2 1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

1-x/2B

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1-x/2B

M

AN US

CR IP T

Figure 12: Velocity profiles across the width of the channels along the lines B1 at the inlet, B2 in the main branch and B3 on the side branch for the (a) 2D simulation, (b) 3D simulation.

AC

CE

PT

ED

Figure 13: particle concentration showing the center-plane of the geometry with “inlet”, “main”, “side” and the two outlet planes shown.

Figure 14: Axial particle concentration map downstream of the fork (y = 7mm) showing the experimental NMR scan on the left[46] and our CFD simulation on the right, lines marking lateral and span-wise directions are shown on the main and side planes.

18

ACCEPTED MANUSCRIPT

Inlet branch A1

Main branch A2

0.4 0.3 0.2 Simulation Experimental data

0.1 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Side branch A3

0.5

0.5

0.45

0.45

Particle volume fraction

0.5

Particle volume fraction

Particle volume fraction

0.6

0.4 0.35 0.3 0.25 0.2 0.15 0.1

Simulation Experimental Data

0.05 0

1

0

0.1

0.2

0.3

0.4

-z/D+1

0.5

0.6

0.7

0.8

0.9

0.4 0.35 0.3 0.25 0.2 0.15 0.1

0

1

Simulation Experimental data

0.05 0

0.1

0.2

0.3

-z/D+1

0.4

0.5

0.6

0.7

0.8

0.9

1

-z/D+1

Main branch B2

Side branch B3

0.6

0.6

0.5

0.5

0.5

0.4 0.3 0.2 3D CCM+ Inlet Inlet branch data

0.1 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Particle volume fraction

0.6

Particle volume fraction

Particle volume fraction

Inlet branch B1

0.4 0.3 0.2 3D CCM+ Main Main branch data

0.1 0

0

1-x/2B

0.1

0.2

0.3

0.4

0.5

1-x/2B

0.6

CR IP T

Figure 15: Particle concentration profiles on the center lines A1 inlet branch (y = −7mm), A2 main branch (y = 7mm) and A3 side branch(y = 7mm). Lines are from the 3D simulation and points are extracted from Xi and Shapely’s experimental data.

0.7

0.8

0.9

1

0.4 0.3 0.2

3D CCM+ Side Side branch data

0.1 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1-x/2B

AN US

Figure 16: Particle concentration profiles across the width of the channels along the lines (a) B1 at the inlet, (b) B2 in the main branch and (c) B3 on the side branch for the 3D simulation.

6

PT

ED

M

iso-concentration lines which are then split at the junction giving two “half ellipses” in the two branches at the plane 7mm downstream from the split as seen in Fig. 13. This shows agreement with experiment (Fig. 14) although the experiment has some scatter in concentration which is smooth in the simulation. Interestingly the simulation predicts some accumulation of particles in the corners, which can be seen in some parts of the experiment, this is evident more at the outlets in Fig. 13 where the concentration tends back to an more elliptical shape. Studying the plots across the lines A1, A2 and A3, we see in Figs. 15 and B1, B2 and B3 16 that we get a excellent agreement with experiment of the shape of the particle distributions.

Conclusion

AC

CE

We have for the first time implemented the Morris and Boulay constitutive model for monodisperse non-colloidal suspension rheology into an Euler-Euler multiphase finite volume solver. By solving the momentum equations for each phase and splitting the stress between them we then solve the particle concentration diffusion. We have modified the particle drag by introducing the relative viscosity into the standard Schiller-Naumann drag for a particle in a Newtonian liquid. The resulting set of equations has been applied to all of the high quality NMR experimental data that we were able to find in the academic literature resulting in predictions that nearly all the volume fraction and velocity profiles agree accurately with the assumption that the normal stress of the original model is isotropic and acts like a pressure. In future work we could revisit the anisotropy tensor Q, one approach being to rotate the local fluid element is in the direction of flow in each cell by extrapolating the method of Miller et. al.[30]. However we were surprised that 19

ACCEPTED MANUSCRIPT

7

CR IP T

the assumption of an isotropic particle pressure gave such good results in comparison with the experiments studied and suggests that the effect of anisotropy may be subtle, or perhaps more relevant in other flow geometries. Clearly for more accurate validation of these models, it is the turn of the experimentalists to perform accurate NMR simulations in more complex geometries such as contraction-expansion flow to really test for anisotropic normal stress effects. The models utilized so far have a relative viscosity that is only dependent on particle volume fraction, but future models may involve shear-rate or even the full deformation rate tensor; although a preliminary exploration showed that shear thinning or thickening effects were not necessary to explain any of the physics in the experiments studied in this work.

Acknowledgements

AN US

We would like to thank Andrew Splawski and Dr. Alexander Vichansky for useful discussions regarding the formulation of the numerics and code design, also the patience of my colleagues in allowing me to develop this model to completion. Additionally I would like to thank Paul Dawson and Dr. Samir Muzaferija for their help on the finite volume method.

References

M

[1] Y. Ahmad, M. G. Kumar, and A. Singh. CFD simulation of shear induced particle migration. ICMF2010, 30 May-4 June 2010, Tampa, Florida, USA., 2010.

ED

[2] G.M. Yezaz Ahmed and A. Singh. Numerical simulation of particle migration in asymmetric bifurcation channel. Journal of Non-Newtonian Fluid Mechanics, 166(1-2):42–51, 2011.

PT

[3] G. K. Batchelor and J. T. Green. The determination of the bulk stress in a suspension of spherical particles to order c2 . J. Fluid Mech., 56(3):401–427, 1972.

CE

[4] G. K. Batchelor and J. T. Green. The hydrodynamic interaction of two small freely-moving spheres in a linear flow field. J. Fluid Mech., 56(2):375– 400, 1972.

AC

[5] F. Boyer, E. Guazzelli, and O. Pouliquen. Unifying suspension and granular rheology. Physical Review Letters, 107:188301, 2011.

[6] M. E. Cates, J. P. Wittmer, J. P. Bouchaud, and P. Claudin. Jamming, force chains, and fragile matter. Physical Review Letters, 81(9):1841–1844., 1998.

[7] A. Deboeuf, G. Gauthier, J. Martin, Y. Yurkovetsky, and J. Morris. Particle pressure in a sheared suspension: A bridge from osmosis to granular dilatancy. Phys. Rev. Lett., 102:108301, 2009. [8] I. Demirdˇzi´c. On the discretization of the diffusion term in finite-volume continuum mechanics. Numerical Heat Transfer, Part B: Fundamentals: 20

ACCEPTED MANUSCRIPT

An International Journal of Computation and Methodology, 68(1):1–10, 2015. [9] E. C. Eckstein, D. G. Bailey, and A. H. Shapiro. Self-diffusion of particles in shear flow of a suspension. J. Fluid Mech., 79(1):191–208, 1977.

CR IP T

[10] A. Einstein. Eine neue bestimmung der molek¨ uldimensionen. Annalen der Physik, 19(2):289, 1906. [11] Joel H. Ferziger and Milovan Peri´c. Computational Methods for Fluid Dynamics. Springer, 1996. [12] F. Gadala-Maria and A. Acrivos. Shear-induced structure in a concentrated suspension of solid spheres. J. Rheol., 24(6):799–814, 1980.

[13] Thomas Graham. Liquid diffusion applied to analysis. Phil. Trans. R. Soc. Lond., 151:183–224, 1861.

AN US

[14] R. E. Hampton, A. A. Mammoli, A. L. Graham, N. Tetlow, and S. A. Altobelli. Migration of particles undergoing pressure-driven flow in a circular conduit. Journal of Rheology, 41:621, 1997.

[15] M. Han, C. Kim, M. Kim, and S. Lee. Particle migration in tube flow of suspensions. Journal of Rheology, 43:1157–1174, 1999.

M

[16] L. Heymann, S. Peukert, and Nuri Aksel. On the solid-liquid transition of concentrated suspensions in transient shear flow. Rheol Acta, 41:307–315, 2002. [17] E. J. Hinch. A perspective of batchelor’s research in micro-hydrodynamics. J. Fluid Mech., 663:8–17, 2010.

PT

ED

[18] N. J. Inkson, J. Plasencia, and S. Lo. Predicting emulsion pressure drop in pipes through CFD multiphase rheology models. Proceedings of the 10th International Conference on CFD in Oil & Gas, Metallurgical and Process Industries SINTEF, Trondheim, Norway, pages 461–466, 2014. [19] M. Ishii and T. Hibiki. Springer, 2006.

Thermo-Fluid Dynamics of Two-Phase Flow.

CE

[20] J. H. So J. D. Lee and S. M. Yang. Rheological behavior and stability of concentrated silica suspensions. Journal of Rheology, 43:1117, 1999.

AC

[21] A. Karnis, A. L. Goldsmith, and S. G. Mason. The kinetics of flowing dispersions: concentrated suspensions of rigid particles. J.Colloid Interface Sci., 22:531–553, 1966. [22] I. M. Krieger and T. J. Dougherty. A mechanism for non-Newtonian flow in suspensions of rigid spheres. Transactions of the Society of Rheology, 3:137–152, 1959.

[23] Pandurang M. Kulkarni. Suspension mechanics at finite inertia. Technical report, PhD Thesis, City University of New York, 2009. [24] P. Langevin. Sur la th´eorie du mouvement brownien. C. R. Acad. Sci., 146:530–533, 1908. 21

ACCEPTED MANUSCRIPT

[25] D. Leighton and A. Acrivos. The shear-induced migration of particles in concentrated suspensions. J. Fluid Mech., 181:415–439, 1987. [26] D. Lhuillier. Migration of rigid particles in non-Brownian viscous suspensions. Phys. Fluids, 21(2):023302, 2009.

CR IP T

[27] S. Lo. Industrial Computational Fluid Dynamics: September 21-25, 2015. Von Karman Institute for Fluid Dynamics, 2015. [28] S. R. Mathur and J. Y. Murthy. A pressure-based method for unstructured meshes. Numerical Heat Transfer, Part B: Fundamentals: An International Journal of Computation and Methodology, 31:195–215, 1997. [29] J. Mewis and N. J. Wagner. Colloidal Suspension Rheology. Cambridge University Press, 2012.

AN US

[30] R. M. Miller, J. P. Singh, and J. F. Morris. Suspension flow modeling for general geometries. Chemical Engineering Science, 64:4597–4610, 2009.

[31] P. Mirbod. Two-dimensional computational fluid dynamical investigation of particle migration in rotating eccentric cylinders using suspension balance model. International Journal of Multiphase Flow, 80:79–88, 2016. [32] J. F. Morris and F Boulay. Curvilinear flows of noncolloidal suspensions: The role of normal stresses. J. Rheol., 43:1213, 1999.

M

[33] P. R. Nott and J. F. Brady. Pressure-driven flow of suspensions: simulation and theory. J. Fluid Mech., 275:157–199, 1994.

ED

[34] T. Palberg and R. Biehla. Sheared colloidal crystals in confined geometry: a real space study on stationary structures under shear. Faraday Discuss., 123:133–143, 2003. [35] S. V. Patankar. Numerical Heat Transfer and Fluid Flow. Taylor & Francis, 1990.

PT

[36] S. V. Patankar and D. B. Spalding. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. J. of Heat and Mass Transfer, 15(10):1787–1806, 1972.

CE

[37] R. J. Phillips, R. C. Armstrong, R. A. Brown, A. L. Graham, and J. R. Abbott. A constitutive equation for concentrated suspensions that accounts for shear-particle migration. Phys. Fluids. A, 4:30, 1991.

AC

[38] R. J. Poole and R. P. Chhabra. [39] C. M. Rhie and W. L. Chow. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA Journal, 21(11):1525–1535, 1983.

[40] J. Richardson and W. Zaki. Sedimentation and fluidisation: part i. Trans. Inst. Chemical Engineers, 32:35–53, 1954. [41] H. Rusche and R. Issa. The effect of voidage on the drag force on particles, droplets and bubbles in dispersed two-phase flow. Japanese European Twophase Flow Meeting, Tsukuba, Japan, 2000. 22

ACCEPTED MANUSCRIPT

[42] L. Schiller and A.Z. Naumann. Ver.Deut.Ing,, 77:318–320, 1933. [43] Siemens PLM Software. STAR-CCM+ documentation. 2017. [44] M. Tandon. Structured adaptive mesh refinement (SAMR) simulation study of the buoyant plumes. Technical report, Ph.D. Thesis, Department of Chemical Engineering, University of Utah, 2008.

CR IP T

[45] S. A. Vasquez and V. A. Ivanov. A phase coupled method for solving multiphase problems on unstructured meshes. Proc. of ASME FEDSM, pages 734–748, 2000. [46] C. Xi and N. C. Shapley. Flows of concentrated suspensions through an asymmetric bifurcation. J. Rheol., 52:625, 2008.

AC

CE

PT

ED

M

AN US

[47] I. E. Zarraga, D. A. Hill, and David T. Leighton. The characterization of the total stress of concentrated suspensions of noncolloidal spheres in Newtonian fluids. J. Rheol., 44:185, 2000.

23