Flows of Dense Suspensions of Polymer Particles through Oblique Bifurcating Channels: Two Continua Approach

A two-velocity mathematical model is proposed for dense suspension flows through channel bifurcations. Equations agree with thermodynamic laws and they are suitable for both heavy and light particles. The pulsatile mode of injection of particles is considered. In the 2D-case, we address the issue of partitioning particles and study how a loss of particles into the side branch depends on the bifurcation angle. A qualitative agreement with experiment data are established. We capture the Zweifach–Fung effect. We treat polymer particles as a phase enjoying the rheology of the Bingham viscoplastic material. We prove that the polymer particle distribution between two branches correlates with the averaged-in-time Bingham number in these branches.


Introduction
Suspension flows through bifurcating channels are relevant to many natural and technological processes. Their application concerns medicine, biotechnology and microfluidic devices. Micromoulding of thermoplastic polymer is a developing process with great potential for producing low-cost microfluidic devices. Among different micromoulding techniques, micro-injection moulding is one of the most promising processes suitable for manufacturing polymeric disposable microfluidic tubes with branches. The review paper [1] presents the main significant developments that have been achieved in different aspects of micro-injection moulding of microfluidic devices. Aspects covered include device design, machine capabilities, mould manufacturing, material selection and process parameters. Problems, challenges and potential areas for research are highlighted. One more application of the present paper concerns the control of proppant flow through perforations in the wellbore. Such a problem is important in hydraulic fracturing technology [2].
We develop a thermodynamically consistent two-velocity model of a granular fluid, taking into account the difference in the rheology of the carrier fluid and the solid phase consisting of particles. There is no assumption on dilute suspensions, and the model allows one to describe such processes as the deposition of particles. The validation of the model was carried out on a benchmark experiment known as the Boycott effect [3], implying that an inclination of a vertical vessel enhances the sedimentation of particles in suspensions. Historically, suspensions have been considered within the framework of a single-phase continuum with a viscosity dependent on particle concentration. Currently, two-phase approaches are being used more and more. They are usually based on methods of averaging either over a volume or over an ensemble of particles. In this paper, when deriving equations, we follow a different method, which was originally developed by Khalatnikov and Landau for the thermodynamics of superfluid helium 2He [4]. According to this method, the forces of interfacial interaction for reversible processes are first uniquely determined, since the energy conservation makes the entire system of conservation laws overdetermined. Then, the dissipative forces of phase interaction are determined by matching the conservation laws with the general principles of irreversible de Groot-Mazur thermodynamics [5].
In a number of works on two-phase flows, additional forces are introduced into the momentum equations, such as the force of Archimedes, Stokes, Suffman, Magnus. Even forces that depend explicitly on time are used (the Basset-Boussinesq force) [6]. In this work, instead of additional forces, we use the generalized Fick diffusion law for the particle mass concentration flux vector. This approach is consistent with thermodynamics and proved to be successful in sedimentation problems if the gravitational component is taken into account in Fick's law. The generalized Fick's law takes into account not only the concentration gradient, but also the gradients of temperature, pressure, and the modulus of the phase velocity difference. The Stokes resistance force is taken into account through the coefficient of interfacial friction, which satisfies the known correlations. The diffusion approach we use imposes limitations on the applicability area. The particles must be sufficiently small.
In the present paper, the problem of the flow when there is a branch in the main channel is considered. The pulsatile regime, when the injection of particles alternates with the injection of a clean fluid, is studied in detail. In a number of papers [7], suspensions of polymer particles are treated as non-Newtonian fluids. Here, we follow the same approach and the particulate phase is represented by the rheology of a viscoplastic Bingham fluid. We establish that the loss of particles into the branching channel is considered as the smaller, the greater the angle of branching from the direction of the main channel. The dependence of the flow on the viscosities and the yield stress of the granular phase is addressed.
Recently, interest in the flow of suspensions through channel bifurcations has been growing in microfluidics and biology due to the problem of particle separation [8]. Partitioning of erythrocytes, leukocytes, and platelets at vascular junctions plays an important role in determining microvascular hematocrits and has important physiological implications. One of the advantages of continuum models is that solvers known in computational fluid mechanics can be applied to them. In addition, these solvers can be easily adapted to complex geometries and non-Newtonian fluid rheology.
A number of approaches have been developed, which use external forces to affect the fluid/particle behavior in microchannel-based separation systems. Such forces can be magnetic [9], optical [10], acoustic [11], chemical [12], thermal [13], gravitational [14,15], electrical and electronic [16,17]. The most popular separation method is based on gravitation. However, this technique is not as effective in microchannel devices because of the low Reynolds number, laminar flow and interface tension [18]. In macro-scale separation systems, gravitational approaches are based on the forces of inertia to help speed up the separation process, an example of this being centrifugal sorting [19]. This has led to an interest in adopting new microscopic separation methods, in which viscous forces, shear strain rate, interface tension and the geometrical effects-often used for inertial focusing-play an important role. Various devices have been designed that make use of the effects of hydrodynamics and viscous forces. Examples include microfilter devices [20,21] and microchannel devices [22,23]. Furthermore, the bending channel [24,25], channel constriction [26,27], and bifurcated channels [28,29] have been explored. In the bending channel design, inertial effects are exploited whereby a centrifugal force is produced on the particles when passing through the curved section of the channel and the skimmed liquid is likely to be formed in the central region. Flow velocity is a key factor for the effectiveness of this design. A channel constriction is used to force particles to move to the channel center.
Fluid separation in T-shaped microchannels is studied in many papers, see, for example [30]. Xi and Shapley studied the flow of concentrated suspension in an asymmetric T-junction bifurcation of rectangular channels with nuclear magnetic resonance imaging [31]. They observed that the particles are almost equally partitioned between the downstream branches, and this indicates the migration of particles across the dividing streamlines near the bifurcation section.
In the present paper, we present a modeling analysis for dense suspension flows through oblique bifurcating channels. Such an issue has been addressed in [32] for steady flows and in the case of neutrally buoyant particles, and in the low Reynolds number regime. We focus on the difference between fluid and particle velocities, interphase friction, viscous forces, non-Newtonian rheology and geometrical effects. Particularly, we consider in detail the pulsatile mode of injection of particles and study how the loss of particles into the side branch depends on the bifurcation angle. We prove that the partitioning of particles occurs in agreement with the Zweifach-Fung effect, stating that particles prefer the highflow-rate branch. Our method is suitable for both heavy and light particles. Qualitative agreements with experiment data are established. Paying emphasis on pulsatile injection of polymer particles, we establish that more particles fall into the branch with a lower mean Bingham number. A numerical algorithm is developed within the framework of the FreeFem++ package. In the incompressible case, a modified SIMPLE method is applied in order to ensure the divergence-free condition for the average volumetric velocity. The stability test was carried out by mesh refinement. The approach has been validated in the previous authors' paper concerning the Boycott sedimentation effect [33].

Mathematical Model
We consider a joint flow of two continua when an arbitrary volume V contains a fluid (index f ) and a granular phase (index s). Volume, mass and pressure of the fluid and the granular phases is denoted by V f , m f , p f and V s , m s , p s , respectively. It is assumed that the granular phase is a mixture of dry particles and a carrier fluid, such as proppant gel. In this case, V s = V M + V p and m s = m M + m p where V M is the mud volume, V p is the volume of dry particles, m M is the mud mass, and m p is the mass of dry particles.
The particles are "frozen" in the carrier liquid, i.e., the granular phase is characterized with just one speed v s , one viscosity and one stress tensor. In what follows, the indexes f and s stand for fluid and solid phases, respectively.
We pass to quantities assigned to the unit volume: Here, c = ρ p /ρ is the particle mass concentration and φ j is the volume fraction of the j-phase with j = f , p, M. It follows from the above definitions that the partial densities ρ j are related to the material densitiesρ j by the following equations Generally, the phase pressures p s and p f are different. However, as in [34], we assume that p s = p f = p. Such a hypothesis works well when the surface tension at the boundaries separating the phases are negligible.
Let v i , T i , l, k stand for the velocity, the viscous part of the stress tensor, the particle concentration flux vector and the interphase friction coefficient, respectively.
We introduce the tensor notations. Given two vectors a and b, we define the scalar product a · b = a i b i . The tensor product a ⊗ b is a matrix, such that (a ⊗ b) ij = a i b j . The matrix A * stands for the adjoint matrix of A, i.e., (A * ) ij = (A) ji . The i-th component of the vector div A is defined by the formula (div A) i = ∂A ik /∂x k .
If we neglect the rotation of particles and thermal effects, then one can derive from [34] the following mathematical model for six unknowns functions ρ s , ρ f , p, c, v s , v f : where p = p(ρ) is the prescribed state equation, g is the gravitation vector, and Now, we discuss rheology. Given a velocity field v(x), we introduce the corresponding rate of strain tensor D and its deviatoric part D d , as follows The fluid phase is considered to be a viscous Newtonian fluid. It implies that with the constant η f staying for shear viscosity. As for the compression viscosity, it is assumed negligible. The solid phase contains particles, with c being its mass concentration. Given a vanishing number δ and the volume fraction of the solid phase φ s , we define rheology of the solid phase by the regularized Bingham equation: Here, is the viscosity given by the Krieger-Douhgerty empirical closure [35], with φ * s and η 0 s being the maximal reference value of φ s and the consistency, respectively. As for for the yield stress, we take it by the correlation formula proposed in [36]. A mathematical proof is provided in [37] to ensure that Equation (7) is a good approximation for the Bingham fluid. One more rheological equation is the Fick law [33]: Due to the mass conservation laws (5), Equations (2) and (3) reduce to Such equations are of use in calculations performed below.
Let us formulate a hypothesis of incompressibility. We assume that the mud volume fraction φ M is negligible and the densities of materials ρ f , ρ p and ρ M are constants. Then, it follows from (1) that where R 0 =ρ s /ρ f . Observe that the total density ρ and the partial densities ρ j are not constant in contrast to the densities of materials. By the incompressibility assumption, one can derive easily the following formulas: The functions r s (c) and r f (c) are dimensionless partial densities. One more consequence of the incompressibility assumption is that the volumetric mean velocity is divergence-free: Equation (4) is equivalent to whereṽ is the mass mean velocity. Thus, we derived a mathematical model for four unknown functions p, c, v s v f , obeying the Equations (11), (12), (14) and (15). The parameters η s , η f , k, γ j are assumed to be know functions of the mass concentration c. Under the incompressibility hypothesis, pressure is no longer a thermodynamic parameter and does not satisfy the equation of state. It is now included in the list of unknown functions as in the case of Navier-Stokes models of a viscous incompressible fluid. Densities can be restored from equalities (13).
The diffusion coefficients γ j vanish when any phase disappears. As for the friction, we use the correlation formula proposed in [38] where d p is the particle diameter and C D is the particle/fluid friction: Now, we apply the described governing equations for 2D flows of suspensions in an oblique bifurcating channel ( Figure 1) representing a branching system, in which the parent branch divides into two daughter branches (one branch follows the inlet, termed as the main branch; and another bifurcate follows at an angle α with the main branch, termed as the side branch). This type of bifurcating channel is often encountered in the industry, nature, and human body, and one of the important tasks is to find out the bulk suspension and particle partitioning in the daughter branches for a better understanding of the flow behavior. On the (x, y) plane, the parent branch is directed along the x-axes, Figure 1. The side branch is at the angle α relative to the main branch. Let H, V, l 0 , t 0 andp stand for the reference values of the channel width, suspension velocity, particle concentration flux, time and pressure. We pass to dimensionless variables and choose the reference values in such a way that Flows are defined by the following demensionless numbers: Dimensionless stress tensors appear: (17) In the absence of any phase, the parameters γ 1 ,γ 2 , γ 3 disappear. This is why we impose the formulas We introduce the mean mass velocityṽ = cv s + (1 − c)v f . Let the differential operators stand for the material derivatives.
When omitting primes, we find that the functions where In what follows, we neglect the gravitation effect since we restrict ourselves to 2D horizontal flows.
The flow domain Ω is depicted in Figure 1. Let n stand for the unit outward normal vector to ∂Ω. The boundary conditions are formulated as follows. At the inlet boundary, The outlet boundary conditions are given as follows: At the impenetrable boundaries, we set The initial conditions are The functions standing in the right-hand sides of Equations (25)- (27) are assumed known.
To calculate the material derivative of a function f (x, t), we use the following approximation: where X m = X(t m ) with X(t) being the solution of the Cauchy problem Within FreeFem++, there is a procedure called "convect" which allows to solve this Cauchy problem and to determine the value f (X m (x), t m ): The material derivatives d f /dt and d s /dt are approximated in a similar way. Let us describe an algorithm.
Step 1. To tackle nonlinearity, we apply iterations as far as the concentration c is concerned. Given c m at the time level t m , we find c m+1 at the time level t m+1 by iterations. Given function c iter , we define the next iteration c iter+1 as a solution to the problem We take c m as the first iteration.
Step 2. We calculate the relative discrepancy Steps 1 and 2 are repeated until the condition E c < 0.01 is met. With the last iteration c f inal at hands, we define concentration at the level t m+1 by the equality c m+1 = c f inal .
Next, we solve the hydrodynamic part of equations. To this end, we use a modified SIMPLE method which is one of predictor-corrector approaches. At the prediction step, the guess velocities v * s and v * f are calculated, starting from the current pressure field p S . As an initial value for the field p S , we use p m . Then, a correction q to the current pressure is calculated with the help of the guess velocities. The final correction step consists of improving the guess velocities and the current pressure. The process is iterated until it converges, i.e., until the correction to pressure becomes small enough. There is a nonlinearity related to the convective terms in the momentum Equations (28) and (29). To tackle such nonlinearity, we use the simple iteration method.
Step 3. With the function c m+1 , p m , v m f and v m s being known, we look for p m+1 , v m+1 f and v m+1 s by iterations. To this end, we construct a finite sequence of guess pressures p S (S = 1, 2, · · · ) starting from p S = p m when S = 1. We associate with p S finite sequences of current velocities v iter f and v iter s (iter = 1, 2, · · · ), such that Given p S , v iter f , v iter s , we define v * f and v * s by solving the following system of equations: Step 4.
Given v * f and v * s , we find q iter as a solution of the equation Step 5. Now, we correct the current velocities, as follows: Step 6. We calculate the relative discrepancy of velocities If E v ≥ 0.05, we substitute v iter f and v iter s in (35) and (36) by v iter+1 f and v iter+1 s and repeat the procedure until the inequality E v < 0.05 is met. Then, iterations are recognized as convergent and we use the final iterations to define v SI MPLE Step 7. We calculate the relative discrepancy of pressure If E p ≥ 10 −3 , we substitute p S in (35) and (36) by p S+1 = p S + q f inal and repeat steps 3-7 until the inequality E p < 10 −3 is satisfied.
For the time level t m+1 , we define v m+1 provided that the condition E p < 10 −3 is satisfied.

Results
Our principle goal is to estimate the loss of particles into the side branch. We conclude that the greater the bifurcation angle, the smaller the particle loss. To justify this conclusion, we study the pulsatile mode of particles injection when the particle influx alternates with pure fluid influx. A pulsating injection finds many applications in the production of various materials [39].
We denote the total length of the channel AC ( Figure 1) by l. We assume that the side branch has the length l/2 and bifurcation of the horizontal parent branch is localized at point L, such that |BL| = l/2. Figure 2 depicts calculated snapshots of mass concentration c for the pulsating mode of particles injection, which corresponds to the following boundary and initial conditions: where t 1 = 50, t 2 = 60, x 1 = 0.01 · l, x 2 = 0.02 · l and the particle injection starts at t = 0 and stops at t = 60. Exactly at t = 60, the front of the particles reaches the location of bifurcation. In fact, the pulsating injection given by (41) and (42) is a regularization of the ideal pulsating injection, with the numbers 0.95 and 0.05 substituted with 1 and 0, respectively. Such a regularization facilitates calculations.  Figure 2 corresponds to data chosen as follows. We assume that all the branches have the same width W. We set To meet an agreement with no slip boundary conditions at rigid boundaries, we require that the inlet velocities v inlet f (t, y) and v inlet s (t, y) have the Poiseuille-like profile versus the vertical variable y [40].
Let us discuss findings and implications concerning the results in Figure 2. First of all, the partitioning of particles occurs at the bifurcation of the channel and one can observe the propagation of the concentration wave both in the side and main branches. We emphasize that the particle mass concentration c obeys the inequalities 0 < c < 1 in spite of the fact that it is governed not by the common transport equation but satisfies a complicated diffusion equation with a matrix diffusion coefficient. Observe that we do not use any numerical cutting tricks to ensure the inclusion 0 < c < 1. The concentration front is not blunted and there is a cusp (inverted V-shape) in the middle in agreement with experimental data [41]. For comparison, we note that the calculations in [32] do not reveal a cusp in the concentration front. The following explanation for the front with cusp can be given. The particle concentration at the center is higher compared to the near-wall regions because the particles migrate from the region of a high shear rate (wall) to the region of a low shear rate (center). Within the framework of our model, estimates show that the transverse phase velocities are very small and the lateral migration of particles is mainly due to diffusion. Since the role of diffusion is significant, there are limitations on the scope. The results are valid for sufficiently small particles.
Let us introduce particle mass fluxes through the cross-sections KL (inlet), MN (outlet) and LN (branch): where the dimensionless total mass velocity j is equal to r f v f + r s v s . We observe that such mass fluxes include both a convective and diffusive particle discharge. By choosing the data (43), we performed a calculation of the introduced fluxes on a pulsation time interval (0, T) great enough to ensure that Q Partitioning of particles between the branches is shown in Figure 3 for Q p in (t), Q p out (t), Q p br (t) and in Figure 4 for M p in , M p out , M p br . One can see that particle mass loss into the side branch occurs. Such a result explains that particle separation in branched channels can really happen due to hydrodynamic forces only. The front of particles reaches the channelbranching location at a moment close to t = 100. This is why the mass flow peaks at the cross sections of branches near the bifurcation location occur on the interval 100 < t < 120. The mathematical model developed here enables one to optimize the separation effect by the variation of geometric, fluidic and other data. Figure 4 is a result of integration of the functions at Figure 3 over the time interval 0 < t < T, T = 400. Figure 5 shows dynamics of the mean mass concentration of particles through the inlet (KL), outlet (MN) and branch (EF) cross-sections. Almost equal fluxes for the great times into the daughter branches bifurcating at the angle 30 • are due to equality of their widths. Figure 2 shows the qualitative partitioning of particles at the branching point of the channel with the bifurcation angle at 30 degrees, while Figures 3-5 depict the quantitative partitioning. Although the side and main branches have the same width, the flow of particles into the side branch is slightly less. Apparently, this is the effect of inertia. In most papers, attempts to explain the separation of particles are associated with the study of the behavior of streamlines at the bifurcation zone [42]. Our approach is different since we consider unsteady flows. Moreover, we take into account the difference in the velocities of solid particles and fluid. In such a case, the notion of streamline is not applicable. Figure 6 proves that an increase in the bifurcation angle results in a decrease in the loss of particles into the side branch. A similar conclusion is derived by simulation in [32]. We     (43). The red, blue and green lines are for the inlet cross-section KL, for the outlet cross-section MN, and for the branch cross-section EF, respectively, (see Figure 1).  (41) and (42), we conclude that there is a correlation between the relative values of mass f and p given at Figure 7. One can observe that an increase in f causes the value of p to increase; this is also in agreement with the available experiment data [43,44]. Although a quantitative comparison with these experiments is not possible due to the 2D assumption of the present simulations, the results in Figure 7 fairly reproduce the experimental trends. We also should note that these experimental results concern very dilute suspensions with c 0.02.
We verified that inequality holds for the data used in this paper. Such a result agrees with the Zweifach-Fung effect [45], amounting to the fact that more particles fall into the branch where more fluid enters. As far as the Zweifach-Fung effect is concerned, we would prefer not only to validate our model, but to address some questions which were not captured by simulations in other publications [45]. First of all, we are interested in the pulsatile mode of injection of particles, while other authors study stationary regimes. The next question concerns particle density. Our method is suitable for both heavy and light particles. All known calculations of flows in the branching channels concern neutrally buoyant particles. We fill this gap. Within our model, the principle polymer parameters are the viscosity and the yield stress of the solid phase since we treat this phase as a viscoplastic Bingham fluid. Calculations for α = 90 • reveal that pressure can be assumed constant along any cross-section of every branch. Hence, the pressure gradient can be measured along the center line of the branch, |∇ p| |∂p/∂s|, with s standing for the branch-length reckoning from the bifurcation point (L for the side branch and N for the main branch on Figure 1). On the other hand, it follows from Figure 8 that pressure decreases in s at any time instant. Let |∇ br p| and |∇ out p| stand for the pressure drop along the side and main branches: where l br = LF and l out = |ND| are the branch lengths. One can see that l out < l br since |LN| + |ND| = l/2 (see Figure 1). We introduce the following dimensionless branch characteristic in dimension variables: Passing to dimensionless variables and omitting the prime indexes, we obtain that where w k is the branch width in dimensionless variables. The parameter b n coincides with the Bingham number Bn for the pressure driven flows in a simple channel [46]. This is why we also call here b N the reduced branch Bingham number. The reduced branch Bingham number b N depends on the branch width and the pressure drop along the branch. Let us consider variation of w br for α = 90 • , keeping w out constant and equal to 1. One can see from Figure 8 that a change in the side branch width implies a change in the pressure drop along the branch. Since both the branches have almost the same lengths, one can also conclude from Figure 8 that the pressure drops |∇ br p| and |∇ out p| in both the branches are almost equal and do not depend on the variation of the side branch width. Thus, the branch Bingham number depends mainly on the branch width.
Let us perform the results of the calculations of the time-average values b N of the reduced Bingham number, Table 1. We remark that the pulsation interval 0 < t < T is chosen in such a way that all the channel is filled with pure fluid, at t = T. Calculations reveal that T = 400 if w br = w out = w in = 1. The terminal value T increases if w br decreases.
When we address the particles partitioning in relation to the change in w br , as in Table 1, we arrive at the results depicted in Figures 9 and 10. One can see that the concentration wave becomes slower when we reduce w br from 0.5 to 0.2. The total particle mass flux also reduces. As far as the polymer particles are concerned in the pulsatile injection mode, we arrive at the following conclusion: more particles fall into the branch with lower mean Bingham number. One of the implications of this result is the stoppage effect. It implies that there is no flow in the branch if the pressure gradient along the branch is small enough or the channel is very narrow. Indeed, Figure 11 shows that the solid phase velocity in the side branch is two orders of magnitude less than the solid phase velocity in the main branch in the case when the side branch is five times narrower than the main branch. Due to the regularization of Equation (7), the velocity of the solid phase in the side branch can be considered zero. The stoppage effect is well known in the Bingham fluid theory in more general context. If the pressure gradient applied in fully-developed Newtonian Poiseuille flow is suddenly set to zero, the velocity decays to zero exponentially, i.e., the theoretical stopping time is infinite [47]. This is not the case for viscoplastic or yield-stress materials [48]. The stopping time is finite.  Figures 11 and 12 show the phase velocity profiles at different cross-sections of both the branches at the terminal time instant T. One can see that there is a correlation between the Bingham number and velocity. For w br = 0.5 and w br = 0.2, both the solid and fluid phases flow slower in the side branch than in the main branch of the width w = 1. The side branch velocity decreases significantly if the width of this branch decreases. One more observation is that the velocities of both phases are almost equal due to the interphaseresistivity effect. Note, that the velocity in the main branch is indifferent to a decrease in width of the side branch. It is well known that non-Newtonian fluids in channels have a blunt Poiseuille profile [47]. The calculation results concerning the velocity field are shown in Figures 11 and 12. Thus, calculated profiles are consistent with the Poiseuille-like flows and validate the numerical scheme developed in the present paper.
To show the role of the yield stress, we performed calculations with τ * = 0, Figures 13 and 14. One can see that under the conditions (43), velocities in the side branch are doubled but velocities in the main branch are unchanged.

Remark 1.
Let us comment on the validation of the numerical simulation performed in the present paper. First of all, our calculations capture the benchmark experimental Zweifach-Fung effect, stating that more particles fall into the branch where more fluid enters. A qualitative agreement with 3D experiment data on the particles partitioning is established (see Figure 7 and [43,44]) by calculations relative to the 2D branching tube. We use the same numerical algorithm as in the recent work on sedimentation [33]. The difference is only in the boundary conditions. As applied to sedimentation, our calculations also explain the Boycott effect, amounting to the fact that an inclination of a vertical vessel enhances the sedimentation of particles in suspension. It is well known that non-Newtonian fluids in channels have a blunt Poiseuille profile. Figures 11 and 12 agree with such a flow property.
In order to illustrate the convergence of the method, we performed a series of calculations for different meshes with n vertices on the border AB, 4n, vertices on the borders LF, NE, BL, ND, and 8n vertices on the border AC for n = 8, 16, 32 and 64.
This gives meshes with a total number of N = 952, 3448, 15,058 and 56,544 vertices, respectively.
We denote the numerical solution f = (v s , v f , p, c) corresponding to n by f n . Convergence manifests itself through the following estimates where r n = 1.6, 1.99 for n = 8, 16, respectively. This estimate indicates that the solution is grid-independent. To guarantee convergence and an acceptable computation time, we use the mesh with at least 10 3 triangular elements.

Conclusions
As in a number of studies, we consider a dense suspension of polymer particles to be a non-Newtonian fluid. Since we assume that the particles and the carrier fluid have different densities and velocities, we use a two-continuum model. The first phase of solid particles is described by the rheology of the Bingham fluid and the second phase is a viscous Newtonian fluid. There is one more rheological feature of the solid phase of particles which is associated with the generalized Fick's law for the particle concentration flux vector. The fact is that in the full model, in addition to ordinary diffusion, barodiffusion and thermal diffusion are taken into account. In addition, the concentration flux vector also takes into account a component that depends on the gradient of the modulus of differences in phase velocities. The model agrees with the basic principles of thermodynamics and is validated through capturing the Boycott sedimentation effect. Starting from the chosen model, we address flows in a branching channel with the rather arbitrary bifurcation angle, which is reckoned from the inlet direction. We study the issue of flow partitioning and estimate the loss of particles into the side branch during the pulsatile injection of particles. We prove that the greater the bifurcation angle, the smaller the loss of particles. Our calculation of the loss agrees qualitatively with experiment data. We prove that the partitioning of particles occurs in agreement with the Zweifach-Fung effect, stating that particles prefer the branch with a higher fluid flow rate. We establish that more particles fall into the branch with a lower mean Bingham number. The results are applicable to the technology of producing microfluidic devices consisting of tubes with branches. One more application concerns the simulation of the proppant particle loss in perforations during the proppant delivery to a hydro-fracture. Now that the model has been tested on the Zweifach-Fung benchmark effect, it becomes possible to use it to study the flows of suspensions of polymer particles in confluences.