World Journal of Mechanics
Vol. 1  No. 5 (2011) , Article ID: 8161 , 12 pages DOI:10.4236/wjm.2011.15029

Numerical Simulation of Non-Newtonian Pseudo-Plastic Fluid in a Micro-Channel Using the Lattice Boltzmann Method

Hossein Hamedi, Mohammad Hassan Rahimian

Department of Mechanical Engineering, University of Tehran, Tehran, Iran


Received July 24, 2011; revised August 29, 2011; accepted September 12, 2011

Keywords: Lattice Boltzmann method, Non-Newtonian fluid, Power-Law model


In this paper, the power-law model for a non-Newtonian (pseudo-plastic) flow is investigated numerically. The D2Q9 model of Lattice Boltzmann method is used to simulate the micro-channel flow with expansion geometries. This geometry is made by two squared or trapezoid cavities at the bottom and top of the channel which can simulate an artery with local expansion. The cavities are displaced along the channel and the effects of the displacements are investigated for inline structures and staggered ones (anti-symmetric expansion). The method is validated by a Poiseuille flow of the power-law fluid in a duct. Validation is performed for two cases: The Newtonian fluid and the shear thinning fluid (pseudo-plastic) with n = 0.5. The results are discussed in four parts: 1) Pressure drop; it is shown that the pressure drop along the channel for inline cavities is much more than the pressure drop along the staggered structures. 2) Velocity profiles; the velocity profiles are sketched at the centerline of the cavities. The effects of pseudo-plasticity are discussed. 3) Shear stress distribution; the shear stress is computed and shown in the domain. The Newtonian and non-Newtonian fluids are discussed and the effect of the power n on shear stress is argued. 4) Generated vortices in the cavities are also presented. The shape of the vortices is depicted for various cases. The results for these cases are talked over and it is found that the vortices will be removed for flows with n smaller than 0.5.

1. Introduction

During the last two decades great attention has been paid to the Lattice Boltzmann method (LB). The LB method has a remarkable ability to solve single phase, multiphase, single component, and multi-component problems in complex geometries. Also the LB method is a strong method in simulating complex fluid flows such as non-Newtonian Fluids.

Non-Newtonian fluid plays an important role in various industries. The flow of petrochemical materials in tubes and blood in the arteries are a simple example of these problems. Blood is a type of a non-Newtonian fluid that should be simulated with a suitable model. Power law shear-thinning (pseudo-plastic) fluid is a well presented model for simulating blood in the arteries. Shearthinning fluid is a type of a generalized non-Newtonian fluid. In this kind of fluid the shear stress is a function of a second invariant of rate of deformation tensor which is famous for its simplicity. In many cases there isn’t an exact solution for non-Newtonian fluids, so they should be solved with a suitable numerical method.

In recent years many non-Newtonian cases have been solved with the LB method. G. Drazer et al. modeled a power-law for a non-Newtonian fluid using this method [1]. They tested the accuracy of the method for a shear thinning and shear thickening power-law fluid in parallel and reentrant geometries. They also found an excellent agreement with the solution obtained by the finite element method.

J. Boyd et al. used LB method to simulate a powerlaw fluid [2]. A second-order accurate LBM for sheardependent non-Newtonian flow was proposed. This method avoids time consuming derivations of the velocity data to calculate the shear flow. S. P. Sullkivan et al. simulated the flow of non-Newtonian fluids through complex random porous media in both two and three dimensions. The power-law model was used for shear-thinning fluids in their study [3].

J. Psihogios et al. applied LB method to simulate the flow of non-Newtonian shear-thinning fluids in threedimensional digitally reconstructed porous domains. They showed that the LB method combined with digital reconstruction is a powerful tool for the study of nonNewtonian flow in porous media [4].

M. Yoshino et al. simulated power-law and Carreau fluid flows in a reentrant corner geometry and in a three-dimensional porous structure [5]. Mahmud Ashrafizaadeh et al. applied K-L and Casson and Carreau-Yasuda models to the blood flow in a channel and compared the models to each other. They showed the capability of the LB method for such complex fluid flows [6].

The nature of the artery walls isn’t rigid. Therefore, assuming geometry with fixed walls isn’t exact. Using moving boundaries instead of fixed walls solves this problem. P. H. Kao et al. studied the application of various curved boundary treatment schemes to the simulation of fixed and moving curved boundary problems [7]. Y. Sui et al. studied a hybrid immersed-boundary method with multi-block Lattice Boltzmann in order to simulate the moving boundaries interacting with incompressible viscous fluid [8]. The results demonstrate the flexibility of their model and show it as an alternative approach for the simulation of fluid-structure interactions. J. Wu et al. present a new version of the immersed boundary-lattice Boltzmann method (IB-LBM) for simulating incompressible viscous flows around moving objects [9]. It is believed that their method has a potential to effectively simulate incompressible viscous flows around moving objects.

The oscillating pressure in an artery causes expansion in it. With some information about the flow field, fatty deposits in the arteries can be reduced. First, we focus on a micro-channel with local expansion and a fixed wall which is a simple type of blood flow in the artery. Three geometries are chosen for this simulation. Two symmetric squared cavities, two trapezoid cavities, and two staggered cavities are studied. We used a Lattice Boltzmann FORTRAN code for this work. The method is validated by a Poiseuille flow in a channel. Finally the accuracy of LB method is shown. For future works we can extend our code for oscillating flow with moving boundaries in a tube to exactly simulate the blood flow in an artery.

2. Background

The most widely used general viscous constitutive relation is the power-law model. In this model, the local viscosity is a function of a second invariant of the rate of deformation tensor [10]:

where n is the power-law exponent and m is the proportional constant. When n is one, µ shows the viscosity of the Newtonian fluid. As a result, m corresponds to the viscosity of the Newtonian fluid. Note that n must be greater than zero. The exponent n can be interpreted as below:

n < 1     Shear-thinning (pseudo-plastic) fluid

n = 1     Newtonian fluid

n < 1     Shear-thickening (dilatant) fluid

Pseudoplastic, or shear-thinning fluids have a lower apparent viscosity at higher shear rates, and are usually solutions of large, polymeric molecules in a solvent with smaller molecules. Dilatant, or shear-thickening fluids increase in apparent viscosity at higher shear rates.

The rate of the deformation tensor, D, is the symmetric part of the velocity gradient tensor, L.

The velocity gradient tensor in the Cartesian coordinates (x, y and z) is:

As a result, the rate of the deformation tensor can be calculated from:

According to the two dimensional geometry, the II2D is obtained by the following equations:

Consequently, the local viscosity can be calculated:

It is obvious that the local viscosity is a function of velocity derivatives.

3. Numerical Scheme

The Lattice Boltzmann method is a suitable way to solve the complex fluid flows numerically, such as the NonNewtonian fluid flow. Here we introduce the method and apply it to our case.

We apply the D2Q9 model, as it is two dimensional and contains 9 directions. Particle positions are confined to the nodes of the lattice. Variations in momenta that could have been due to a continuum of velocity directions and magnitudes and varying particle mass are reduced (in the simple 2-D model we focus on here) to 8 directions, 3 magnitudes, and a single particle mass. Figure 1 shows the Cartesian lattice and the velocities ea where a = 0, 1,···,8 is the direction index and e0 = 0 denotes particles at rest. This LB method classification scheme was proposed by Qian et al. [11] and is used widely. Because particle mass is uniform (1 mass unit or mu in the simplest approach), these microscopic velocities and momenta are always effectively equivalent. The lattice unit (lu) is the fundamental measure of length in the LBM models and time steps (ts) are the lattice time unit.

The velocity for the particles 1, 2, 3 and 4 is 1 lu/ts, and for the particles 5, 6, 7 and 8 is lu/ts.

We introduce fi distribution function in the ith direction. Accordingly, the macroscopic density of the fluid can be defined as the sum of the distribution function in 9 directions:

Figure 1. D2Q9 model.

Similarly, the macroscopic velocities can be derived from the sum of directional distribution functions of microscopic momentum per density:

where ex(i) and ey(i)are the directions of the ith particle in x and y directions, respectively.

The evolution of a lattice Boltzmann proceeds in two steps that take place during each time step. The first step is a streaming step in which the particles move to new sites according to their previous positions and velocities. Next, the particles collide according to collision rules. The BGK (Bhatnagar-Gross-Krook) approximation is used in the simplest LBM [8,10]. Streaming and collision are combined in the equation below:

Collision of the fluid particles is considered as a relaxation towards a local equilibrium and the D2Q9 equilibrium distribution function f eq is defined as:

where the wi:

The relaxation time, τ, is obtained by the cinematic viscosity, ν, according to the equation below:

In the single relaxation time model, it is usual to set τ = 1, but in non-Newtonian cases, relaxation times for each point change by velocity derivatives. As a result, the relaxation times depend on the coordinates, therefore, local relaxation time is obtained.

4. Boundary Conditionss

According to Sukup and Succi [12,13], we have 3 general kinds of boundary conditions:

• Periodic boundary

• Bounceback boundary

• Von-Newman boundary Here, bounceback and Von-Newman boundaries have been used.

Bounceback boundaries are particularly simple and have played a major role in making LB method popular among modelers. These modelers have shown interest in simulating fluids in domains characterized by complex geometries such as those found in porous media. The uniqueness of it is that one simply needs to designate a particular node as a solid obstacle and no special programming treatment is required. In our case, Bounceback boundaries are applied to the walls. In the walls, as described in the figure 2, the distribution function can be obtained in all directions except the ones that have the directions towards the domain.

According to figure 2, we can calculate the f0, f1, f2, f3, f5 and f6, but the f4, f7 and f8 remain unknown. The bounceback rule says that f4, f7 and f8 can be obtained by considering these three equations:

In other words, after the streaming steps, the bounceback boundary condition sets these unknown distribution functions equal to their opposite directions. Conclusively,

Figure 2. Bounceback boundary condition.

it sets f4 equal to f2, f7 equal to f5 and f8 equal to f5.

Von Neumann boundary conditions constrain the flux at inlet and outlet boundaries. We set a constant x-velocity at the inlet and a fully-developed x-velocity at the outlet. The procedure, in which we can obtain the macroscopic values, is explained as below.

In the inlet boundary nodes, since the velocity is known, the following equations can be written:

In the above equations, three distribution functions and one density are unknown: f1, f5, f8 and ρ respectively. In oder to find these four magnitudes, another equation is needed. As proposed by Zou and He [14], the fourth equation can be written by assuming that the non-equilibirium distribution function of the bounceback condition is used in the direction normal to it.

In conclusion, all four unknown magnitudes are obtained.

For the outlet boundary the procedure is like the inlet one. The only difference between these two is that the velocity magnitudes are calculated in a way in which the velocity gradient in the x-direction and the velovity magnitude in the y-direction are both set to zero. The velocity gradient can be calculated with the first order or the second order approximation.

Figure 3 shows all the boundary conditions of the domain.

Figure 3. Boundary conditions.

5. Results and Discussion

5.1. Validation

In order to validate our method, we applied it the Poiseuille flow in the parallel channel geometry. The exact solution of this flow is brought here to compare to the numerical solution.

where G is the pressure gradient along the channel, H is the channel width, and y is the vertical coordinate. Also, we can rewrite the above equation as a function of average velocity:

Note that is the average velocity in the channel. Figures 4-6 show the velocity profile obtained from the exact solution.

To solve the parallel channel numerically, we set the width of the channel to 80 lu and its length to 400 lu. Average velocity is taken to 0.03 m/s and the Newtonian viscosity is set to 0.002. The Reynolds number is defined as:

Here, the Reynolds number is equal to 1. The boundary condition is constant velocity at the inlet and fully-

Figure 4. X-velocity profile for n = 0.5 (Pseudoplastic fluid).

Figure 5. X-velocity profile for n = 1 (Newtonian fluid).

Figure 6. Velocity profiles for independency investigation.

developed velocity at the outlet. We studied three cases to verify our code:

• n = 0.5

• n = 1

• n = 2 The results have a good agreement with the exact solution of the problem which shows the accuracy of the method.

Solid lines stand for the analytical solution and dash lines are the numerical results using the Lattice Boltzmann method.

5.2. Grid Independency

In order to prove that the results are independent of the grid size, the method is applied to the geometry of the problem with two different grids.

In the first case, the channel inlet width is divided into 50 lattice units and the channel length is divided into 700 lattice units. Then the results are compared to the geometry with twice the width and length (100 lattice unit width and 1400 lattice unit length). The chosen fluid is pseudoplastic with n = 0.7.

The velocity profiles are presented in figure 7 for both grids. The results match each other completely; therefore the solution is independent of the grid. The channel width is non-dimensionlized in order to compare the results better.

5.3. Convergence and Error Discussion

To check the convergence of the code, the following condition is set:

where, X and Y are respectively the number of nodes in x-direction and y direction. N is the whole number of the nodes in the domain (N = X × Y). V is the velocity and i is the iteration number.

The convergence condition is modified in such a way that the error magnitude reaches specified limits. Here we set the convergence error to 10−10 (i.e. −8 in logarithmic scale).The error is plotted in figures 8-9 for the Newtonian and non-Newtonian cases, respectively. As seen in Newtonian case solution converges smoothly while in non-Newtonian flow some oscillations occur.

Figure 7. the error reduction for the Newtonian case.

Figure 8. the error reduction for the non-Newtonian case (n = 0.7).

Figure 9. Non-dimensional pressure drop for the inline case.

5.4. Discussions

The results are discussed in four sections: pressure drop, velocity profile, shear stress and generated vortices.

5.4.1. Pressure Drop Discussion

In order to compare pressure drop along the channel, relative pressure drop is computed. The relative pressure drop is defined as below:

where P and Pin are the average pressure (normal to the channel axis) for each point and the inlet pressure, respectively.

The relative pressure drop is plotted against the channel length in figures 9-11, for inline geometry, staggered geometry and trapezoid geometry, respectively. The relative pressure drop is plotted for different shearthinning fluids with the n equal to 0.5 to 1.

According to this figures, the differences between Newtonian and shear-thinning fluids are distinct. As the fluid becomes more viscous (shear-thinning), the pressure drop decreases rapidly. In the cavity area, the pressure increases because of velocity decrease.

At the inlet of the channel, there is a boundary layer and flow isn’t developed; therefore pressure drop isn’t linear. Pressure gradient for the staggered cavities is lower than the inline ones. Hence, in inline cavities there is a sudden expansion and contraction in the area of the

Figure 10. Non-dimensional pressure drop for the staggered case.

Figure 11. Non-dimensional pressure drop for the trapezoid case.

channel, while in staggered cavities, changing occurs gradually.

5.4.2. Velocity Profile Discussion

The dimensionless velocity is defined as below:

where U is the uniform inlet velocity and V is the flow velocity.

The velocity profile is sketched along the mid-line of the cavities. These lines are shown in figure 12.

In figures 13-16, the dimensionless velocity profile for inline geometry, the staggered geometry (upper cavity) and (lower cavity), and trapezoid cavity are presented respectively. Figures 17-19 show the streamlines and velocity contours in the domain.

According to the velocity profiles, it can be said that because of the vortices, velocity in the cavities has a negative magnitude. As seen the magnitude of negative

Figure 12. Velocity profile sections.

Figure 13. Non-dimensional velocity for the shear-thinning and inline case.

Figure 14. Non-dimensional velocity for the shear-thinning and staggered case―first cavity (upper cavity).

Figure 15. Non-dimensional velocity for the shear-thinning and staggered case―second cavity (lower cavity).

velocity is small compare to the main flow velocity. Therefore the flow will be trapped in the cavities. But for the case n = 0.5, this negative magnitude isn’t seen. Therefore as n decreases and flow become more viscous, the vortices become weaker. This will be discussed in the vortices discussion section.

The maximum velocity occurs in the centerline for all geometries, but with a different magnitude. The staggered geometry has the highest peak value and the trapezoid one has the least value. This is mainly because of the change happening in the area of the channel and the vortices created inside the cavities, their sizes and shapes. The velocity magnitude in the cavities for in-

Figure 16. Non-dimensional velocity for the shear-thinning and inclinedcavity.

clined geometry is higher than the other two cases.

5.4.3. Shear Stress Discussion

The shear stress contour is computed from the power law equation. As mentioned before, shear stress is obtained from the power law formulation:

where, II2D is the second invariant of the deformation rate tensor, u is x-velocity, and v is y-velocity of the flow.The shear stress contours are presented in figures 22-24, respectively for inline, staggered, and inclined geometries. Considering these figures, shear stress is higher for the lower n numbers. The shear stress contours become wavier in the staggered geometry. In the inclined cases, the shear stress inside the cavity is less than the other two cases.

The shear stress profile is depicted along the lines shown in figure 11. These profiles are demonstrated in figures 20 and 21. In all cases, the magnitude of shear stress increases as we move from center line toward the cavities. Maximum shear stress occurs in the entrance of cavities, then it decreases until it reaches near the wall and increases again at the wall of the cavity. But in the inclined cavity and for n less than 0.7, increasing doesn’t happen in the cavity wall.

Considering the comparison of different shear-thinning fluids, it is seen that shear stress has a higher magnitude for a lower power of n.

The peak of the shear stress moves toward centerline of the cannel as n increases.

Figure 17. Stream lines for the inline cases (from left to right: n = 0.5, n = 0.7 and n = 1).

Figure 18. Stream lines for the staggered cases (from left to right: n = 0.5, n = 0.7 and n = 1).

Figure 19. Stream lines for the inclined cases (from left to right: n = 0.5, n = 0.7 and n = 1).

5.4.4. Generated Vortices Discussion

The streamlines are sketched in figures 17-19 for inline, staggered, and inclined geometries, respectively. The size of the vortex created by Newtonian fluid (n = 1) is greater than the shear-thinning fluids. As n decreases, the size of the vortex becomes smaller, so the center of the vortex moves left and to the upper corner of the cavity. So, reverse flow is lower for the shear-thinning fluid inside the cavity. We can conclude that in the non-Newtonian fluid, flow enters the cavity area more than the Newtonian fluid does.

In the staggered geometry, the condition is the same. The difference between the two is that the flow gets somewhat wavier in the centerline of the channel.

Figure 20. Shear stress profile for the shear-thinning and inline case.

Figure 21. Shear stress profile for the shear-thinning and inclined case.

Figure 22. Shear stress contours for the inline cases. n = 0.5, n = 0.7, and n = 1 (left to right).

Figure 23. Shear stress contours for the staggered cases. n = 0.5, n = 0.7, and n = 1 (left to right).

Vortices cause the flow to become trapped inside the cavities and hence some fatty deposit will be formed in blood flow.

In the trapezoid case, the fluid can easily flow inside the cavity and thus a vortex isn’t created. Figure 21 demonstrates how the vortex is removed. In the Newtonian fluid case, the vortex becomes smaller and moves toward the wall and as a result, the flow enters the cavity area.

Figure 24. Shear stress contours for the inclined cases. n = 0.5, n = 0.7, and n = 1 (left to right).

For n = 0.7, this vortex is small but isn’t yet completely disappeared. In the case n = 0.5, has the cavity has no vortex and fluid flows smoothly in the cavity area.

Therefore, by making the flow more pseudo-plastic, vortices will be removed.

6. Conclusions

The non-Newtonian pseudo-plastic (shear-thinning) fluid flow is investigated numerically. The D2Q9 Lattice Boltzmann method is used to simulate the blood flow in an artery with local expansion. The two-dimensional micro channel is used as a geometry of the problem. Three cases are investigated: inline cavities, staggered cavities, (antisymmetric) and inline trapezoid cavities.

The method is validated by simulation of flow between parallel plates (Poiseuille flow). Two cases are confirmed: Newtonian and shear-thinning. The validation results are completely similar to the analytical results. Therefore, the method is verified.

The results are discussed in four sections: pressure drop, velocity profiles, shear stress, and generated vortices. Pressure drop is compared between inline and staggered geometries. The inline cases have higher pressure drop than the staggered ones. The reason for this difference is that the expansions in the staggered geometry happen gradually and in two steps.

Velocity profiles are presented for inline, staggered, and trapezoid cases. The negative velocity isn’t seen in the shear-thinning fluid with n = 0.6 and or lower. For the trapezoid geometry, the flow becomes more uniform than other geometries.

The shear stress profile is sketched in the centerline of the cavities. It is realized that the maximum shear stress become smaller for the trapezoid case.

The vortices are shown inside the cavities. In the trapezoid geometry they are smaller. For the cases with n = 0.5 and or lower, the vortices are discarded and there were no reverse flow inside the cavities.

It is proved that by making the fluid flow more pseudo-plastic, vortices are removed. Therefore, the problem of fatty deposit inside the artery will be solved. To have success in real cases, local injection of pseudoplastic blood is proposed.

Finally, it is understood that the Lattice Boltzmann method is a suitable method for simulation of the nonNewtonian pseudoplastic fluid in a microscale channel.


  1. S. Gabbanelli, G. Drazer and J. Koplik, “Lattice Boltzmann Method for Non-Newtonian (Power-Law) Fluids,” Physical Review E, Vol. 72, 2005, p. 046312. doi:10.1103/PhysRevE.72.046312
  2. J. Boyd, J. Buick and S. Green, “A Second-Order Accurate Lattice Boltzmann Non-Newtonian Flow Model,” Journal of Physics A: Mathematical and General, Vol. 39, 2006, pp. 14241-14247. doi:10.1088/0305-4470/39/46/001
  3. S. P. Sullivan, L. F. Gladden and M. L. Johns, “Simulation of Power-Law Fluid Flow through Porous Media Using Lattice Boltzmann Techniques,” Journal of NonNewtonian Fluid Mechanics, Vol. 133, No. 2-3, 2006, pp. 91-98. doi:10.1016/j.jnnfm.2005.11.003
  4. J. Psihogios, M. E. Kainourgiakis, A. G. Yiotis, A. Th. Papaioannou and A. K. Stubos, “A Lattice Boltzmann Study of Non-Newtonian Flow in Digitally Reconstructed Porous Domains,” Transport in Porous Media, Vol. 70, No. 2, 2007, pp. 279-292. doi:10.1007/s11242-007-9099-2
  5. M. Yoshino, Y. Hotta, T. Hirozane and M. Endob, “A Numerical Method for Incompressible Non-Newtonian Fluid Flows Based on the Lattice Boltzmann Method,” Journal of Non-Newtonian Fluid Mechanics, Vol. 147, No. 1-2, 2007, pp. 69-78. doi:10.1016/j.jnnfm.2007.07.007
  6. M. Ashrafizaadeh and H. Bakhshaei, “A Comparison of Non-Newtonian Models for Lattice Boltzmann Blood Flow Simulations,” Computers and Mathematics with Applications, Vol. 58, No. 5, 2009, pp. 1045-1054. doi:10.1016/j.camwa.2009.02.021
  7. P. H. Kao and R. J. Yang, “An Investigation into Curved and Moving Boundary Treatments in the Lattice Boltzmann Method,” Journal of Computational Physics, Vol. 227, No. 11, 2008, pp. 5671-5690. doi:10.1016/
  8. Y Sui, Y. T. Chew, P. Roy and H. T. Low, “A Hybrid Immersed-Boundary and Multi-Block Lattice Boltzmann Method for Simulating Fluid and Moving-Boundaries Interactions,” International Journal for Numerical Methods in Fluids, Vol. 53, No. 11, 2007, pp. 1727-1754. doi:10.1002/fld.1381
  9. J.Wu, C. Shu and Y. H. Zhang, “Simulation on Incompressible Viscous Flows around Moving Objects by a Variant of Immersed Boundary-Lattice Boltzmann Method,” International Journal for Numerical Methods in Fluids, Vol. 62, 2009, pp. 327-354.
  10. Ch. W. Macosko, “Rheology Principles, Measurements and Aplications,” 1st Edition, Wiley-VCH, 1994.
  11. Y. H. Qian and S. Chen, “Finite Size Effect in Lattice-BGK Models,” International Journal of Modern Physics C (IJMPC), Vol. 8, No. 4, 1977, pp. 763-771. doi:10.1142/S0129183197000655
  12. M. C. Sukop and D. T. Thorne, “Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers,” Springer, Heidelberg, 2006.
  13. S. Succi, “The Lattice Boltzmann Equation for Fluid Dynamics and Beyond,” Clarendon Press, Oxford, 2001.
  14. Q. Zou and X. He, “On Pressure and Velocity Flow Boundary Conditions and Bounceback for the Lattice Boltzmann BGK Model,” Physics of Fluids, Vol. 9, 1997, pp. 1591-1598. doi:10.1063/1.869307