Numerical simulation of the effects of non-Newtonian property of the solvent on particle suspension rheology

To investigate the effects of the non-Newtonian properties of solvents on particle suspension rheology (non-Newtonian property and viscosity), we analyzed suspension flows with different solvent characteristics. We used the power-law model to express the non-Newtonian properties of solvents. The power-law index of the suspension was denoted as nsus. and that of the solvent as n sol . to distinguish between their non-Newtonian properties. Besides, to evaluate the non-Newtonian property of the particle suspension, we estimated its power-law index n sus . by comparing its particle velocity with the velocity profile of a power-law fluid. As a result, regardless of the different area fraction, the power-law index of the particle suspension was lower than that of the solvent for n sol . = 0 . 6 , 0.8, and 1.0. On the other hand, for n sol . = 0 . 4 and ϕ ≥ 1 . 57 , the power-law index of the particle suspension was higher than that of the solvent. Regarding the particle position, particles for n sol . = 0 . 4 were more concentrated on the wall side than those of the suspensions with n sol . = 0 . 6 , 0.8, and 1.0. Thus, this study confirms that a change in the non-Newtonian property of a suspension due to interactions between particles and power-law fluids is an essential factor when discussing the rheology of suspensions.


Introduction
Rheology deals with the deformation and flow of materials and is an interdisciplinary science that covers a wide range of disciplines, including physics, chemistry, engineering, medicine, biology, agriculture, pharmacy, food science, and domestic science.A dispersed system of solid particles in a liquid is called a suspension, and many suspensions are found in our daily lives (e.g.paint, ink, coating, cosmetics, and blood) because the rheology of suspensions (viscous properties) can be adjusted easily in response to many requirements by changing the content of suspended particles.The rheology of such suspensions is determined primarily by the interactions between particles and fluids.The most important factors are particle concentration (volume fraction), shape, particle-particle interaction, microstructure -which is generally defined by the relative position of the suspended particles -and solvent properties.Thus, the rheology of suspensions is complex, as it is determined by many factors.However, understanding these properties is important to solve various problems because it may improve the quality of industrial products.
Viscous fluids are classified into Newtonian and non-Newtonian fluids.Newtonian fluids exhibit constant viscosity independent of the shear rate, whereas the viscosity of non-Newtonian fluids depends on the shear rate.Shear-thinning fluids are characterized by decreasing apparent viscosity with increasing shear rate, and this phenomenon is called thixotropy.On the other hand, shear-thickening fluids have the opposite properties of shear-thinning fluids; they are characterized by increasing apparent viscosity with increasing shear rate.Such non-Newtonian fluids are frequently treated in the fields of chemical engineering and bioengineering.Blood is a non-Newtonian fluid, which is a suspension of blood cells (e.g.erythrocytes, leukocytes, and platelets) suspended in plasma.Erythrocytes are the majority of blood cell components, and their movement, deformation, and aggregation are thought to contribute to the non-Newtonian properties of blood.In a previous study, Mueller et al. 1 demonstrated that the dispersion of solid particles affected the non-Newtonian properties of the suspension even though the solvent was a Newtonian fluid.In addition, although plasma, the solvent of blood, is generally treated as a Newtonian fluid, plasma exhibits non-Newtonian properties depending on the content of plasma proteins. 1However, it is still unclear how the solid particles in the suspension affect the non-Newtonian characteristics of the suspension when a non-Newtonian fluid is used as the solvent of the particle suspension.
In addition, intrinsic viscosity has been recognized as an important factor when discussing the rheology of suspensions.In studies on suspension viscosity, many viscosity equations have been proposed, and Einstein presented Einstein's viscosity equation for dilute suspensions, which is given by equation (1) 2 : where h represents the relative viscosity, h 0 denotes the viscosity without particles, h eff denotes the effective viscosity (or apparent viscosity), and f denotes the solid volume fraction.h ½ denotes intrinsic viscosity, which is known to be h ½ = 2 for two-dimensional (2D) circular particles and h ½ = 2:5 for three-dimensional spherical particles.The applicable conditions for this formula are shown below; the suspension is dilute (low concentration), the particles are sufficiently small and homogeneously dispersed, the solvent of the suspension is Newtonian, and the inertial forces on particles are negligible.Therefore, when the solvent of the suspension is non-Newtonian, or when the particles flow, where inertial forces are not negligible, Einstein's viscosity equation may not be applicable.
Regarding the study on the inertial migration of particles, Segre´and Silberberg 3 showed that particles were concentrated at equilibrium positions in an inertial flow in a pipe.Fukui and Kawaguchi 4 considered particles' behavior in an inertial field and their effects on the microscopic particle arrangement and macroscopic suspension rheology.Hu et al. 5 also found that the solvent properties of a suspension affected the equilibrium position of particles.However, the trajectory and equilibrium position of particles in an inertial flow field are governed by the lift force of the particles, and the relationship between the non-Newtonian properties of the solvent and the lift force of the particles has not yet been clarified.
Regarding the effect of the non-Newtonian properties of the solvent on the viscosity of the suspension, Trofa and D'Avino 6 showed that power-law index affected intrinsic viscosity using fractal-shape aggregates in shear flow.Tanaka et al. 7 found that the non-Newtonian properties of the solvent affected the relative and intrinsic viscosity of the suspension in a flow field in which the particles are homogeneously distributed.Hu et al. 8 also found that the way particle train was formed depended on power-law index, Reynolds number and particle spacing.However, assuming an inertial flow field, the effect of the non-Newtonian property of the solvent on the viscosity of the suspension is unclear.
Based on the above background, the purpose of this study is to investigate the effects of the non-Newtonian property of the solvent on the rheology (non-Newtonian property and viscosity) of the suspension in a flow field with inertia.Then, we focused on relative viscosity and power-law index ratio n sus: =n sol: n sus: is power-law index of suspension, and n sol: is power-law index of solvent.In this study, we do not treat colloids.

Computational model
In this analysis, we introduced a two-way coupling scheme and performed a pressure-driven flow simulation in a 2D Poiseuille flow.Figure 1 depicts the numerical model used in this study.The initial positions of particles were randomly arranged, and the particles were non-homogeneously dispersed.The power-law model 9 was used to represent the non-Newtonian properties of the solvent, and the power-law indices n were set to 0.4-1.2.The suspension is composed of rigid circular buoyant particles in shear-thinning (n\1), Newtonian (n = 1:0), and shear-thickening (n.1) fluids.The virtual flux method 10,11 was applied to represent particles, the no-slip condition was applied to wall surfaces, and the periodic boundary condition 12 was applied to boundaries in the flow direction.Brownian motion was not considered in this analysis.The confinement C was defined as C = 2r=2l using the particle size 2r and channel width 2l, and numerical simulations were performed for C = 0:05, 0:10 and Reynolds number Re = 4, 32, 128.The number of particles N p was set to N p = 1-64, and the area fraction f varied from 0.20 to 12.6, accordingly.The computational domain was L x 3 L y = D 3 2l(0\x\D, À l\y\l).Table 1 shows the purpose and conditions of each analysis.In this analysis, the particle diameter was defined as D p , where D p = 30Dx in Case 1, D p = 20Dx in Case 2, and D p = 40Dx in Case 3, and we used a sufficient grid resolution that had little influence on the results of each analysis.To confirm the reliability of the program for Case 1, particles were ideally homogeneously distributed, and we compared our results with those obtained from Einstein's viscosity equation and investigated the effects of solvent characteristics on relative viscosity.These simulations were performed until non-dimensional time T = 10, 20, 100 correspond to Cases 1, 2, and 3, which represent sufficient time in terms of evaluation of particle suspension rheology.In this study, the particle-wall interaction was worked through fluid, and collision between a particle and a wall was not observed during the simulation.

Governing equation for fluids
In this study, incompressible fluid was used as the working fluid.The governing equations were the normalized lattice Boltzmann method 13 (2D9V model) with the incompressible formulation.The normalized lattice Boltzmann method maintains the algorithmic simplicity and high efficiency of the original lattice Boltzmann method while providing computational stability in higher Reynolds number regions and reducing the amount of memory used in the analysis.In the normalized lattice Boltzmann method, the pressure distribution function p a can be written in terms of the equilibrium part p eq a and non-equilibrium part p neq a of the pressure distribution function as in equation (2).
where e a denotes the discrete velocity vector of particles, l represents the relaxation time.p neq a and l are obtained using equations ( 3) and ( 4), respectively.
where d ij denotes Kronecker's delta, v a denotes the weight coefficient, P neq ij denotes the non-equilibrium part of the stress tensor, n represents the kinematic viscosity, and _ g represents the shear rate.Note that c is the advection velocity of particles, defined as c = dx=dt, using the time interval dt and the space interval dx.Furthermore, c s denotes the speed of sound, defined as c s = c= ffiffi ffi 3 p .

Power-law model
In this study, the power-law model 9 was used to represent the non-Newtonian property of the solvent.The  power-law model is frequently used in industry because of its ease of analytical treatment and its ability to represent a wide range of fluid flow properties.In the power-law model, apparent viscosity n _ g ð Þ is obtained using equation (5).
where m represents the power-law consistency, which can be derived from the Reynolds number Re, the characteristic velocity u 0 , the characteristic length D, and the power-law index n using the following equation: Equation ( 5) is calculated using m derived from equation (6).The shear rate _ g is obtained from the second invariant D II of the strain rate tensor.The second invariant is obtained from the value of the strain rate tensor by the following equation: The theoretical equation for the velocity profile in a 2D channel used in this study is obtained by the following equation: Governing equations for suspended particles The force acting on a particle is obtained by the following equation: where p represents the pressure on a small area of the particle, t represents the viscous stress on a small area of the particle, and Dl represents the small area of particles.The drag and lift coefficients, C L and C D , can be derived using F x and F y .
where r represents the density of the fluid, r represents the particle diameter.Because neutrally buoyant suspended particles were assumed in this analysis, the densities of fluid and particles were set equal.Suspended particles receive fluid force, so we used two-way-coupling scheme.Then we considered Newton's second law for translation and the angular equation of motion for rotation as follows: where F p is the external force vector, m p is particle mass, x p is the position vector, T p is torque, I is the moment inertia, u p is the particle angle.Then, equations ( 12) and ( 13) were solved using an explicit thirdorder Adams-Bashforth method: where u p is the velocity vector of suspended particles u p = u px , u py À Á .

Relative viscosity
The mass flow rate Q 0 of a Poiseuille flow with a welldeveloped boundary layer of a particle-free fluid can be obtained from the following equation: where h 0 denotes the viscosity of the particle-free fluid, and Dp 0 denotes the pressure difference for a channel length L without particles.The mass flow rate Q, in the case of a suspension containing particles, can be obtained similarly from the following equation: where h eff denotes the effective viscosity, and Dp denotes the pressure difference for a channel length L in the case of a suspension.Because Q = Q 0 in this study, the relative viscosity of the suspension, h, is obtained by the following equation:

Results and discussion
Validation: Homogeneous flow (Case 1) In this section, to validate the viscosity estimation using the pressure drop ratio expressed in equation ( 19) and the power-law fluid properties, we analyzed an ideal suspension flow with homogeneously dispersed particles.
As initial conditions, the velocity distribution and particle positions for each area fraction at T = 0 are depicted in Figure 2. To satisfy the condition for applying Einstein's viscosity equation (homogeneously dispersed, negligible inertia), particles were ideally placed at equal intervals, and Re = 4 was used in the analysis.
Figure 3 shows the relationship between relative viscosity and the area fraction for a Newtonian fluid (n = 1:0).The solid red line indicates the relative viscosity obtained using Einstein's viscosity equation, and the dotted black line denotes a straight line approximated using the least-squares method.The relative viscosity increased with increasing area fraction f.Furthermore, the slope of the approximated line yielded the intrinsic viscosity h ½ = 2:09, and the error with the theoretical value h ½ = 2:00 obtained from equation (1) was 4.5%.For more details, Table 2 shows the relative viscosity values obtained from the simulation and Einstein's viscosity equation and the difference between the two values.We found that for f ł 12:6, the error from the theoretical value obtained from Einstein's viscosity equation was less than 1.0%.In summary, for all area fractions, the difference between the relative viscosity obtained from the analysis and that from Einstein's viscosity equation was less than 1%, demonstrating the validity of our computational results.
Figure 4 shows the relationship between relative viscosity and the area fraction for each power-law fluid   (n = 0:8, 1:0, and 1:2).Blue circles indicate the relative viscosity for n = 0:8, black triangles indicate the relative viscosity for n = 1:0, and red rectangles indicate the relative viscosity for n = 1:2.Each solid line represents a straight line approximated using the least-squares method.The relative viscosity increased as the power-law index increased.This was consistent with the trend of Pal. 14 In addition, the intrinsic viscosity for each powerlaw index was estimated to be 1.57, 2.09, and 2.74 for the case of shear-thinning (n = 0:8), Newtonian (n = 1:0), and shear-thickening fluids (n = 1:2), respectively.Thus, intrinsic viscosity increased with increasing power-law index, which is a physically reasonable result.
Effect of solvent property on lift coefficient: Single particle (Case 2) In this section, to clarify the effect of the non-Newtonian property of the solvent on the microstructure -generally defined by the relative position of the suspended particles -of the suspension, we first focused on a single particle and investigated the lift force acting on the particle.A particle was allowed to flow downstream while keeping its y-position fixed, and the initial y-position varied from 0.0 to 0.6.Figure 5 depicts the relationship between the lift coefficient and normalized y-axis position.Note that the center of the channel and the channel wall correspond to y=l = 0:0 and y=l 6 1:0, respectively.The blue circles indicate the lift coefficient for n = 0:8, the black triangles indicate the lift coefficient for n = 1:0, and the red rectangles indicate the lift coefficient for n = 1:2.
The solid lines connecting the plots were obtained by cubic spline interpolations.The black line indicates the lift coefficient of zero, which represents the mechanical equilibrium position, where the wall repulsive force and the shear-velocity-dependent force are balanced. 15,16As illustrated in Figure 5, the y-axis position, where the lift coefficient becomes zero, moves toward the channel wall with decreasing power-law index.This was a similar trend to that of Hu et al. 5 For y=l.0:4, the difference in the lift coefficient increased when the particle approached the upper wall of the channel.This indicated that the wall repulsive force became stronger when the particle approached the upper wall of the channel, and the wall repulsive force was more dominant for the lift force compared with the shearvelocity-dependent force.We found that the power-law index had a significant effect on the wall repulsive force.
Next, focusing on the channel centerline around y=l = 0:1, we confirmed that the lift coefficient decreased with decreasing power-law index.We thought that this was because the power-law index affected the shear-velocity-dependent force.
To examine the effects of the power-law index on the lift force attributed to the shear rate, we focused on the velocity distribution of each power-law fluid.Figure 6 displays the theoretical values of the velocity distributions for each power-law index obtained from equation (9).Note that the center of the channel and channel wall correspond to y=l = 0:0 and y=l = 6 1:0, respectively.Figure 6 shows that the velocity distribution near the channel centerline was more blunted with decreasing power-law index.In other words, the shear rate decreased with decreasing power-law index near the channel centerline.Therefore, the relative velocity difference between the fluid and particle decreased with decreasing power-law index near the channel centerline, and it may reduce the lift coefficient.These facts suggest that the non-Newtonian property of the solvent of a suspension has a significant effect on the lift force acting on a particle.

Effect of solvent properties on the microstructure of suspensions: Non-homogeneous flow (Case 3)
To investigate the effects of the non-Newtonian property of the solvent on the microstructure, viscosity, and non-Newtonian properties of the suspension, we analyzed suspension flows with different solvent properties.Note that the initial positions of particles were randomly placed as in Figure 1.Generally, initial random positions affect particle flow patterns and the resultant particle arrangement, and these simulations were performed and repeated 10 times.
To evaluate the microstructure of the suspension, Figure 7 shows the PDF s, which represents the probability of particle existence, and the relationship between the particle positions in the y-direction after a sufficient time elapsed.The probability of particle existence represents the dispersion state of particles in the y-axis position.Note that the probability of particle existence was evaluated by dividing the y-direction into 20 parts.The solid black line represents the PDF when particles are homogeneously distributed.Figure 7 shows that the particles of the suspension with n = 0:4 more concentrated on the wall side than those of suspensions with n = 0:6 and 0:8. Figure 7 also shows that there are two distinct peaks for all power-law indices.These are because the equilibrium position of particles move closer to the wall surface as the power-law index decreases according to the previous study, 5 and the particles were concentrated at an equilibrium position due to the Segre´and Silberberg effect 3 when the inertia effects were not negligible.This trend was shown in another previous study using inertial multi particles flow in Newtonian fluid. 4Furthermore, Figure 5 shows lift coefficient near the center of the channel becomes lower as power-law index n decreases, and that near the wall becomes lower as power-law index n increases.Therefore, particles of suspension whose power-law index of solvent was low were hard to migrate from the channel centerline to the wall.
The heights of the two peaks decreased as the powerlaw index decreased.To investigate the cause of this phenomenon, we focused on the PDF value at the channel centerline.Figure 8 shows the relationship between the PDF and power-law index at the channel centerline.The PDF value at the channel centerline was derived from the average of the PDF values at y=l = 6 0:05.From Figure 8, we confirmed that the number of particles at the channel centerline increased as the power-law index decreased.This was caused by the decrease in the lift coefficient of particles with decreasing power-law index near the channel centerline, as illustrated in Figure 5. Therefore, the number of particles near the channel centerline increased with decreasing power-law index, whereas the number of particles concentrated at the equilibrium position (the position of the peak) on the upper and lower walls decreased.These findings suggested that the number of particles near the channel centerline increases as the power-law index decreases, with three possible peaks.
Regarding the effect of the particle position on the viscosity of the suspension, relative viscosity increases when the particle is close to the channel walls. 17,18rom Figure 5, we confirmed that the equilibrium position, where the lift coefficient becomes zero, shifted toward the channel walls as the power-law index decreased.Therefore, the particles may increase the relative viscosity when the power-law index decreases.On the other hand, for a homogeneous dispersion state (an equal distance between particles and walls for each  power-law index), relative viscosity decreases as the power-law index decreases, as depicted in Figure 4.
From the above, when the power-law index of the solvent decreases, the non-Newtonian property of the solvent decreases relative viscosity, whereas the aggregation position of the particles increases relative viscosity.In this analysis, we calculated 100 nondimensional time, and took the average of relative viscosity from 90 to 100 non-dimensional time.Figure 9 shows the relationship between relative viscosity and the area fraction.The data are the mean 6 1 SE.Note that the particles were non-homogeneously suspended because of inertial migration.From Figure 9, we confirmed that the relative viscosity of the suspension decreased as the power-law index decreased.Therefore, this suggests that the solvent properties had more significant effects on relative viscosity than the microstructure.In addition, the intrinsic viscosity for each power-law was estimated to be 0.63, 0.87, 1.02, and 1.21 for n = 0:4, 0:6, 0:8, and 1:0, respectively.Thus, intrinsic viscosity increased as the power-law index increased.
Figure 10 shows the time history of relative viscosity for each area fractions.The data are the mean 6 1 SE.In any area fractions, relative viscosity increases as time passes in shar thinning fluid.This tendency becomes stronger as power-law index n decreases, or the area fraction becomes higher.This is because inertial migration was caused and the number of particles near the wall increased as shown in Figure 7 as time passed.Thus, relative viscosity increased.Then, we consider the cases power-law index n is lower or the area fractions are higher.When power-law index n was lower, the number of particles near the wall increased as shown in Figure 7.In the case the area fractions were higher, the number of particles that contributed to the increase of relative viscosity increased.This phenomenon can be explained in the section ''Validation: Homogeneous flow (Case 1).''Therefore, we consider that the changes of the microstructure because of inertia resulted in the increase of relative viscosity.
When a non-Newtonian fluid is used as the solvent of a suspension, we investigated the effects of solid particles in the suspension on the non-Newtonian property of the whole suspension.Figure 11 shows the relationship between the power-law index ratio of the whole suspension n sus: and the solvent n sol: and the area fractions.The data are the mean 6 1 SE.The power-law index of the suspension is denoted as n sus: and that of the solvent as n sol: to distinguish between the non-Newtonian properties of the whole suspension and the solvent.Note that n sol: and the previously mentioned statement n were set equal.To evaluate the non-Newtonian property of the particle suspension, we estimated the power-law index n sus: of the particle suspension by comparing its particle velocity with the velocity profile of a power-law fluid.The value of n sus: was obtained using the least absolute value method such that the following cost function is minimized.
where u pk denotes the velocity of each particle in the flow direction, and u y pk À Á denotes the velocity obtained from equation (9).Note that the power-law index ratio n sus: =n sol: represents the change in the non-Newtonian property of the particle suspension due to the interactions between particles and power-law fluids.Figure 11 shows that for each power-law index, the value of n sus: =n sol: approaches 1.00 as the area fraction  decreases.This indicates that the properties of the whole suspension and the solvent become closer when the particle contribution decreases, which is physically reasonable.For n sol: = 0:6, 0:8, and 1.0, the power-law index ratio decreased as the area fraction increased.On the other hand, for n sol: = 0:4, the power-law index ratio increased as the area fraction increased.Table 3 shows the power-law index ratio n sus: =n sol: under each condition.It was confirmed that, for n sol: = 1:0, the power-law index ratio n sus: =n sol: was less than 1 regardless of area fractions, that is, the power-law index of the whole suspension was lower than that of the solvent.This is consistent with the trend of others. 1 Furthermore, for n sol: = 0:6 and 0:8, the power-law index ratio n sus: =n sol: was less than 1 regardless of area fractions, that is, the power-law index of the whole suspension was lower than that of the solvents, which is similar to the case in a Newtonian fluid (n sol: = 1:0).Therefore, the thixotropic property of the suspension flow was stronger.On the other hand, for n sol: = 0:4 and f ø 1:57, the power-law index ratio n sus: =n sol: exceeded 1, that is, the power-law index of the whole suspension was higher than that of the solvent, that is, the thixotropic properties of the suspension were weakened.Fukui and Kawaguchi 4 showed that for n sol: = 1:0, n sus: =n sol: increased as the Reynolds number decreased, and it exceeded 1 for Re = 2 and f = 4:07.This result was caused by the increase in the number of near-wall particles due to margination. 19In this analysis, as illustrated in Figure 7, particles were more concentrated on the wall side for n sol: = 0:4, and it decreased the shear rate near the wall, which results in n sus: =n sol: .1 for f ø 1:57 and n sol: = 0:4.Therefore, we discovered that there may be a critical point on the particle y-axis position that has a significant effect on the shear rate near the wall.

Conclusion
In this study, we assumed an inertial flow field and investigated the effects of the non-Newtonian properties of solvents on the rheology (non-Newtonian property and viscosity) of suspensions with nonhomogeneously dispersed solid particles.Based on the analysis results of the microstructures of the suspensions, we confirmed that the particles of the suspension with n sol: = 0:4 were more concentrated on the wall   side than those of the suspensions with n sol: = 0:6, 0:8, and 1:0.We also found that the number of particles near the channel centerline increases as the power-law index decreases and that three peaks may emerge.Furthermore, we confirmed that the relative viscosity of the suspension decreased as the power-law index decreased.This suggests that the solvent properties have more significant effects on relative viscosity than the microstructure.Regarding the non-Newtonian properties of the suspensions, for n sol: = 0:6 and 0:8, the power-law index ratio n sus: =n sol: was less than 1 regardless of area fractions, which was similar to the case of a Newtonian fluid (n sol: = 1:0).On the other hand, for n sol: = 0:4 and f ø 1:57, the power-law index ratio n sus: =n sol: exceeded 1, that is, the power-law index of the whole suspension was higher than that of the solvent.Thus, this study confirms that a change in the non-Newtonian property of a suspension due to interactions between particles and power-law fluids is an essential factor when discussing the rheology of suspensions.power-law index of suspension n sol:

Appendix
power-law index of solvent u pk velocity of each particle in the flow direction

Figure 1 .
Figure 1.Schematic view of computational model.The rigid buoyant circular particles are suspended in a pressure-driven flow through a 2D channel.The particles are randomly distributed in the channel.

Figure 2 .
Figure 2. Normalized axial velocity distribution and particle position for different area fractions f at nondimensional time T = 0.The rigid buoyant circular particles were suspended in a pressure-driven flow through a 2D channel.Note that the particles were homogeneously distributed in the channel to satisfy the assumption in Einstein's formula: (a) f = 0.79%, (b) f = 3.14%, (c) f = 7.07%, and (d) f = 12.7%.

Table 1 .
Purpose and conditions of the simulation.

Table 2 .
Comparison of relative viscosity.Relative viscosity obtained from our simulation h, that obtained from Einstein's formula h th and the difference between these two values.

Table 3 .
Power-law index ratio of the whole suspension n sus: and the solvent n sol: .