## Abstract

We present a lattice Boltzmann approach for the simulation of non-Newtonian fluids. The method is illustrated for the specific case of dilute polymer solutions. With the appropriate local equilibrium distribution, phase-space dynamics on a lattice, driven by a Bhatnagar–Gross–Krook (BGK) relaxation term, leads to a solution of the Fokker–Planck equation governing the probability density of polymer configurations. Results for the bulk rheological characteristics for steady and start-up shear flow are presented, and compare favourably with those obtained using Brownian dynamics simulations. The new method is less expensive than stochastic simulation techniques, particularly in the range of small to moderate Weissenberg numbers (W*i*).

## 1. Introduction

In the last two decades, the lattice Boltzmann method (LBM) has emerged as an important tool to simulate the hydrodynamics of Newtonian fluids [1–4]. In recent years, the LBM has been extended beyond the Newtonian realm in a few important cases [5,6]. However, it may be safely said that, from an application point of view, the success story of Newtonian fluids is yet to be repeated in the case of non-Newtonian fluids. One of the important reasons underlying this lack of success is the inability to incorporate elements of an internal microstructure in the lattice Boltzmann framework. Almost all attempts to incorporate viscoelasticity in the lattice Boltzmann framework have revolved around the addition of a body force to the generalized Boltzmann/BhatnagarGrossKrook (BGK) description [7,8]. However, it was shown by Wagner [9] that this straightforward approach is not expected to give the correct nonlinear rheology (the nonlinearity inherent in the convected derivatives that are an integral part of the constitutive equations for polymer solutions). An alternate, more appealing, approach, based on the lattice Fokker–Planck equation, was initiated by Moroni *et al*. [10]. In this approach, the discrete lattice representation of the Fokker–Planck collision term, expanded in a truncated sequence of Hermite polynomials, allows for the implementation of a LBM. However, an analogous approach for polymer molecules would be limited to weak flows, since a similar expansion of the configuration distribution function would require the inclusion of a prohibitively large number of terms with an increasing departure from equilibrium (strong flows).

The present work is inspired, in part, by the aforementioned approach of a direct discretization of the Fokker–Planck equation. However, we have found a simpler approach of discretizing this equation. The crucial new idea is to recognize that the Fokker–Planck is, in general, an advection–diffusion equation in the relevant configurational space, and the application of the LBM for such equations is well developed (e.g. [11–13]). In what follows, we use these facts to develop a lattice Boltzmann approach to polymer dynamics. In our approach, unlike [10], the discrete lattice representation is employed for the velocity degrees of freedom while the distribution function of interest exists in configuration space. Since the former are almost always close to equilibrium even for the strongest flows typically encountered in applications (the momentum relaxation time being *O*(10^{−6}s) or smaller; see [14]), our approach allows for the simulation of strong departures from configuration equilibrium with a relatively inexpensive discrete representation of the velocity distribution function.

The work is organized as follows. In §2, we review the lattice Boltzmann approach to solving the advection–diffusion equation. In §3, the lattice Boltzmann implementation of the Fokker–Planck equation for a finitely extensible nonlinear elastic (FENE) dumb-bell will be outlined; the dilute solution rheology follows directly from the configuration distribution function for a single dumb-bell. In §4, the material functions obtained, for both steady and start-up simple shear flow, via the present lattice Boltzmann approach will be compared with the results of Brownian dynamics simulations. Finally, in §5, we summarize the main results of the study, and briefly comment on future extensions of the present work.

## 2. The lattice Boltzmann method for the advection–diffusion equation

The LBM to solve the advection–diffusion equation relies on the fact that breaking a conservation law, via the incorporation of a relaxational term, leads to the appearance of diffusion in the hydrodynamic limit. Our presentation of the advection–diffusion LBM broadly follows Karlin *et al.* [11]. We will work with the D3Q27 lattice, obtained as a tensor product of the one-dimensional D1Q3 lattice, for which the discrete velocities and the associated weights are
2.1
and
2.2
where *m* is the mass, *T*_{0} is the temperature and *k*_{B} is the Boltzmann constant; for such a lattice, the choice of the *H*-function [15,16] is given by
2.3
The discrete-velocity local equilibrium is the minimizer of the *H*-function under the constraints of fixed mass and momentum density,
2.4

For the D3Q27 model, using the factorization property of the lattice, the explicit expression for the equilibrium reads [16]
2.5
where we note that the equilibrium distribution is completely parametrized in terms of mass density *ρ* and velocity **u**.

For such a discrete velocity model, we start with the modified discrete Boltzmann–BGK equation for the distribution function *f*_{i} as
2.6
where the density *ρ* is computed from *f* using equation (2.4) and **v** is some imposed velocity (rather than being given by the first moment of *f*). As a result, we only have the mass conservation equation,
2.7
Momentum is no longer a conserved variable, and the momentum density, *ρ***u**, in equation (2.7) may now be obtained using a Chapman–Enskog expansion of the form . It may then be shown that the continuity equation in the long-time long-wavelength limit reduces to the following advection–diffusion equation [11]:
2.8
where the diffusion coefficient is *D*=*τc*^{2}/3. Note that ** v** is the advection velocity in the configuration space of interest (for instance, the velocity of the end-to-end vector for a FENE dumb-bell; see §3), and the underlined term represents the leading-order error in the formulation. Although both the terms on the right-hand-side of equation (2.8) are

*O*(

*τ*), the ratio of the error to the diffusive term is

*O*(

*v*/

*c*)

^{2}, and the simulation would reproduce the advection–diffusion dynamics of interest accurately provided Ma≪1, where Ma=

*v*/

*c*is the Mach number.

Equation (2.6) can be approximated along the characteristic to obtain the discrete space–time evolution equation as
2.9
where the effective relaxation parameter is *β*=*δt*/(2*τ*+*δt*) and
2.10

Thus, once **v** is given, equation (2.9) can be used to solve the advection–diffusion equation
2.11

## 3. Lattice Boltzmann method for polymer dynamics

Herein, we use the procedure outlined in the previous section to obtain the non-equilibrium configuration distribution function, *ρ*(** Q**,

*t*), of a polymer molecule modelled as a FENE dumb-bell [17], an approximation known to remain qualitatively accurate even in complex flows. The non-dimensional Fokker–Planck equation in an imposed homogeneous flow is given by Bird

*et al*. [18], 3.1 where, in accordance with usual convention, we have used

**in place of**

*Q***x**to denote the configuration space variable—in this case, the dumb-bell end-to-end vector. In equation (3.1), is the dimensionless nonlinear restoring force for the FENE connector, given by 3.2 where is the dumb-bell length at full extension. Lengths here have been scaled with the equilibrium extension, (

*kT*/

*H*), and time with

*τ*

_{R}=(

*ζ*/4

*H*), the dumb-bell relaxation time in the Hookean regime;

*ζ*is the bead friction coefficient and

*H*the Hookean spring constant. Finally, , the Weissenberg number, is a measure of the departure from equilibrium, with being a representative magnitude of the velocity gradient tensor (

*κ*^{†}) that characterizes the strength of the imposed flow.

Equation (3.1) is evidently an advection–diffusion equation with the imposed velocity, as defined in §2, being given by
3.3
With this choice of ** v**, and the use of the numerical procedure summarized in equations (2.9) and (2.10), solving the BGK equation on the phase-space lattice leads precisely to the advection–diffusion dynamics in configuration space described in equation (3.1).

The accessible ** Q**-space for the FENE dumb-bell is a sphere of radius

*b*

^{1/2}. The simulation domain used in the LBM implementation is a cube, however, and must extend outside the accessible region in order to include a significant fraction of the accessible

**-space (figure 1). The singular FENE potential is, therefore, regularized outside a sphere (of radius**

*Q**L*

_{c}) that lies just inside the true spherical domain (); a Padé potential replaces the true FENE potential in the region outside this sphere.

## 4. Results and discussion

The deviatoric stress in a polymer solution may be written in the form ** σ**=

*σ*_{S}+

*σ*_{P}, where

*σ*_{S}is the Newtonian solvent stress and

*σ*_{P}is the additional stress owing to the polymer molecules. The expression for the deviatoric polymeric stress contribution is given by (see [19]) 4.1 where

*n*is the number density of polymer molecules, and the angular brackets denote an ensemble average. In the dilute limit, when each polymer molecule (modelled as a FENE dumb-bell) is effectively isolated from the others, the average in equation (4.1) is with respect to

*ρ*(

**,**

*Q**t*), the configurational probability density associated with the single dumb-bell end-to-end vector in the imposed flow, governed by equation (3.1) (see §3).

Herein, we present results for the shear viscosity and the two normal stress coefficients, as a function of W*i*, for steady and start-up simple shear flow. As is well known, these three material functions completely characterize the flow behaviour of a non-Newtonian (simple) fluid in a viscometric flow configuration (of which simple shear is a special case), and are defined by
4.2
and
4.3
The transpose of the velocity gradient tensor is given by
4.4
where is a constant for steady simple shear, and for start-up shear, *H*(*t*) being the Heaviside function. In what follows, the material functions obtained from the solution of equation (3.1), using the LBM approach described earlier, are compared with the results from Brownian dynamics simulations (see [20]).

### (a) Steady simple shear flow

Figures 2 and 3 show the growth of the shear viscosity (*η*) and the first normal stress coefficient (*ψ*_{1}) during steady simple shear; also shown are the corresponding analytical asymptotes in the limit , obtained by Warner [17]. The numerical results are in good agreement with the analytical predictions for small W*i*, and exhibit a characteristic shear thinning behaviour for Weissenberg numbers of order unity and larger; a behaviour that arises from the dumb-bells aligning with the flow direction. The results for both *η* and *ψ*_{1} begin to deviate from the Brownian dynamics simulations for large W*i* possibly because of the breakdown of the Chapman–Enskog expansion. Note that the dumb-bell model under consideration neglects hydrodynamic interactions between the beads, and the second normal stress coefficient (*ψ*_{2}) is, therefore, identically zero at steady state.

### (b) Start-up simple shear flow

Figure 4 shows the growth of the shear viscosity (*η*^{+}) and the first normal stress coefficient () for W*i*=1. The results are again in good agreement with those obtained from Brownian dynamics simulations. Figure 5 shows the variation in the second normal stress coefficient (*ψ*_{2}) as a function of time for W*i*=1. Note that *ψ*_{2} goes through a maximum before decreasing to zero for long times (the latter being consistent with *ψ*_{2}=0 for steady shear). The non-zero *ψ*_{2} for intermediate times is a characteristic feature of the FENE model and is accurately captured by our numerical scheme. For the commonly used pre-averaged versions of the FENE constitutive equation (FENE-p and FENE-CR; see [18,21]), the components of ** Q** along the gradient and vorticity directions obey identical stochastic differential equations, and

*ψ*

_{2}is therefore identically zero for all times.

## 5. Outlook

In this paper, we have presented an LBM approach to polymer dynamics, and validated it for steady and start-up simple shear flow by comparing the results obtained with those known from Brownian dynamics simulations. The agreement is good for low to moderate W*i*, and we have found that, at least for this range of W*i*, the convergence is extremely fast. This behaviour is visible in figure 6, which shows that less than 5 per cent error is observed for 25 grid points up to W*i*=2.0. Furthermore, the simulation seems to be converging very fast for low W*i* (<1). In order to extend this agreement to higher W*i*, we might require a restructuring of the underlying Chapman–Enskog expansion as well as a proper treatment of boundary conditions. For the present case, the latter seems to be the major source of error at high W*i*. Minimizing this error might be particularly relevant for extensional flows in order to be able to capture the coil–stretch transition. The proposed approach already complements existing ones based on stochastic simulation techniques; for instance [5], where it has been shown that the lattice Boltzmann formulation has a better scaling as a function of the number of beads *N* in the limit of than Brownian dynamics. Since W*i* is the measure of the imposed deterministic motion to the random diffusive motion, the size of the averaging ensemble, needed to obtain bulk properties to a given degree of accuracy, diverges as , and stochastic simulation techniques become increasingly inefficient in this limit. One, therefore, expects a critical W*i* below which the present LBM approach would become more efficient. The proposed scheme gives reasonably good statistics in the limit of W*i*≤1.25, but is quite memory extensive when compared with Brownian dynamics. A precise estimate of the range of applicability again requires a more detailed examination. For the problems such as turbulent drag reduction, where the feedback of the polymer stresses on the flow field is crucial, the present solver may be coupled with a hydrodynamic solver to determine macroscopic fields in an iterative manner. The latter solver may again be based on the LBM.

## Acknowledgements

We are grateful to the Department of Science and Technology (DST), India, for providing funding for this work via a Ramanujan Fellowship grant.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society