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

1. Introduction

Fluid flows with suspended particles are found in many industrial applications such as concrete pumping [1] or extrusion and injection molding processes with filled polymers [2] [3], thus understanding the fundamentals of these flows through numerical modeling is essential. It is more important to investigate the behavior of suspension flows with generalized rheological models to have a broader range of application. We focus on steady state flows in circular tube by assuming the non-Newtonian fluids and present the simple approach of numerical modeling the particle volume fraction and velocity distributions by using the simplified form of shear-induced diffusion equation. We then present the influence of physical parameters as well as the rheological parameters of Carreau-Yasuda model on the concentration and velocity profiles. Understanding the rheology of electrode slurries is one of the applications crucial as it profoundly influences the final coating microstructure. Issues such as high slurry viscosity, which increases pressure and limits coating speed, elasticity leading to instabilities and defects, and excessive flow causing slumping and poorly structured coatings, highlight the complexities in optimizing these formulations. Due to the diverse solvent systems and components involved, pinpointing the origins of these adverse rheological properties can be challenging [4]. Additionally, precise measurement of extensional properties using miniature rheometers has revealed significant variability in slurry behavior across different formulations and binder systems [5]. The paper by Phillips et al. [6] is the fundamental reference used in this paper. They derived the shear-induced diffusion equation for the particle volume fraction ($\varphi$) and present the corresponding no-flux boundary condition. The difference between our work and theirs is that they assume Newtonian fluid in Couette and Poiseuille flows, whereas, we consider non-Newtonian flows governed by the Carreau-Yasuda [7] model. In 2003, Chen et al. [8] studied the shear-induced particle migration in non-Newtonian flows of nickel-powder-filled polymers. They facilitated the power- law model to describe the viscosity and shear rate relation of the flow. In addition, they use the Arrhenius law to add temperature dependence to the viscosity model. It is nontrivial to solve the coupled differential equations of velocity and particle volume fraction distributions. Chen et al. [9] presented the numerical scheme to get the solutions of particle volume fraction and velocity fields for non-Newtonian flows for both steady and unsteady cases. They used the Newton-Raphson iterative method to solve the velocity fields and the finite difference method with semi-implicit scheme (known as Crank-Nicolson [10]) to solve the diffusive equation. In contrast with their paper, we use non-dimensionalization of the equations and some changes of variables techniques presented in Wang [11] to reduce the system of differential equations into a system of algebraic equations. The solution is obtained by solving numerically the nonlinear algebraic system which is computationally less intensive than the differential equations-based methods. Moreover, the Carreau-Yasuda model is more general as compared with the Power-Law used. Several other papers have been studied the different fluid flows in various processes considering the particle migration. Choi et al. [12] studies the pipe flow of pumped concrete assuming the shear-induced particle migration. They used the Bingham fluid model to describe the viscosity - shear rate relation. The commercial CFD Fluent software was used to get the numerical solutions. The paper by Siqueira and Carvalho [13] presented a numerical investigation of particle migration in a planar extrudate flow of suspensions of some hard spheres. In particular, they studied the Newtonian flow in parallel-plates and then the extruded flow in the ambient atmosphere. The numerical solutions have been obtained by the stabilized finite element method together with the elliptic mesh generation method to compute the position of the free surface. Rebouças et al. [14] studied the fully-developed suspension flow in a circular tube and the effects of particle and shear-induced viscosity. Their work is similar to ours in the sense that they also presented the numerical solution of shear-induced diffusion equation with the assumption of non-Newtonian flow and investigated the effects of average particle volume fraction on the concentration and velocity profiles. Moreover, the results of both concentration and velocity profiles for different bulk concentration ($\overline{\varphi }$) of the suspensions and the parameter $K_c/K_{\eta }$ agree with our numerical solutions. The difference between their approach and ours is that they used the procedure to overcome the singularity of the model near the centerline of the tube. Also they used Cross model to describe the fluid viscosity whereas we use the Carreau-Yasuda model that has more parameters (with extra variable a). Another difference is that we present the effects of the whole range of parameters including the pressure gradient, average particle volume fraction, all 5 parameters of Carreau-Yasuda and the model parameter $K_c/K_{\eta }$ (thereafter $\alpha$). In addition to numerical studies, Rueda [15] presented a complete experimental study of the rheological characteristics of suspension flows. In particular, the author evaluates the influence of various parameters on the viscosity.

2. Model and Methods

2.1. Governing Equations of Fluid

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

$$\nabla ·\mathbf{v}=0$$

(1)

the conservation equation of linear momentum

$$\nabla ·\tau +\nabla p=0$$

(2)

and the Newton’s Law of viscosity

$$\tau =\eta \left(\dot{\gamma }\right)\dot{\gamma }$$

(3)

2.2. Particle Poiseuille Flow Through a Cylindrical Tube

Let $\varphi (r,t)$ denote the particle fraction. Consider the Krieger & Dougherty [16] model with the following terms on particle volume induced relative viscosity

$$\begin{array}{c}\hfill \begin{array}{cccc}\eta (\varphi ,\dot{\gamma })\hfill & =\hfill & {\eta }_r\left(\varphi \right){\eta }_s\left(\dot{\gamma }\right)\hfill & -\mathrm{viscosity}\mathrm{of}\mathrm{concentrated}\mathrm{suspension}\\ {\eta }_r\left(\varphi \right)\hfill & =\hfill & {(1-\varphi /{\varphi }_m)}^{-1.82}\hfill & \mathrm{relative}\mathrm{viscosity}\\ {\eta }_s\left(\dot{\gamma }\right)\hfill & =\hfill & {\eta }_{\infty }+({\eta }_0-{\eta }_{\infty }){(1+{\left(\lambda \dot{\gamma }\right)}^a)}^{{\displaystyle \frac{n-1}{a}}}\hfill & -\mathrm{viscosity}\mathrm{of}\mathrm{the}\mathrm{media}/\mathrm{fluid}\\ {\mathbf{N}}_c\hfill & =\hfill & -K_c{a}^2\varphi \nabla \left(\dot{\gamma }\varphi \right)\hfill & -\mathrm{flux}\mathrm{by}\mathrm{varying}\mathrm{collision}\\ {\mathbf{N}}_{\eta }\hfill & =\hfill & -K_{\eta }{a}^2\dot{\gamma }{\varphi }^2{\displaystyle \frac{\nabla \eta }{\eta }}\hfill & -\mathrm{flux}\mathrm{by}\mathrm{spatially}\mathrm{varying}\mathrm{viscosity}\\ {\mathbf{N}}_b\hfill & =\hfill & -D\nabla \varphi \hfill & -\mathrm{accounts}\mathrm{for}\mathrm{Brownian}\mathrm{diffusion}\mathrm{of}\mathrm{particles}\end{array}\end{array}$$

(4)

and convective diffusion of particles satisfying the conservation equation for the particle volume fraction

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \varphi }{\partial t}}=-\nabla ({\mathbf{N}}_c+{\mathbf{N}}_{\eta }+{\mathbf{N}}_b)\end{array}$$

(5)

According to Phillips et al. [6], ${\mathbf{N}}_b$ can be neglected for very large Peclet number $\mathrm{Pe}\equiv \dot{\gamma }{a}^2/D$. Then replacing the flux terms in (5) with the equations from (4) gives us

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \varphi }{\partial t}}-a^2\nabla \left[K_c\varphi \nabla \left(\dot{\gamma }\varphi \right)+K_{\eta }\dot{\gamma }{\varphi }^2{\displaystyle \frac{\nabla \eta }{\eta }}\right]=0\end{array}$$

(6)

Since we are considering the flow in cylindrical tube, the divergence is $\nabla =\left({\displaystyle \frac{\partial }{\partial r}}+{\displaystyle \frac{\partial }{\partial z}}\right)$ in general form. No-flux boundary condition at the wall which is represented as

$$\begin{array}{c}\hfill \mathbf{n}·\left[K_c\varphi \nabla \left(\dot{\gamma }\varphi \right)+K_{\eta }\dot{\gamma }{\varphi }^2{\displaystyle \frac{\nabla \eta }{\eta }}\right]=0\end{array}$$

(7)

So, combining (6) and (7), we get the generalized system of differential equations of velocity profile and particle volume fraction distribution, which is

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\displaystyle \frac{\partial p}{\partial z}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r\dot{\gamma }\eta \right)& =& 0\hfill \\ \hfill \\ \hfill {\displaystyle \frac{\partial \varphi }{\partial t}}-a^2\nabla \left[K_c\varphi \nabla \left(\dot{\gamma }\varphi \right)+K_{\eta }\dot{\gamma }{\varphi }^2{\displaystyle \frac{\nabla \eta }{\eta }}\right]& =& 0\hfill \end{array}\end{array}$$

(8)

with the no-flux boundary condition at the wall

$$\begin{array}{c}\hfill \mathbf{n}·\left[K_c\varphi \nabla \left(\dot{\gamma }\varphi \right)+K_{\eta }\dot{\gamma }{\varphi }^2{\displaystyle \frac{\nabla \eta }{\eta }}\right]=0\end{array}$$

(9)

2.3. Steady State Velocity and Particle Distributions

In steady state flow, ${\displaystyle \frac{\partial \varphi }{\partial t}}=0$ and we assume the r-direction, so the second equation in the coupled system (8) becomes

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial }{\partial r}}\left[K_c\varphi {\displaystyle \frac{\partial }{\partial r}}\left(\dot{\gamma }\varphi \right)+K_{\eta }\dot{\gamma }{\varphi }^2{\displaystyle \frac{1}{\eta }}{\displaystyle \frac{\partial \eta }{\partial r}}\right]=0\end{array}$$

(10)

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

$$\begin{array}{c}\hfill K_c\varphi {\displaystyle \frac{\partial }{\partial r}}\left(\dot{\gamma }\varphi \right)+K_{\eta }\dot{\gamma }{\varphi }^2{\displaystyle \frac{1}{\eta }}{\displaystyle \frac{\partial \eta }{\partial r}}=0\end{array}$$

(11)

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

$$\begin{array}{c}\hfill {\displaystyle \frac{\dot{\gamma }\varphi }{{\dot{\gamma }}_w{\varphi }_w}}-{\left({\displaystyle \frac{{\eta }_w}{\eta }}\right)}^{K_{\eta }/K_c}=0\end{array}$$

(12)

Now we want to simplify the system. By rearrangement of the first part of the coupled system (8) we get

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial }{\partial r}}\left(r\dot{\gamma }\eta \right)=-{\displaystyle \frac{\partial p}{\partial z}}r\end{array}$$

(13)

and integration of (13) gives us

$$\begin{array}{c}\hfill \dot{\gamma }\eta =-{\displaystyle \frac{\partial p}{\partial z}}{\displaystyle \frac{r}{2}}+C\end{array}$$

(14)

applying the boundary case (at $r=0$), we get the constant of integration $C=0$. So, we have the simplified equation

$$\begin{array}{c}\hfill \dot{\gamma }\eta =-{\displaystyle \frac{\partial p}{\partial z}}{\displaystyle \frac{r}{2}}\end{array}$$

(15)

Summarizing the coupled ODE for steady-state case, we have the following system of differential equations

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \dot{\gamma }\eta +{\displaystyle \frac{\partial p}{\partial z}}{\displaystyle \frac{r}{2}}& =& 0\hfill \\ \hfill \\ \hfill {\displaystyle \frac{\dot{\gamma }\varphi }{{\dot{\gamma }}_w{\varphi }_w}}-{\left({\displaystyle \frac{{\eta }_w}{\eta }}\right)}^{K_{\eta }/K_c}& =& 0\hfill \end{array}\end{array}$$

(16)

2.4. Nondimensionalization and Numerical Details

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

$$\begin{array}{c}\hfill \begin{array}{cccc}\hfill \tilde{r}& =& r/R& \in [0,1]\hfill \\ \hfill s& =& \dot{\gamma }/{\dot{\gamma }}_w& \in [0,1]\hfill \\ \hfill \tilde{\eta }& =& \eta /{\eta }_w& \in [1,+\infty ]\hfill \end{array}\end{array}$$

(17)

Substituting the new variables (17) into (16) gives us the system in nondimensional form as

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\eta }_w{\dot{\gamma }}_{w}s\tilde{\eta }+{\displaystyle \frac{dp}{dz}}{\displaystyle \frac{\tilde{r}R}{2}}& =& 0\hfill \\ \hfill \\ \hfill s{\displaystyle \frac{\varphi }{{\varphi }_w}}-{\left({\displaystyle \frac{1}{\tilde{\eta }}}\right)}^{K_{\eta }/K_c}& =& 0\hfill \end{array}\end{array}$$

(18)

At the wall we have $s=1$, $\tilde{\eta }=1$, and $\tilde{r}=1$, and the first equation in (18) becomes,

$${\eta }_w{\dot{\gamma }}_w=-{\displaystyle \frac{dp}{dz}}{\displaystyle \frac{R}{2}}=\left({\displaystyle \frac{P_0-P_L}{2L}}\right)R$$

(19)

where ${\displaystyle \frac{dp}{dz}}=(P_0-P_L)/L>0$ is the magnitude of the pressure gradient. Suppose we use the Krieger-Dougherty equation for ${\eta }_r\left(\varphi \right)$ and the Carreau-Yasuda model for ${\eta }_s\left(\dot{\gamma }\right)$ as in (4), we have

$${\eta }_w={\eta }_w({\dot{\gamma }}_w,{\varphi }_w)={\eta }_r\left({\varphi }_w\right){\eta }_s\left({\dot{\gamma }}_w\right)$$

(20)

Combining Eqs. 19 and 20 gives

$${\eta }_r\left({\varphi }_w\right){\eta }_s\left({\dot{\gamma }}_w\right){\dot{\gamma }}_w=\left({\displaystyle \frac{P_0-P_L}{2L}}\right)R$$

(21)

Here, ${\varphi }_m$ is taken to be $0.68$. Given ${\varphi }_w$, $(P_0-P_L)/L$, and R, we can obtain ${\dot{\gamma }}_w$ by solving the above nonlinear equation (Eq. 21) for ${\dot{\gamma }}_w$. Once ${\dot{\gamma }}_w$ is obtained, then it is straightforward to obtain ${\eta }_w$, which is used as the unit for viscosity. Using the identity shown in Eq. 19, we can rewrite Eq. 18 as

$$\tilde{r}=s\tilde{\eta }(s,\varphi )$$

(22)

$$\tilde{\eta }={\left(s\varphi /{\varphi }_w\right)}^{-\alpha }$$

(23)

where $\alpha =K_c/K_{\eta }$ and is taken to be $0.66$ in this work, following Phillips et al. [6]. Again, using the Krieger-Dougherty equation for ${\eta }_r\left(\varphi \right)$ and the Carreau-Yasuda model for ${\eta }_s\left(\dot{\gamma }\right)$, we have

$$\tilde{\eta }=\eta /{\eta }_w={\eta }_r\left(\varphi \right){\eta }_s\left(\dot{\gamma }\right)/{\eta }_w={\eta }_r\left(\varphi \right){\tilde{\eta }}_s\left(s\right)$$

(24)

where

$${\tilde{\eta }}_s\left(s\right)={\tilde{\eta }}_{\infty }+\left({\tilde{\eta }}_0-{\tilde{\eta }}_{\infty }\right){\left[1+{\left(\tilde{\lambda }s\right)}^a\right]}^{(n-1)/a}$$

(25)

Here, ${\tilde{\eta }}_0={\eta }_0/{\eta }_w$, ${\tilde{\eta }}_{\infty }={\eta }_{\infty }/{\eta }_w$, and $\tilde{\lambda }=\lambda {\dot{\gamma }}_w$. The $\tilde{\lambda }$ term is the time constant $\lambda$ in dimensionless form, using ${\dot{\gamma }}_w^{-1}$ as the characteristic time scale. Combining Eqs. (22) and (23), we have

$$f_1(\tilde{r},s,\varphi )=\tilde{r}-s^{1-\alpha }{\left(\varphi /{\varphi }_w\right)}^{-\alpha }$$

(26)

and combining Eqs. (22)-(24), we have

$$f_2(\tilde{r},s,\varphi )=\tilde{r}-s{\eta }_r\left(\varphi \right){\tilde{\eta }}_s\left(s\right)$$

(27)

From the above two equations, Eqs.(26) and (27), one can easily obtain the volume fraction $\varphi \left(\tilde{r}\right)$, the shear rate distribution $s\left(\tilde{r}\right)$, and the volume fraction as a function of the shear rate, $\varphi \left(s\right)$. The velocity distribution can be obtained using the method of Wang (2021, under review). Using $R{\dot{\gamma }}_w$ as the unit for velocity, i.e., $\tilde{u}=v/\left(R{\dot{\gamma }}_w\right)$, we have

$${\displaystyle \frac{d\tilde{u}}{d\tilde{r}}}=-s$$

(28)

Once the shear rate distribution $s\left(\tilde{r}\right)$ is obtained, it is straightforward to obtain $\tilde{u}\left(r\right)$ by solving the above ODE numerically.

If we try to solve the system for s, complications in the numerical computations arise, for instance, ${\tilde{\eta }}_s\left(s\right)$ is highly nonlinear. Thus, following the idea from the paper of Wang (2021), it is wise to solve for $\tilde{r}$ instead of s. This means, given $s=[0,\dots ,1]$, we want to solve for $\tilde{r}$ and $\varphi$ which can be easily obtained by root finding ($\varphi$) the following equation

$$s^{1-\alpha }{\left(\varphi /{\varphi }_w\right)}^{-\alpha }-s{\eta }_r\left(\varphi \right){\tilde{\eta }}_s\left(s\right)=0$$

(29)

and substituting the solution back to one of the equations (26) or (27) to find $\tilde{r}$. Because ${\varphi }_w$ is required to solve the system, with given $\overline{\varphi }$, it can be estimated by using the concept of binary search or so-called "interval halving method". Since ${\varphi }_w$ is always less than the $\overline{\varphi }$, we know that the true ${\varphi }_w$ is greater than 0 and less than $\overline{\varphi }$. By using this information, we can find the the corresponding ${\varphi }_w$ that is shown as an intermediate step in Figure 1.

3. Results and Discussion

3.1. Verification and Validation of Numerical Results

Now we want to generate some output by using the algorithm shown in Figure 1. 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 2 and Figure 3 show the numerical solutions of the system.

The solutions of particle volume fraction ($\varphi$) and the dimensionless shear rate (s) distributions are shown in Figure 2a,b. The dimensionless velocity ($\tilde{u}$) profile in Figure 3a is computed from the integration of the simple ODE (28). The iterative behavior of the inner algorithm that finds the best approximation of ${\varphi }_w$ is shown in Figure 3b.

To verify the numerical model and the corresponding Python code used to solve for $\varphi$ and s, we back-substitute these solutions into the system (Eqs. 26, 27 labeled as $f_1$ and $f_2$, respectively) and check the errors. As we can notice, the errors are considerably small with the order of ${10}^{-16}$ (Figure 3).

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. Authors of the analytical solution of the equation governing $\varphi$ at steady state presented the explicit formula for $K_c/K_{\eta }=0.66$ as

$$\varphi \left(\tilde{r}\right)={\displaystyle \frac{{\varphi }_m}{1+i\tilde{r}}}$$

(30)

where $i={\displaystyle \frac{{\varphi }_m-{\varphi }_w}{{\varphi }_w}}$. This explicit formula describes the particle volume fraction distribution at steady state Newtonian case. We compare our numerical algorithm with the analytical solution for different $\overline{\varphi }$ and the result is shown in Figure 4.

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 power-law model as following

$$\eta (\dot{\gamma },\varphi )={(1-{\displaystyle \frac{\varphi }{{\varphi }_m}})}^{-1.82}K{\dot{\gamma }}^{n-1}$$

(31)

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 Eq. 15 we have

$$\dot{\gamma }={\left[{\displaystyle \frac{dp}{dz}}{\displaystyle \frac{r}{2K{\eta }_r\left(\varphi \right)}}\right]}^{1/n}$$

(32)

and substituting this formula into Eq. 12 gives us

$${\displaystyle \frac{\varphi }{{\varphi }_w}}\dot{\gamma }={\displaystyle \frac{{\dot{{\gamma }_w}}^{1+(n-1)K_{\eta }/K_c}}{{\dot{\gamma }}^{(n-1)K_{\eta }/K_c}}}{\left[{\displaystyle \frac{{\eta }_r\left({\varphi }_w\right)}{{\eta }_r\left(\varphi \right)}}\right]}^{K_{\eta }/K_c}$$

(33)

$${\displaystyle \frac{\varphi }{{\varphi }_w}}{\dot{\gamma }}^{1+(n-1)K_{\eta }/K_c}={\dot{{\gamma }_w}}^{1+(n-1)K_{\eta }/K_c}{\left[{\displaystyle \frac{{\eta }_r\left({\varphi }_w\right)}{{\eta }_r\left(\varphi \right)}}\right]}^{K_{\eta }/K_c}$$

(34)

$${\displaystyle \frac{\varphi }{{\varphi }_w}}{\displaystyle \frac{{\dot{\gamma }}^{1+(n-1)K_{\eta }/K_c}}{{\dot{{\gamma }_w}}^{1+(n-1)K_{\eta }/K_c}}}={\left[{\displaystyle \frac{{\eta }_r\left({\varphi }_w\right)}{{\eta }_r\left(\varphi \right)}}\right]}^{K_{\eta }/K_c}$$

(35)

$${\left[{\displaystyle \frac{{\eta }_r\left({\varphi }_w\right)}{{\eta }_r\left(\varphi \right)}}\right]}^{{\displaystyle \frac{1}{n}}(1-K_{\eta }/K_c)}={\displaystyle \frac{{\varphi }_w}{\varphi }}{\displaystyle \frac{1}{{\tilde{r}}^{\frac{1}{n}+\frac{n-1}{n}{K}_{\eta }/K_c}}}$$

(36)

$${\left[{\displaystyle \frac{{\varphi }_m-{\varphi }_w}{{\varphi }_m-\varphi }}\right]}^{{\displaystyle \frac{1}{n}}(-1.82)(1-K_{\eta }/K_c)}={\displaystyle \frac{{\varphi }_w}{\varphi }}{\displaystyle \frac{1}{{\tilde{r}}^{\frac{1}{n}+\frac{n-1}{n}{K}_{\eta }/K_c}}}$$

(37)

Letting $-1.82(1-K_{\eta }/K_c)=1$ with $K_{\eta }/K_c\approx 1.549$ yields the analytical solution for the fully-developed steady-state power-law fluid as

$$\varphi ={\displaystyle \frac{{\varphi }_w}{{\tilde{r}}^{\frac{1}{n}(1.549n-0.549)}}}{\left({\displaystyle \frac{{\varphi }_m-\varphi }{{\varphi }_m-{\varphi }_w}}\right)}^{1/n}$$

(38)

Figure 5 shows both analytical solution from the above derivation and our numerical solution for different power-law index n. Interestingly, when $(1.549n-0.549)=0$, that is $n\approx 0.3544$, $\tilde{r}$ vanishes and the solution of Eq. 38 is constant for any value of $\tilde{r}$. For values $n⪆0.3544$ the concentration profiles are higher near the centerline and lower near the wall of the circular tube, however if $n⪅0.3544$ we observe very small particle volume fraction near the centerline and it increases towards the wall. From this special case, we can conclude that for $K_{\eta }/K_c\approx 1.549$ and $n\approx 0.3544$ the particle volume fraction distribution is exactly flat at $\varphi =0.40$. Generally speaking, when $(1+(n-1)K_{\eta }/K_c)=0$ in Eq. 38 the concentration profile is $\varphi =\overline{\varphi }$ along the radial position.

3.1.3. Other Relative Viscosity Models

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

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 ($\overline{\varphi }$), pressure drop (${\displaystyle \frac{dp}{dz}}$), the Carreau-Yasuda model parameters, and the $K_c/K_{\eta }$. 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 7a, the concentration profile moves upward as $\overline{\varphi }$ increases. This observation is as expected as the result in Figure 7b where the relation between ${\varphi }_w$ and $\overline{\varphi }$ lies below the diagonal line which agrees with the theory that $0<{\varphi }_w<\overline{\varphi }$. We can also see that the velocity profile decreases with the increasing average particle volume fraction $\overline{\varphi }$ (Figure 7c). Correspondingly, the average velocity decreases as well (Figure 7d). In general, the plots in Figure 7 suggest that the particle volume fraction and velocity distribution is highly dependant on the average particle volume fraction.

3.3. Influence of Pressure Gradient

Next, we want to see the influence of pressure drop on the output, so we vary it from $0.2\times {10}^6$ Pa to $1.4\times {10}^6$ Pa, and keep everything else constant (as in Table 1).

The behavior of the concentration profile is not uniform, since, in Figure 8a, it decreases near the centerline ($\tilde{r}=0$) but increases near the wall ($\tilde{r}=1$) as the pressure drop increases. Clearly, the ${\varphi }_w$ vs. $\overline{\varphi }$ curve moves upwards as shown in Figure 8b, and the velocity profile increases as well as the average velocity (Figure 9a,b) with the increasing pressure drop. Overall, we can see that seemingly the pressure drop affects the velocity profile, but has weak relation with the particle volume fraction.

3.4. Influence of Carreau-Yasuda Parameters

Our next analysis examines the effect of the viscosity at infinite shear rate, ${\eta }_{\infty }$, which was varied from 100 to 700 Pa·s.

We can notice the higher value of ${\eta }_{\infty }$ leads to a slight increase near the center and slight decrease of particle volume fraction distribution ($\varphi$) near the wall in Figure 10a. The ${\varphi }_w$ and $\overline{\varphi }$ relation curve moves downwards as shown in Figure 10b. The velocity decreases gradually as well as the average velocity with the increase of ${\eta }_{\infty }$ (Figure 10c,d).

Somehow strong impact is observed on the velocity profile, but the concentration profile stays almost unaffected. Same behavior is shown in Figure 11 when we increase the ${\eta }_0$ from 1000 Pa·s to 7000 Pa·s. No considerable influence is made by the change of dimensionless parameter a as can be seen from Figure 12. In particular, the velocity profile decreases by increasing a, but stops changing regardless the change in the parameter (Figure 12c) which is same for the average velocity.

Interestingly, the time constant $\lambda$ of Carreau-Yasuda model affects the output of $\varphi$ and v curves to a high extent. When we vary the $\lambda$ parameter from 1.2 s to 5 s, in Figure 13a, the $\varphi$ distribution decreases near the centerline ($\tilde{r}=0$) and increases near the wall ($\tilde{r}=1$). However almost no change is observed about $\tilde{r}=0.6$.

The relation curve between ${\varphi }_w$ and $\overline{\varphi }$ moves upward as shown in Figure 13b, and the velocity profile as well as the average velocity linearly increase with the increasing $\lambda$ as demonstrated in Figure 13c,d).

The parameter n was varied from 0.1 to 1 and the concentration profile is less near the wall and is more near the center as shown in Figure 14a. 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 the wall. The velocity profile and average velocity decrease as n goes towards 1 (Figure 14c,d).

3.5. Influence of Other Model Parameters

The fitting parameter, $\alpha =K_c/K_{\eta }$, took the values from 0.4 to 0.95. This change affected the output considerably (Figure 15 and Figure 16). In particular, the particle volume fraction distribution is higher near the center and is lower near the wall. There is a point about $\tilde{r}=0.6$ where $\varphi$ does not change. The velocity profile and average velocity increase gradually as $\alpha$ increases.

4. Conclusions

In this study, fully-developed fluid flows with filled particles in 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 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 simplicity of our approach, the influence of particle and flow distributions by the model and the rheological parameters were observed. In particular, both particle concentration and velocity profiles showed strong dependancy on the relation with the rheology parameters $\lambda$, n and model parameter $\alpha$ and $\overline{\varphi }$. The particle volume fraction distribution was weakly varied by the change in pressure gradient ${\displaystyle \frac{dp}{dz}}$ and rheology parameters ${\eta }_{\infty }$, ${\eta }_0$ and a, but they strongly affected the velocity profile.

For future research, the circuit network method [19] offers a promising avenue for further investigation. 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.