Enhanced Numerical Modeling of Non-Newtonian Particle-Laden Flows: Insights from the Carreau–Yasuda Model in Circular Tubes

TL;DR

This paper improves the modeling of particle movement in complex fluids using the Carreau–Yasuda model, leading to better predictions for industrial applications.

Contribution

The study introduces the Carreau–Yasuda model for non-Newtonian particle-laden flows, capturing shear-thinning and Newtonian behavior.

Findings

The Carreau–Yasuda model significantly alters particle migration dynamics compared to simpler models.

Parametric studies show how flow conditions and model parameters affect particle concentration and velocity profiles.

Validation with experimental data confirms the model's accuracy in predicting suspension behavior.

Abstract

Particle-laden flows in non-Newtonian fluids are encountered in a variety of industrial applications, such as concrete pumping and battery electrode slurry processing, where accurate prediction of particle migration is essential for performance and product quality. This work investigates fully developed suspension flows in circular tubes, combining the shear-induced diffusion framework of Phillips et al. with the Krieger–Dougherty relative viscosity and the Carreau–Yasuda constitutive model. Unlike previous studies that generally rely on Newtonian or simple non-Newtonian rheology models, we employ the Carreau–Yasuda model, a more sophisticated constitutive relation that captures both shear-thinning behavior and Newtonian plateau regimes. By applying nondimensionalization and variable transformations, we reduce the governing coupled differential equations to a system of nonlinear algebraic equations, which allows for efficient computation of both particle concentration and velocity profiles. A systematic parametric study was conducted to evaluate the influence of several factors, including the pressure gradient, average particle concentration, and the five parameters of the Carreau–Yasuda model. Additionally, the migration parameter α=Kc/Kη was considered. The results reveal how the increased rheological complexity of the Carreau–Yasuda model significantly alters the migration dynamics when compared to simpler models. These novel findings have direct implications for optimizing industrial processes involving highly loaded suspensions, offering more accurate predictions of particle behavior under varying flow conditions. For the validation of our findings, experimental data in the literature was used.

Figures

Full text

1. Introduction

Fluid flows with suspended particles are encountered in many industrial applications, such as concrete pumping [1] and extrusion and injection molding processes with filled polymers [2,3]. Understanding the fundamentals of these flows through numerical modeling is therefore essential. In particular, it is important to investigate the behavior of suspension flows using generalized rheological models, which offer a broader range of applicability.

In this work, we focus on steady-state flows in circular tubes with non-Newtonian fluids. We present a simplified approach to numerically modeling the particle volume fraction and velocity distributions based on a reduced form of the shear-induced diffusion equation. We then examine the influence of both physical parameters and the rheological parameters of the Carreau–Yasuda model on concentration and velocity profiles.

One important application of such studies is the rheology of electrode slurries, which has a critical impact on the final coating microstructure. Challenges include high slurry viscosity, which increases pumping pressure and limits coating speed; elasticity, which can lead to instabilities and defects; and excessive flow, which may result in slumping and poorly structured coatings. These complexities highlight the importance of optimizing slurry formulations. However, the wide range of solvent systems and components involved makes it difficult to identify the exact origins of these unfavorable rheological properties [4]. Furthermore, precise measurements of extensional properties using miniature rheometer have revealed significant variability in slurry behavior across different formulations and binder systems [5].

The work of Phillips et al. [6] provides a fundamental reference for this study. They derived the shear-induced diffusion equation for the particle volume fraction ( $[eqn]$ ) and established the corresponding no-flux boundary condition. A key difference between our work and theirs is that Phillips et al. assumed Newtonian fluids in Couette and Poiseuille flows, whereas we consider non-Newtonian flows governed by the Carreau–Yasuda model [7].

Following the seminal work of Phillips et al. [6], a substantial body of literature has focused on the theoretical and numerical modeling of shear-induced particle migration without direct experimental validation, largely due to the limited availability of comprehensive experimental data for complex suspension flows. Such pure modeling studies have nonetheless played a central role in advancing the understanding of suspension dynamics. For instance, Kang and Mirbod [8] employed direct numerical simulations to investigate shear-induced migration in Poiseuille and circular Couette flows of semi-dilute and concentrated Brownian suspensions with shear-thinning behavior. Dbouk et al. [9] studied particle migration using the Suspension Balance Model and finite-volume simulations, examining the influence of particle stress formulations on migration predictions. More recently, Kauzlarić et al. [10] proposed a non-local extension of the Phillips model, demonstrating how alternative modeling assumptions affect concentration profiles. These works illustrate that modeling-based investigations, even in the absence of new experimental data, constitute a well-established and widely accepted approach in the study of shear-induced particle migration.

Earlier, Chen et al. [11] investigated shear-induced particle migration in non-Newtonian flows of nickel-powder-filled polymers. They employed the power-law model to describe the viscosity–shear rate relation and incorporated temperature dependence using the Arrhenius law. Solving the coupled differential equations for velocity and particle volume fraction proved to be nontrivial. Subsequently, Chen et al. [12] developed a numerical scheme to obtain these solutions for both steady and unsteady cases. Their method combined the Newton–Raphson iterative approach for velocity fields with a finite-difference semi-implicit scheme (Crank–Nicolson [13]) for the diffusion equation.

In contrast, we applied nondimensionalization and variable transformations, following techniques presented by Wang [14], to reduce the governing system of differential equations to a system of nonlinear algebraic equations. The solutions were then obtained numerically, which is computationally less intensive than differential equation–based approaches. Moreover, the Carreau–Yasuda model provides a more general description compared to the power-law model.

The originality of the present work lies in (i) the incorporation of the Carreau–Yasuda constitutive model within the shear-induced diffusion framework, enabling the treatment of shear-thinning fluids with finite viscosity plateaus; (ii) the reduction of the governing equations to a system of nonlinear algebraic equations through nondimensionalization and variable transformation, resulting in a computationally efficient solution strategy (see Appendix A); and (iii) a systematic parametric investigation encompassing physical parameters, particle migration coefficients, and all rheological parameters of the Carreau–Yasuda model (see Results and Discussion section and Appendix B).

Several related studies have examined particle migration in different flow configurations. Choi et al. [15] studied concrete pumping in pipes using the Bingham fluid model for the viscosity–shear rate relation and commercial CFD software (Fluent) for numerical solutions. Siqueira and Carvalho [16] investigated particle migration in planar extrudate flows of suspensions of hard spheres. They analyzed Newtonian flows between parallel plates and subsequent free-surface extrudate flows, employing a stabilized finite element method along with elliptic mesh generation to compute the free surface. Rebouças et al. [17] studied fully developed suspension flows in circular tubes, focusing on the effects of particle and shear-induced viscosity. Their work is closely related to ours, as they also solved the shear-induced diffusion equation for non-Newtonian flows. They examined the effects of the average particle volume fraction on concentration and velocity profiles, and their results agree qualitatively with our findings. However, their approach addressed the singularity near the tube centerline and relied on the Cross model for fluid viscosity, whereas we adopt the more general Carreau–Yasuda model, which includes an additional parameter (a). Furthermore, we systematically explore the influence of a broader set of parameters, including the pressure gradient, average particle volume fraction, all five Carreau–Yasuda parameters, and the model parameter  $[eqn]$  (hereafter  $[eqn]$ ).

In addition to numerical studies, Rueda [18] presented a comprehensive experimental investigation of the rheological properties of suspension flows, examining how particle size, concentration, and fiber interactions influence viscosity and flow regimes. His results highlight the sensitivity of suspension rheology to microstructural features, providing valuable experimental context for the numerical analysis carried out in the present work. In particular, these findings reinforce the importance of understanding how rheological parameters and particle migration mechanisms shape the concentration and velocity profiles studied here.

It is also important to emphasize that our approach has been verified against existing models in the literature. Specifically, back-substitution of the computed particle volume fraction and shear-rate distributions into the governing equations yields negligible residual errors, confirming the correctness of the formulation and implementation. Moreover, qualitative agreement with the results of Rebouças et al. [17] further supports the validity of our model. To further validate our numerical framework, we compared our predicted velocity profiles against experimental data from Preziosi et al. [19] for low-viscosity emulsions flowing in a confined capillary. The comparison shows good agreement, with a mean relative error of less than 8%, thereby demonstrating the applicability of our Carreau–Yasuda-based model to real experimental systems.

The remainder of this paper is organized as follows. Section 2.1, Section 2.2 and Section 2.3 present the mathematical formulation of the problem, including the shear-induced diffusion model and the Carreau–Yasuda constitutive relation. Section 2.4 describes the nondimensionalization, numerical scheme, and iterative algorithm used to obtain the solutions. Section 3 discusses the numerical results, highlighting the effects of physical and rheological parameters on particle concentration and velocity profiles. Finally, Section 4 summarizes the main findings and outlines the implications of this modeling framework for future studies and applications.

2. Model and Methods

2.1. Governing Equations

The flow is governed by equation of continuity for an incompressible fluid

$$[eqn]$$

the conservation equation of linear momentum

$$[eqn]$$

and Newton’s Law of viscosity

$$[eqn]$$

2.2. Particle Poiseuille Flow Through a Cylindrical Tube

We consider the flow of a particle-laden non-Newtonian suspension in a circular tube under steady pressure-driven conditions (Poiseuille flow). The central modeling challenge is to describe how particles migrate across streamlines due to hydrodynamic interactions and spatial variations in viscosity. To this end, we adopt a shear-induced diffusion framework originally proposed by Phillips et al. [6], which balances the competing effects of particle collisions, viscosity gradients, and Brownian motion. Figure 1. is used as a schematic illustration of the anticipated results of our numerical simulation.

Let  $[eqn]$  denote the particle volume fraction. To describe the rheology of the suspension, we combine a concentration-dependent relative viscosity model with a shear-rate-dependent constitutive law for the carrier fluid. Consider the Krieger & Dougherty [20] model, which is typically applied for the particle size range between nano- and micrometers, for laminar flow with the following terms on particle-volume-induced relative viscosity:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Fluxes  $[eqn]$  and  $[eqn]$  represent shear-induced migration mechanisms introduced by Phillips et al. [6], while  $[eqn]$  accounts for random particle motion due to Brownian effects. The symbol b in the particle flux terms  $[eqn]$  and  $[eqn]$  denotes the particle radius (m). Unless otherwise indicated,  $[eqn]$  and  $[eqn]$  are taken as dimensionless migration coefficients (Phillips et al. formulation) so that the combination  $[eqn]$  and  $[eqn]$  has dimensions of  $[eqn]$  (as needed for a diffusive flux).

The Brownian diffusion coefficient D has units of  $[eqn]$ , and the Peclet number is defined as

$$[eqn]$$

This number  $[eqn]$  measures the relative strength of shear-induced transport compared to Brownian motion. Thus, the conservation equation for the particle volume fraction becomes

$$[eqn]$$

In this work we adopt the steady-state assumption  $[eqn]$ , meaning that the Brownian flux  $[eqn]$  is negligible compared to the shear-induced fluxes  $[eqn]$  and  $[eqn]$ . For clarity, we explicitly state the range of Peclet numbers for which this approximation is valid. The Peclet number quantifies the competition between shear-driven migration and Brownian diffusion. For colloidal particles, Brownian contributions are significant when  $[eqn]$  and become negligible only when  $[eqn]$ – $[eqn]$ , as reported in Phillips et al. [6] and subsequent experimental studies on shear-induced migration. For the particle sizes and shear rates relevant to the systems modeled here (typically  $[eqn]$ – $[eqn]$  and  $[eqn]$ – $[eqn]$ ), the resulting Peclet numbers fall in the range  $[eqn]$ – $[eqn]$ . Thus, Brownian diffusion is several orders of magnitude weaker than the hydrodynamic migration mechanisms, and neglecting  $[eqn]$  is fully justified for the parameter regime considered in this study. Thus, replacing the flux terms in (10) with Equations (4)–(9) gives us

$$[eqn]$$

For flows in a cylindrical tube, the gradient can be written in  $[eqn]$  in general form. At the tube wall, particle conservation requires that the net flux is normal to the boundary vanishes, leading to the no-flux condition

$$[eqn]$$

So, combining (11), (12), and the governing Equations (1)–(3), we get the generalized system of differential equations of velocity profile and particle volume fraction distribution, which is

$$[eqn]$$

with the no-flux boundary condition at the wall

$$[eqn]$$

We note that the phrase “Newton’s law of viscosity” in the text is used to refer to the kinematic relation between shear stress and shear rate. Strictly speaking, Newton’s law requires a constant viscosity (Newtonian fluid),  $[eqn]$  with  $[eqn]$  const. In this work we instead adopt the generalized Newtonian form

$$[eqn]$$

where the scalar viscosity  $[eqn]$  depends on the local shear rate and particle concentration through the constitutive relations (Krieger–Dougherty for the concentration dependence and Carreau–Yasuda for the shear-rate dependence). Thus there is no inconsistency: the momentum equation, which is the first equation in (13), is written for the generalized-Newtonian stress  $[eqn]$ .

Regarding the pressure gradient appearing first equation in (13), in the fully developed pipe flows considered here, the axial pressure gradient  $[eqn]$  is spatially uniform and acts as the external driving force. In our computations it is treated as the prescribed driving parameter (related to the imposed pressure drop  $[eqn]$  over the pipe length L), but equivalently one can prescribe a total flow rate and determine the corresponding  $[eqn]$ . The wall force balance then yields

$$[eqn]$$

which connects the (dimensional) pressure gradient to the wall shear rate used for nondimensionalization.

2.3. Steady State Velocity and Particle Distributions

In steady state flow,  $[eqn]$ , and we assume the r-direction, so the second equation in the coupled system (13) becomes

$$[eqn]$$

which can be integrated over r, and applying the no-flux boundary condition to compute the constant of integration, we get

$$[eqn]$$

which means no flux everywhere in the suspension in the steady-state scenario. This equation can be integrated with the boundary condition at the wall (subindex w), so we get the following simplification

$$[eqn]$$

Now we want to simplify the system. For steady, axisymmetric, fully developed flow in a cylindrical tube (with  $[eqn]$ ,  $[eqn]$ ), the radial momentum balance reduces to  $[eqn]$ . Hence the pressure depends only on the axial coordinate,  $[eqn]$ . By rearrangement of the first part of the coupled system (13), we get

$$[eqn]$$

and integration of (18) gives us

$$[eqn]$$

Applying the regularity condition at  $[eqn]$  gives  $[eqn]$ . Then, we have the simplified equation

$$[eqn]$$

Summarizing the coupled ODE for steady-state case, we have the following system of differential equations (see details in Appendix C)

$$[eqn]$$

2.4. Nondimensionalization and Numerical Details

In order to further simplify the coupled system (21), we introduce the nondimensional variables as follows

$$[eqn]$$

Substituting the new variables (22) into (21) gives us the system in nondimensional form as

$$[eqn]$$

Note that, in our formulation, the quantity  $[eqn]$  appears only as a constant driving term, and the equation of (23) is nondimensional relative to the other variables because the combination  $[eqn]$  naturally scales with the characteristic shear stress  $[eqn]$ . Since we have  $[eqn]$ ,  $[eqn]$ , and  $[eqn]$  at the wall, the first equation in (23) becomes

$$[eqn]$$

where  $[eqn]$  is the magnitude of the pressure gradient. Suppose we use the Krieger–Dougherty equation for  $[eqn]$  and the Carreau–Yasuda model for  $[eqn]$  as in Equation (6), we have

$$[eqn]$$

Combining Equations (24) and (25) gives

$$[eqn]$$

The value  $[eqn]$  is used in the Krieger–Dougherty relation because it represents the maximum volume fraction attainable in random close packing (RCP) of spherical particles. This is the limit beyond which the viscosity diverges due to geometric crowding, preventing further flow. This value is physically realistic for disordered suspensions under shear and aligns with experimental observations of dense colloidal and granular flows. Using  $[eqn]$  ensures that the model remains consistent with the assumption of a non-crystalline, randomly packed suspension. Given  $[eqn]$ ,  $[eqn]$ , and R, we can obtain  $[eqn]$  by solving the above nonlinear Equation (Equation (26)) for  $[eqn]$ . Once  $[eqn]$  is obtained, then it is straightforward to obtain  $[eqn]$ , which is used as the unit for viscosity. Using the identity shown in Equation (24), we can rewrite Equation (23) as

$$[eqn]$$ $$[eqn]$$

where  $[eqn]$  and is taken to be  $[eqn]$  in this work, following Phillips et al. [6]. Again, using the Krieger–Dougherty equation for  $[eqn]$  and the Carreau–Yasuda model for  $[eqn]$ , we have

$$[eqn]$$

where

$$[eqn]$$

Here,  $[eqn]$ ,  $[eqn]$ , and  $[eqn]$ . The  $[eqn]$  term is the time constant  $[eqn]$  in dimensionless form, using  $[eqn]$  as the characteristic time scale. Combining Equations (27) and (28), we have

$$[eqn]$$

and combining Equations (27)–(29), we have

$$[eqn]$$

From the above two equations, Equations (31) and (32), one can easily obtain the volume fraction  $[eqn]$ , the shear rate distribution  $[eqn]$ , and the volume fraction as a function of the shear rate,  $[eqn]$ . The velocity distribution can be obtained using the method of Wang [14] as follows. Using  $[eqn]$  as the unit for velocity, i.e.,  $[eqn]$ , we have

$$[eqn]$$

Since the method of Wang [14] suggests  $[eqn]$ , the velocity profile  $[eqn]$  is obtained by integrating the corresponding ordinary differential equation numerically in the radial direction. In this work, we employ the classical fourth-order Runge–Kutta method with a uniform radial grid and the boundary condition  $[eqn]$  at the tube centerline. This approach provides stable and accurate integration for smooth profiles and is widely used for ODEs of this type.

If we try to solve the system, consisting of Equations (31) and (32), for s, complications in the numerical computations arise, for instance,  $[eqn]$  is highly nonlinear. Thus, following the idea from the paper of Wang [14], it is wise to solve for  $[eqn]$  instead of s. This means, for any s∈ $[eqn]$ , we want to solve for  $[eqn]$  and  $[eqn]$ , which can be easily obtained by root finding ( $[eqn]$ ) in the following equation

$$[eqn]$$

and substituting the solution back to one of Equations (31) or (32) to find  $[eqn]$ . Because  $[eqn]$  is required to solve the system, with given scalar  $[eqn]$ , it can be estimated by using the concept of binary search or the so-called “interval halving method”. Since  $[eqn]$  is always less than the  $[eqn]$ , we know that the true  $[eqn]$  is greater than 0 and less than  $[eqn]$ . By using this information, we can find the the corresponding  $[eqn]$  that is shown as an intermediate step in Figure 2.

3. Results and Discussion

3.1. Verification and Validation of Numerical Results

The Carreau–Yasuda rheology enhances particle migration toward the tube centerline, producing stronger core accumulation than in Newtonian suspensions. We aimed to generate some output by using the algorithm shown in Figure 2. The Python programming language was used to implement the numerical algorithm. We assume that the fluid viscosity is described by the Carreau–Yasuda model. The following table represents the constant values we chose to use as a base in the computations. Figure 3 and Figure 4 show the numerical solutions of the system.

Remark **1.**The Carreau–Yasuda parameters  $[eqn]$ ,  $[eqn]$ , λ, a, and n in Table 1 are adopted from the representative example used by Wang et al. [14], where the reference values  $[eqn]$ ,  $[eqn]$ , and  $[eqn]$  are taken directly from their study; in Wang et al.,  $[eqn]$ , but a nonzero value introduces no analytical or numerical difficulty, so we adopt  $[eqn]$  as a physically plausible lower-shear viscosity plateau. The maximum packing fraction  $[eqn]$  represents a typical limit for concentrated suspensions. The ratio  $[eqn]$  is set to  $[eqn]$  following the original migration-coefficient values proposed by Phillips et al. [6]. The remaining parameters, including the imposed pressure gradient and radius, are chosen to provide a physically reasonable flow regime for illustrating the numerical behavior of the coupled model in Figure 3 and Figure 4.

The solutions of particle volume fraction ( $[eqn]$ ) and the dimensionless shear rate (s) distributions are shown in Figure 3a,b. The dimensionless velocity ( $[eqn]$ ) profile in Figure 4a is computed from the integration of the simple ODE (33). The iterative behavior of the inner algorithm that finds the best approximation of  $[eqn]$  is shown in Figure 4b.

To verify the numerical model and the corresponding Python code used to solve for  $[eqn]$  and s, we back-substituted these solutions into the system (Equations (31) and (32), labeled as  $[eqn]$  and  $[eqn]$ ) and confirmed that the residual errors are on the order of  $[eqn]$ . This provides confidence in the numerical accuracy.

Unless otherwise specified, we focus on numerical accuracy and parametric trends, while physical interpretations follow the framework outlined at the beginning of this section.

3.1.1. Newtonian Fluid

In this section, the validation of our numerical algorithm in comparison with Phillips et al.’s [6] analytical solution is presented. The authors of the analytical solution of the equation governing  $[eqn]$  in the steady state presented the explicit formula for  $[eqn]$  as

$$[eqn]$$

where  $[eqn]$ . This explicit formula describes the particle volume fraction distribution in the steady-state Newtonian case. We compared our numerical algorithm with the analytical solution for different  $[eqn]$  and the results are shown in Figure 5.

We can see that our numerical solution perfectly agrees with their results.

3.1.2. Power-Law Fluid

Next, let us consider the viscosity defined by the power-law model as follows

$$[eqn]$$

where K is the fitted parameter and n is the power-law index. We incorporate this equation into our numerical algorithm to solve for concentration profiles for different values of power-law index n. To validate our numerical approach, let us derive the analytical solution in a similar way presented by Phillips et al. [6]. From Equation (20) we have

$$[eqn]$$

and substituting this formula into Equation (17) gives us

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

By imposing the relation

$$[eqn]$$

we obtain the following analytical expression for the fully developed, steady-state power-law fluid:

$$[eqn]$$

The choice in Equation (43) is a specific mathematical condition used to simplify the steady-state shear-induced migration equation and obtain an explicit analytical expression for the particle volume fraction in a power-law fluid. Physically, this corresponds to a situation where the migration flux due to gradients in viscosity ( $[eqn]$ ) and the collision-induced flux ( $[eqn]$ ) are related by this particular ratio. While such a precise ratio may not exactly correspond to a specific experimental system, it provides a convenient benchmark for comparison with numerical results and illustrates the influence of the Carreau–Yasuda rheology and particle interactions on the steady-state profile.

Figure 6 shows both the analytical solution from the above derivation and our numerical solution for different power-law indices n. Interestingly, when  $[eqn]$ , that is  $[eqn]$ ,  $[eqn]$  vanishes, and the solution of Equation (43) is constant for any value of  $[eqn]$ . For values  $[eqn]$  the concentration profiles are higher near the centerline and lower near the wall of the circular tube; however, if  $[eqn]$  we observe a very small particle volume fraction near the centerline, and it increases towards the wall. From this special case, we can conclude that for  $[eqn]$  and  $[eqn]$ , the particle volume fraction distribution is exactly flat at  $[eqn]$ . Generally speaking, when  $[eqn]$  in Equation (43), the concentration profile is  $[eqn]$  along the radial position.

3.1.3. Comparison with Experimental Data

To validate the predictive capability of our numerical model, we compared the simulated velocity profile against experimental measurements reported by Preziosi et al. [19] for a low-viscosity emulsion flowing in a confined capillary. In their work, emulsions consisting of silicone oil droplets dispersed in a Boger fluid (high elasticity, HE) were studied under confined Poiseuille flow. Velocity profiles were obtained via high-speed microscopy and droplet tracking. We fitted the Carreau–Yasuda model parameters to the rheological data of the HE fluid provided by Preziosi et al. [19] (see Figure 7a) and used the following parameters in our simulation:

• $[eqn]$  Pa·s (zero-shear viscosity);
• $[eqn]$  Pa·s (infinite-shear viscosity);
• $[eqn]$  (Carreau–Yasuda transition parameter);
• $[eqn]$  s (relaxation time);
• $[eqn]$  (power-law index);
• $[eqn]$  Pa/m (pressure gradient);
• $[eqn]$  m (tube radius);
• $[eqn]$  (maximum particle volume fraction used by Preziosi et al.);
• $[eqn]$  (migration parameter  $[eqn]$ );
• $[eqn]$  (average particle volume fraction).

Figure 7b shows the comparison between the experimentally measured velocity profile and the profile predicted by our model. The simulation reproduces the experimental trend with good agreement, particularly in the central region of the capillary where the velocity profile exhibits a blunted, plug-like shape due to droplet accumulation. The mean relative error between the predicted and experimental velocities is 7.78%, indicating that the Carreau–Yasuda model combined with the shear-induced migration framework captures the essential rheological and migration behavior of the emulsion. Discrepancies near the wall may be attributed to experimental uncertainties, the simplified assumption of a homogeneous suspension, and possible entry effects not accounted for in our fully developed flow model. Nonetheless, the overall consistency supports the use of our approach for predicting confined suspension flows with elastic matrices.

3.1.4. Other Relative Viscosity Models

Now let us test our numerical algorithm for various relative viscosity models other than Krieger–Dougherty [20]. The models of Chong (1971) [21] and Quemada (1977) [22] were compared and demonstrated the results of velocity and concentration profiles as seen in Figure 8.

3.2. Influence of Average Particle Volume Fraction

Now we analyze the influence of the parameter change on the output. In particular, we will alter one parameter at a time while keeping everything else fixed (as in Table 1) and observe the change in the output. The parameters that will be altered are the average particle volume fraction ( $[eqn]$ ), pressure drop ( $[eqn]$ ), the Carreau–Yasuda model parameters, and the  $[eqn]$ . First, let us consider the effect of average volume fraction on the output results.

We change the values of the average particle volume fraction from 0.1 to 0.4 while keeping other parameters fixed (as in Table 1). In Figure 9a, the concentration profile moves upward as  $[eqn]$  increases. This observation is as expected given the result in Figure 9b where the relation between  $[eqn]$  and  $[eqn]$  lies below the diagonal line, which agrees with the theory that  $[eqn]$ . We can also see that the velocity profile decreases with the increasing average particle volume fraction  $[eqn]$  (Figure 9c). Correspondingly, the average velocity decreases as well (Figure 9d). In general, the plots in Figure 9 suggest that the particle volume fraction and velocity distribution are highly dependent on the average particle volume fraction.

3.3. Influence of Pressure Gradient

Next, we wanted to determine the influence of pressure drop on the output, so we varied it from  $[eqn]$  Pa/m to  $[eqn]$  Pa/m and kept everything else constant (as in Table 1).

The behavior of the concentration profile is not uniform since, in Figure 10a, it decreases near the centerline ( $[eqn]$ ) but increases near the wall ( $[eqn]$ ) as the pressure drop increases. Clearly, the  $[eqn]$  vs.  $[eqn]$  curve moves upwards as shown in Figure 10b, and the velocity profile increases as well as the average velocity (Figure 11a,b) with the increasing pressure drop. Overall, we can see that seemingly the pressure drop affects the velocity profile but has weak relation to the particle volume fraction.

3.4. Influence of Carreau–Yasuda Parameters

Our next analysis examines the effect of the viscosity at infinite shear rate,  $[eqn]$ , which was varied from 100 to 700 Pa·s.

We can notice that the higher value of  $[eqn]$  leads to a slight increase near the center and slight decrease in particle volume fraction distribution ( $[eqn]$ ) near the wall in Figure 12a. The  $[eqn]$  and  $[eqn]$  relation curve moves downwards as shown in Figure 12b. The velocity decreases gradually as well as the average velocity with the increase in  $[eqn]$  (Figure 12c,d).

A strong impact is observed on the velocity profile, whereas the concentration profile is nearly unaffected. The same behavior is shown in Figure 13 when we increase the  $[eqn]$  from 1000 Pa·s to 7000 Pa·s. No considerable influence results from the change in the dimensionless parameter a as can be seen from Figure 14. In particular, the velocity profile decreases by increasing a, but stops changing regardless of the change in the parameter (Figure 14c), which is same for the average velocity.

Interestingly, the time constant  $[eqn]$  of the Carreau–Yasuda model affects the output of  $[eqn]$  and v curves to a high extent. When we vary the  $[eqn]$  parameter from 1.2 s to 5 s, in Figure 15a, the  $[eqn]$  distribution decreases near the centerline ( $[eqn]$ ) and increases near the wall ( $[eqn]$ ). However almost no change is observed for  $[eqn]$ .

The relation curve between  $[eqn]$  and  $[eqn]$  moves upward as shown in Figure 15 b, and the velocity profile as well as the average velocity linearly increase with the increasing  $[eqn]$  as demonstrated in Figure 13d and Figure 15c.

The parameter n was varied from 0.1 to 1, and the concentration profile is lower near the wall and is higher near the center as shown in Figure 16a. This means that if we assume the Newtonian fluid rheology, this would model the concentration distribution with higher values closer to the center and lower values closer to the wall. The velocity profile and average velocity decrease as n approaches 1 (Figure 16c,d).

3.5. Influence of Other Model Parameters

The fitting parameter,  $[eqn]$ , took the values from 0.4 to 0.95. This change affected the output considerably (Figure 17 and Figure 18). In particular, the particle volume fraction distribution is higher near the center and is lower near the wall. There is a point at about  $[eqn]$  where  $[eqn]$  does not change. The velocity profile and average velocity increase gradually as  $[eqn]$  increases.

The behavior observed in Figure 16a for small  $[eqn]$  and large  $[eqn]$  (red curve) is primarily numerical rather than physical. Equations (31) and (32) and the single-equation form (34) show that for high  $[eqn]$ , the particle concentration  $[eqn]$  can become very large near the y-axis (center of the tube), approaching the maximum packing fraction  $[eqn]$ .

Through constitutive relations such as Krieger–Dougherty, this leads to a rapid increase in the effective viscosity  $[eqn]$ , and the governing relation for the shear rate becomes extremely stiff. In addition, the cylindrical-coordinate terms (effectively containing  $[eqn]$ ) require careful regularization at  $[eqn]$ .

Numerically, the combination of (1)  $[eqn]$  approaching  $[eqn]$  (causing large viscosity), (2) steep gradients amplified by large  $[eqn]$ , and (3) singular  $[eqn]$ , can make the behavior at the axis to produce overflow/NaN or cause the integrator to reject steps near the axis, which appears as missing points for the red curve in the plotted results in Figure 17.

4. Conclusions

In this study, fully developed fluid flows with filled particles in a circular tube were analyzed. The Carreau–Yasuda and Krieger & Dougherty models were used to describe the particle-filled fluid and effective viscosities. respectively. The well-known diffusive flux model, proposed by Phillips et al. [6], was taken as a base of our analysis. The model was simplified by means of the non-dimensionalization technique, and a corresponding numerical method was applied for the solution. The verification was achieved by back-substituting the solution to the model that showed very small error. With the simplicity of our approach, the influence of particle and flow distributions by the model and the rheological parameters were observed. In particular, the particle concentration and velocity profiles depend strongly on  $[eqn]$  and n (Carreau–Yasuda parameters), the migration parameter  $[eqn]$ , and the average volume fraction  $[eqn]$ . The particle volume fraction distribution was weakly varied by the change in pressure gradient  $[eqn]$  and rheology parameters  $[eqn]$ ,  $[eqn]$  and a, but they strongly affected the velocity profile.

Future research can apply this method for determining particle volume fractions in other channels such as rectangular or curved channels. The limitation of our study is that our method is applied to only circular domains and can be extended to non-circular channels. This method could be particularly useful for studying similar fluid dynamics in more complex scenarios, such as coating and extrusion die manifolds, where the behavior of particle-filled fluids in confined geometries presents additional challenges and opportunities for analysis.