Radial distribution function of Lennard-Jones fluids in shear flows from intermediate asymptotics

Radial distribution function of Lennard-Jones fluids in shear flows from intermediate asymptotics

Luca Banetta, Alessio Zaccone Statistical Physics Group, Department of Chemical Engineering and Biotechnology, University of Cambridge, Cambridge, CB3 0AS, United Kingdom Department of Physics “A. Pontremoli”, University of Milan, via Celoria 16, 20133 Milan, Italy Cavendish Laboratory, University of Cambridge, JJ Thomson Avenue, CB3 0HE Cambridge, U.K.

Determining the microstructure of colloidal suspensions under shear flows has been a challenge for theoretical and computational methods due to the singularly-perturbed boundary-layer nature of the problem. Previous approaches have been limited to the case of hard-sphere systems and suffer from various limitations in their applicability. We present a new analytical theory based on intermediate asymptotics which solves the Smoluchowski diffusion-convection equation including both intermolecular and hydrodynamic interactions. The method is able to recover previous results for the hard-sphere fluid and, for the first time, to predict the radial distribution function (rdf) of attractive fluids such as the Lennard-Jones (LJ) fluid. In particular, a new depletion effect is predicted in the rdf of the LJ fluid under shear. This method can be used for the theoretical modelling and understanding of real fluids subjected to flow, with applications ranging from chemical systems to colloids, rheology, plasmas, and atmospherical science.


I Introduction

This work presents a new analytic resolution of the Smoluchowski diffusion equation with shear to analyze the microstructure of strongly sheared colloidal suspensions. The Smoluchowski equation provides a means to determine the pair distribution function, or radial distribution function (rdf) , which gives the probability to find a particle at a certain distance r with respect to a reference (tagged) particle Hansen_book (); Larsen ().

The rdf is usually influenced by contributions which can be divided into Brownian-induced and shear-induced effects. The first class includes diffusion (Brownian motion) and inter-particle interactions, while the second class includes the various effects due the flow field.

It is also important to consider the effect of hydrodynamic interactions within the liquid medium: the presence of a second particle or more particles will contribute both a shear-induced hydrodynamic interaction BatchelorGreen () as well as a lubrication contribution Honig (), if we are under Stokes regime.

The relative importance of shear-induced to Brownian effects is compactly described through the Péclet number which, for particles of equal radius , takes on the following form:


where is the viscosity of the medium, the Boltzmann constant and the absolute temperature: if then the flow field contribution is going to be dominant; as a consequence the Brownian motion is going to be more important if .

The analysis of this problem started with the pioneering work of Smoluchowski Smoluchowski () who evaluated the two asymptotic limits of (i) purely Brownian and (ii) purely convective dynamics, that were used to determine the coagulation rate for the two limits, respectively.

Subsequently, there have been several studies addressing the complex interplay of Brownian-induced and shear-induced contributions to the .

Earlier studies adopted various approximation schemes and focused on hard-sphere suspensions under weak shear flows Ronis (); Hess (); Dhont (); Blaw (). Motivation for these earlier studies came from pioneering experiments on hard-sphere colloids where the structure factor distorted by the shear flow was measured with optical techniques Ackerson (); Clark (). A further motivation comes from rheology: well documented behaviors such as shear thinning of suspensions have their origin in the distortion of the rdf Hess (); Blaw (). The rheology of glassy systems under shear can also be described using the Smoluchowski equation with shear as input for the dynamics Brader (); Cates (). Finally, the microstructure of complex fluids under shear flow is also an essential input for studies which address the dynamics of phase transitions in those systems Onuki ().

A parallel line of works focused on solutions to the Smoluchowski diffusion-convection equation as a means to obtain the coagulation rate of colloids, thus focusing on realistic physico-chemical interactions between the particles VandeVen (); Feke ().

Several works have addressed the same problem of finding the rdf for hard-sphere suspensions that are subjected to strong shear flows: Batchelor and Green Batchelor () found the spherically averaged pair distribution function which depends on the radial distance between the particles. It is characterized by an exponential trend which diverges at the contact between the target and the reference particles. Twenty years later, Brady and Morris MorrisBrady () presented an analytical study on the microstructure of strongly sheared suspensions through perturbative methods and they found a peak in the g( next to the surface of a reference particle whose magnitude is directly proportional to the Péclet number. This result has been reconsidered based on Stokesian dynamics simulations Morris () which found a weaker dependence and highlighted the decrease of the peak with the volume fraction occupied by the particles, with the number density.

Concerning the effect of non-trivial (or non hard-sphere) inter-particle potential on the rdf, which is encountered in most applications, Feke and Schowalter Feke () considered the case of strongly sheared particles interacting via the Derjaguin-Landau-Verwey-Overbeek (DLVO) potential. Their scheme however relies on a numerical evaluation of the g(r) whose starting point is an approximate far-field analytic resolution which does not go beyond the zeroth order, and is therefore of very limited applicability. Further limitations of the approach pf Feke () were discussed in Ref. Zinchenko ().

A systematic analytical method to study the influence of shear flow and physico-chemical interactions on the microstructure of sheared suspensions, is still missing. In the following, we develop a new methodology, based on the rigorous application of intermediate asymptotics, to obtain the rdf as a solution to the Smoluchowski equation with shear flow for the case of Lennard-Jones (LJ) particles. The method provides predictions of the microstructure with non-trivial features due to the interplay between the LJ interaction, the flow field and the Brownian motion. This method can be used in the future to open up new opportunities for elucidating key features of complex fluids, from the rheology of suspensions to problems in atmospheric science. The proposed work goes beyond the current paradigms for hard-sphere systems and presents a solution framework for interacting particles on the example of the 12-6 Lennard-Jones potential which can be easily extended to other important interaction potentials (e.g. Debye-Hueckel) in the future.

Ii Model

The method is based on finding an analytical solution for the Smoluchowski equation with shear which describes the system depicted in Fig.1: we suppose to put a reference particle at the centre of the spherical reference frame and a second particle at a certain position .

A simple shear flow is implemented which causes the second particle to have a certain (relative) velocity with respect to the reference one (or equivalently, this is the velocity of the second particle in the rest frame of the reference particle).

Figure 1: Schematic illustration of a pair of LJ-interacting particles subjected to a simple shear flow where v = (y, 0,0) in Cartesian coordinates.


We model the hydrodynamic disturbance of the flow field around a particle due to the presence of another particle through the adoption of two functions and derived by Batchelor Batchelor () which influence (see Appendix A for more details). Instead, to describe the effect of lubrication forces we use the widely used parameterized function Honig (); ZacconeNess ():


where represents the surface-to-surface distance between the particles.

Next we pack the two effects of interparticle potential and shear flow into an external force term acting on the target particle:


where is the Stokes friction factor.

Iii Derivation

Using the functions introduced in the previous section, we can now write the two-body Smoluchowski equation with shear in vector notation:


is the diffusion coefficient, is the Boltzmann factor with the Boltzmann constant and the absolute temperature.

Next, we make Eq.(4) dimensionless through:


where is the hard-core particle diameter.

The velocity v(r) can be made dimensionless through:


A more detailed discussion of is presented in Appendix A.

Replacing all the previously introduced terms in the Smoluchowski convection-diffusion equation we obtain


Upon replacing the Péclet number we obtain:


Equation (8) is to be solved perturbatively. A perturbative method approach is based on the introduction of a small perturbation parameter , by definition much smaller than unity, which simplifies the analytical treatment of the partial differential equation (PDE) of interest BenderOrszag ().

Focusing on situations where the effect of shear flow is substantial, we fix:


Applying Eq.(9) to Eq.(8) we obtain:


Starting from (10) we apply the linearity of the divergence operators obtaining:


To be as general as possible, we do not neglect the compressibility of the fluid throughout the following manipulations and we consider the divergence of the velocity field to be not null, as a consequence.

Next, we introduce a useful approximation that was proposed in Ref.ZacconePRE2009 (), in order to make the 3D problem analytically solvable. The approximation consists in the application of an angular average, denoted as , over the solid angle, to Eq.(11), leading to:


The result is the following spherically-averaged solution which depends on the radial coordinate only:


Moreover, we use a decoupling approximation: we suppose that the velocity and the pair correlation function are weakly correlated, so that:


A general flow field can be separated into downstream and upstream regions: in the former regions the particle approach each other (compressing sectors of solid angle), so the relative velocity between the particles is negative; instead, in the upstream regions (extensional sectors), the particles move away from each other, so the radial component of the velocity is positive. Within this methodology, one replaces of the actual relative velocity and the flow field divergence with their angular averaged values within downstream and upstream regions.

Here, we will consider, as a result of the averaging procedure, two coefficients: for the downstream and for the upstream zone, which descend from the angular averaging, and are explicitly introduced and defined in Appendix A. The two coefficients define the influence of the angular coordinates on the radial relative velocity and the divergence of the flow field. Thus for the compressing quadrants we have:


and analogous expressions are found for the extensive quadrants.

Next, we arrive at the following expression:


Then, we isolate the dependence on the radial coordinate only, obtaining


Finally, we put the equation in the following final form which is most convenient for the subsequent perturbative treatment:


Since Eq.(18) is a second order differential equation we need two boundary conditions (BCs). The first one is the usual BC in the far-field


The second BC expresses the fact that the radial flux is null when the two particles are in direct contact:


where is a value of radial distance sufficiently close to the reference particle. In our calculations we will take .

From inspection of Eq.(18), it can immediately be seen that the perturbation parameter is linked with the highest order derivative of the ordinary differential equation (ODE). This means that we are dealing with a singular perturbation problem and it must be solved by the application of the boundary layer theory BenderOrszag (). The approach consists of the evaluation of two different power series related to two different regions of the domain. One region is the outer layer (in this case, farther away from the reference particle), where the solution is slowly changing with , and the inner region (closer to the reference particle), usually called boundary layer, where the solution is steeply and very rapidly changing with the radial coordinate.

iii.1 Outer solution

We write the outer solution as a power series in the small parameter :


We will carry out our derivation up to first order in the perturbative series.

The initial step is the assumption of the solution to be characterized by low values of its derivatives. As a consequence, the derivatives become negligible if they are multiplied by a small parameter such as ; this concept can be implemented by imposing in Eq.(18).

It should be emphasized that for the outer solution we obtain a first order boundary value problem, so we are allowed to impose only one of the two BCs Eq.(19) and Eq.(20). Since we know that the boundary layer is close to we choose Eq.(19), thus ending up with


which leads to the following solution (derived with full details in Appendix B):


Equation (23) is formally identical to the solution proposed by Batchelor and Green BatchelorGreen (), an evidence of the good reliability of the proposed method.
Now, we will evaluate the first order term ; starting from Eq.(21) we need to evaluate the first and second derivative of :


Then we replace Eq.(24) and Eq.(21) in Eq.(18) and we group all the terms linked with the same power of the perturbation parameter. Since the perturbation parameter can never be null, the only way to find the n-th order coefficient of Eq.(21) is to impose the coefficient related to the n-th power of to be zero.

In particular, to find , in Eq.(25) we isolate the coefficient which multiplies the first power of :


The boundary condition in Eq.(25) represents the independence of far-field boundary condition of the Péclet number.

It is possible to solve Eq.(25) analytically following the steps reported in Appendix B leading to:


where is a function introduced in order to simplify the structure of the first order term; its full expression has been explicitly defined in Appendix B. Finally, we can obtain the first order angular-averaged outer solution:


iii.2 Inner Solution

We now focus on the small section of the domain where the solution varies dramatically with respect to variations in . The first step is to apply a change of variable called inner transformation. Since it is known that the region where the rdf varies the most is close to the reference particle, we can define the inner coordinate as:


where is the order of magnitude of the width of the boundary layer.

Before turning to the evaluation of the inner solution we must find the relationship between and the perturbation parameter itself; this procedure is called dominant balancing BenderOrszag ().

At first, we need to apply the inner transformation to Eq.(18), giving


The crucial point of the dominant balancing is to find the asymptotic behavior of all the functions in the above ODE with respect to and . For this reason we need to find a set of “gauge functions” representing the asymptotic trend of every term in Eq.(29) as VanDyke (); all the mathematical steps are reported in Appendix C and it has been found that


where , , and are functions supposed to be bounded as , so they will be considered as . It is important to highlight that the terms related to the inter-particle potential have been evaluated considering a 12-6 Lennard-Jones potential.

Upon replacing all the above expressions in Eq.(29) and rearranging we obtain:


The gauge functions representing the orders of magnitude of the terms in Eq.(31) are listed in Table 1.

(i) (ii) (iii) (iv) (v)
O( O() O( O(1) O()
Table 1: coefficients appearing in Eq.(31)

as they appear in the equation from left to right.

The coefficient (iv) multiplying the velocity has to be because the relative velocity between the particles can never be null. The aim of this procedure is to find the pair of coefficients of Eq.(31) which counts the most as , through a ”trial and error” procedure; immediately we can notice that, since both and are much less then unity, (iii) is going to be for sure less dominant than the other terms, so it will be discarded.

Let us now suppose that the pair of terms that dominate are (ii) and (iv): this choice is unreasonable since we have assumed at the beginning of the derivation that . Next, if we assume that (i) and (ii) are dominant as , this means that . In this case (i), (ii), (iv) and (v) become respectively ), O(), O(1) and O(): it is evident that the second last term, considered to be negligible, is actually the most important one as , so even this hypothesis is unreasonable.

Finally, we suppose that (i) and (iv) are the dominant terms, which leads to . This final assumption is the correct one because (i), (ii),(iv) and (v) become , , and O() respectively; in this case the excluded terms are negligible compared to the other two as ; for this reason we will consider the following expression for the boundary layer:


Equation (32) provides an estimate of the width of the boundary layer which is in agreement with the literature Feke (), but it has been derived here in a rigorous way for the first time.

Based upon the previous results, we can write the differential equation for the evaluation of in its final form as:


Since we want the inner solution to be a perturbation series based on the powers of a small parameter, we can adopt the result previously obtained in Eq.(32) to describe the structure of as :


Now it is possible to evaluate the leading order term of Eq.(34); under the same hypothesis made for the outer solution we impose in Eq.(33) ending up with the following boundary value problem:


The solution of Eq.(35) is straightforward and we find:


where and are integration constants to be evaluated later.

Considering the behavior of , we need to replace the first and the second order derivative of the inner solution, defined as


into Eq.(33). Then, we apply the same procedure that we used for the outer solution: we group together all the terms linked with the same power of and we set them separately equal to zero. In the case of the first order term this means


The above equation can be analytically solved through the mathematical steps reported in Appendix D, ending up with the following expression:


iii.3 Integration constants evaluation

To summarize, we have evaluated two different series and which describe the behavior of the solution in two different adjacent sections of the integration domain. The final step to obtain the analytical solution of Eq.(18) is the evaluation of the integration constants , , and present in the inner solution; the full procedure is summarized in Fig. 2.

Figure 2: Block diagram with the fundamental steps for the evaluation of the integration constants within

Since we have four unknown parameters we need four equations to determine them: the first one will be the condition of zero flux at the reference particle surface Eq.(20), while the other three can be obtained from the patching procedure BenderOrszag ().

The general principle is as follows. We start from two solutions which share a common border. If one of the two solutions is known and the other has constants to be evaluated, it is necessary to apply a continuity of order to evaluate the constants.

This principle is suitable for our case since we know the full behavior of the outer solution and we have three remaining conditions to be fixed in order to find the remaining constants. Hence, since we need to fix a second order continuity between and at their shared border, that is . After having obtained the complete structure of the inner solution, we need to group together all the coefficients related to each integration constant; for clarity we will adopt the following mathematical notation to describe the structure of the inner solution:


So, can be written as


To find the integration constants it is necessary to go back to the variable :


Next, we need to write the first and the second order derivative of :


It is possible to evaluate all the integration constants through the resolution of the following linear system which is dictated by the patching procedure,


together with the application of Eq.(20) recalled briefly here


The result of Eq.(45) plus Eq.(20) is the evaluation of the four integration constants , , and as functions of the Péclet number which will provide the final form of .

In the following section we will present the results related to the spherically-averaged rdf which will be given by the patching of the inner solution inside the boundary layer with the outer solution outside the boundary layer:


Iv Results

Our main interest is the analysis of the rdf inside the compression quadrants, which contains the more interesting physics, where and . On the other hand, the only available data to verify this analytical approach are given by simulations of hard-spheres where the rdf is evaluated for .

For this reason and for validation purposes, we include calculations for related to the extensional quadrants, where , into the calculations and then average over the solid angle to obtain the rdf averaged over all sectors. From this point onward we will refer to the rdf averaged over all as , to be distinguished from the rdf averaged only over the compressional sectors, which will be denoted as .

iv.1 Hard Spheres

This case corresponds to setting in the above derivation.

iv.1.1 Comparison with simulation data for Hard Spheres

In Fig.3 we compare the behavior of the ) with the results obtained through Stokesian dynamics simulations Morris () for strongly sheared suspensions at different volume fractions.

Figure 3: Comparison of evaluated from (47) with simulations (symbols) of Hard-Sphere suspensions at Pe = 1000 from Ref.Morris ().

It is seen that, at sufficiently long distances from the surface of the reference particle, the match between the theory and both sets of simulation data is excellent. It is possible to notice that the peak next to decreases with the volume fraction and our theory predicts a peak that is lower than both of them. This feature makes perfect sense because the theory is supposed to be valid for more dilute systems than those of simulations from Morris (). This result is an important improvement, also with respect to previous theories MorrisBrady () where the prediction of the rdf close to the surface of the reference particle was somewhat overestimated.

iv.1.2 Effect of the Péclet number

We are interested in determining how the rdf evolves as a function of the Péclet number.
In Fig. 10 the behaviors evaluated for the Hard-Sphere model with different values of Pe are plotted. We can immediately notice that the peak of the decreases with the Péclet number: this statement is physically meaningful because, as the compressing effect provided by the advective term is getting weaker, the probability to find the second particle close to the reference one decreases. Moreover, upon decreasing the Péclet number, the region where the balance between the advective and the Brownian contribution is not negligible increases: this is also reflected by the increasing width of the boundary layer which is inversely proportional to the Péclet number.

Figure 10: Effect of the Péclet number on the behavior of for the Hard-Sphere model.

iv.2 Radial Distribution Function of the Lennard-Jones fluid under shear flow

Finally, we consider the influence of a non-trivial interaction potential on the rdf of the sheared fluid. In particular, we will focus on the the 12-6 Lennard-Jones (LJ) interaction potential, which forms a paradigm for many liquids:


where is the minimum of the interaction potential between the particles.

It is expected that, for large values of the Péclet number, the effect of the LJ potential is totally negligible since the shear flow contribution is dominant. This feature is clearly demonstrated in Fig.11 where the for the calculation including the LJ potential is exactly the same as for the hard sphere calculation.

Figure 11: Effect of the interaction potential on at high values of the Péclet number (Pe = 1000).

On the other hand, if the Péclet number is sufficiently small, as the one chosen for Fig 12, we can clearly see the effect of the interaction potential at different values of . To this aim, we can divide the radial domain into three different regions.

In the radial region furthest away from the reference particle, we see a significant undershoot and a minimum in , which becomes deeper upon increasing the attraction energy of the LJ potential. Interestingly, this undershoot is absent in the hard sphere case, and is an original prediction of the theory developed here for the LJ fluid. From a physical point of view, this minimum represents a depletion region which is necessary in order to balance (from a continuity point of view) the strong accumulation of particles at shorter distance caused by the synergy between the attractive force of the LJ and the action of the shear flow, both of which are pushing the particles close to the reference one. Indeed, at closer distance we find a broad and pronounced peak which reflects the increased probability of finding particles in that region, due to the action of both the attractive interaction and the shear flow.

In other words, there is a region corresponding to the peak where the particles are pushed towards the reference one so strongly that there must be a nearby region at further distance where particles are depleted. Interestingly, we also note that with the LJ potential the accumulation peak becomes smooth and non-singular, in contrast with the peak of the hard sphere system, which, instead, features a singularity at the point of maximum of the peak. This feature is to be attributed to the softness of the short-range repulsive part of the interaction of the LJ potential.

It is important to highlight that the depth of the undershoot, the height of the peak and the slope of before and after the maximum, all increase with the attraction parameter : a deeper well of potential corresponds to stronger attractive long-range interactions and also, at the same time, to stronger short-range repulsive contributions.

Figure 12: Effect of the interaction potential on at low values of the Péclet number (Pe = 5).

V Conclusion

In this work we have proposed a method to analytically solve the two body Smoluchowski equation with shear flow in a spherical reference frame by means of an intermediate asymptotics method. The result of this work is the evaluation of the radial distribution function (rdf), which describes the microstructure of colloidal suspensions under a simple shear flow. The reliability of the method across a wide range of conditions of Péclet numbers and in analyzing the influence of an attractive interaction potential on the behaviour of the rdf has been demonstrated. The theory predicts an important feature of the rdf of attractive fluids in compressive quadrants: the presence of a pronounced depletion effect resulting in a minimum or undershoot in the rdf at radial separations right after the accumulation peak. This effect is already visible at modest values of the attraction energy, and becomes more and more important upon increasing the attraction. This new effect may have important consequences on the rheology (e.g. viscosity, non-Newtonian behaviors etc) of colloidal suspensions, which will be explored in future studies.

The method presented here is fairly general and can easily be extended in future work to systems of great interest such as plasmas Larsen (); Hansen (); Rosenbluth (), the rdf of denser and ordered systems Yurchenko (), droplet clustering in atmospheric flows Glienke (), and nucleation and crystallization under shear flows Mura ().

Vi Acknowledgements

Prof. Massimo Morbidelli is gratefully acknowledged for many inspiring discussions and for providing motivation to study this problem. L.B. gratefully acknowledges financial support from Synthomer UK Ltd.

Appendix A Angular averaging

In this section we describe the procedure where we describe the angular averaging procedure with which we evaluate and . First we introduce the spherical components of the dimensionless velocity in a simple shear flow whose equivalent Cartesian coordinates are () shear_flow_spherical_coordinates ():


where and are functions representing the effect of the hydrodynamic disturbance along the radial and angular coordinate, respectively. Their values can be taken from the literature Batchelor () and, in order to use them in the present analytical calculations, they are fitted through the following algebraic expressions Melis ():


Our goal is to evaluate the average radial velocity in the area where the particles are approaching each other, which means the ensemble of angular coordinates .

It is found that the above mentioned condition is satisfied, for , , and . Now we apply the angular average obtaining:


Through this procedure we can obtain


To find the upstream region we need to impose , which is given by , and . Applying the same procedure seen before for we obtain: