Thermodiffusion of rigid particles in pure liquids

Thermodiffusion of rigid particles in pure liquids

Physica A 390 (2011) 1221–1233 Contents lists available at ScienceDirect Physica A journal homepage: www.elsevier.com/locate/physa Thermodiffusion ...

278KB Sizes 2 Downloads 70 Views

Physica A 390 (2011) 1221–1233

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Thermodiffusion of rigid particles in pure liquids Daniel Lhuillier Institut Jean le Rond d’Alembert, CNRS (UMR 7190) and UPMC Univ Paris 6, F-75005 Paris, France

article

info

Article history: Received 5 July 2010 Received in revised form 1 December 2010 Available online 21 December 2010 Keywords: Thermodiffusion Colloidal suspensions Irreversible thermodynamics Hydrodynamics

abstract Thermodiffusion of particles suspended in a pure liquid is a thorny problem which has not yet received a solution admitted by all the different communities interested in. We approach the subject with macroscopic tools exclusively, hydrodynamics and irreversible thermodynamics. These tools have proved their relevance for molecular mixtures and the Soret effect, and we here extend them to suspensions of particles with supra-molecular size. In particular, we obtain the momentum balance of the particulate phase from which are deduced all the physical phenomena inducing a migration of the particles relative to the carrier fluid. Focussing on thermodiffusion, we show that the osmotic pressure is irrelevant and that thermodiffusion cannot have but two distinct origins : the temperature dependence of the stress associated with the distorted particle microstructure and a fluid–particle interaction force involving the temperature gradient. For deformable particles, it is well known that the origin of the fluid–particle temperature gradient force is the temperature dependence of the surface tension. For rigid particles, we suggest it stems from the temperature dependence of the small density jump, the carrier liquid displays close to the particle’s surface. © 2010 Elsevier B.V. All rights reserved.

1. Suspensions of particles vs. molecular mixtures The main differences between suspension of particles in a carrier fluid and mixture of solute molecules in a solvent stem from the supra-molecular size of the solute. If the size of the solute exceeds about 0.1 µm, it is called a particle and the particles can be considered as pieces of a first continuous medium. If the carrier fluid is not a rarefied gas, it can be considered as a second continuous medium embedding the particles. Hence, a suspension is made of two juxtaposed continuous media. However, it is clear that the transformation of wet air (a molecular mixture) into fog (a suspension) does not occur abruptly for a particular drop size and that the above-mentioned size is an order of magnitude only, suggesting an approximate condition for considering a new type of mixture. That mixture of two juxtaposed continuous media is then transformed into a single continuum of two interpenetrating media by means of averaging. The averaging process has far-reaching consequences. While the concept of sharing the total mass (associated with mass fractions) is clear for both molecular mixtures and suspensions, the concept of sharing the total volume (witnessed by volume fractions) is justified for suspensions but not for molecular mixtures because the true volume of a molecule is not defined unambiguously contrary to the true volume of a particle. As an example, the total mass ρ per unit volume can be written in all cases as

ρ = np mp + nf mf

(1)

with nk , the molecular number density and mk , the molecular mass. But it is only when dealing with a suspension that one can write

ρ = φρp + (1 − φ)ρf ,

E-mail address: [email protected]. 0378-4371/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2010.12.010

(2)

1222

D. Lhuillier / Physica A 390 (2011) 1221–1233

where φ = np θp is the particle volume fraction, θp the mean particle volume, 1 −φ the fluid volume fraction and ρk the mass per unit volume of pure phase k. Another remarkable difference between suspensions and molecular mixtures concerns the averaged velocity. Common to both suspensions and molecular mixtures is the mass-averaged velocity v defined as v = cvp + (1 − c )vf ,

(3)

where c = np mp /ρ is the mass fraction of the particles (or solute) while vk is the velocity of phase k. There is a second mean velocity, which is volume-weighted and unambiguously defined for suspensions only u = φ vp + (1 − φ)vf .

(4)

The role played by u in suspensions is best seen when considering mass conservation which is expressed for both molecular mixtures and suspensions as

∂ nk m k + ∇ · (nk mk vk ) = 0. ∂t

(5)

For a suspension of rigid particles suspended in a incompressible fluid mass conservation transforms into volume conservation expressed as

∂φk + ∇ · (φk vk ) = 0. (6) ∂t Because φp + φf = 1 an immediate consequence is ∇ · u = 0 while in general ∇ · v ̸= 0. Despite φ is relevant to suspensions only, one must acknowledge the existence, for molecular mixtures, of quantities looking like volume fractions. De Groot and Mazur [1] introduced θ¯p and θ¯f which have the dimension of a volume and are related to the pressure-derivative of the chemical potentials. Because of the Gibbs–Duhem relation these two quantities are such that np θ¯p + nf θ¯f = 1 which is quite similar to the relation obeyed by the volume fractions. But as was suggested by Felderhof [2] it is only when the size of a solute molecule far exceeds that of a solvent molecule (in other words for suspensions) that θ¯p ≈ θp and consequently np θ¯p ≈ φ . Besides the relevance of quantities like φ and u, suspensions of particles differ from molecular mixtures by the relative ease with which (a) the microstructure of the particles is distorted by non-uniform flows, (b) the relative motion between fluid and particles is created. For a suspension of deformable particles what is meant by microstructure is mainly the particle shape while in case of rigid particles the microstructure is represented by the relative position of the particles characterized by the probability g (R) of observing a particle at position R when another particle is known to be located at the origin. For the solute molecules in a molecular mixture, that probability is so to say frozen at its equilibrium value while in suspensions that probability can be modified by any non-uniform flow, even a weak one. In case there exists direct inter-particle forces (Van der Waals or DLVO forces for example) the flow-induced deformation of the microstructure will have far-reaching consequences on the modeling of the stress associated with these colloidal forces and, as a consequence, on the relative motion induced by this inter-particle stress. More generally, there will be many reasons for observing a relative motion in a suspension. In molecular mixtures the relative velocity occurs in the diffusion flux J = np mp (vp − v) only and its role is neglected everywhere else. In suspensions the relative velocity is upgraded to the status of an internal variable (a variable which vanishes at equilibrium but plays a role out of equilibrium because of its long equilibration time) and one introduces a momentum balance for the particulate phase besides the momentum conservation for the whole suspension. The momentum exchange between the particles and the carrier fluid is linked to the relative velocity. Similarly one can introduce heat exchanges linked to a temperature difference between the two phases and a volume exchange linked to a pressure difference. However to simplify the issue we will assume in what follows that the suspension is depicted by a single temperature T and a single pressure p. That simplified description is justified for rather small and almost incompressible particles (say with a size less than 10 µm) because they have a very small equilibration time for pressure and temperature. Within this simplified model we will emphasize on the various forces exerted on the particles, and in particular the forces associated with a temperature gradient. The paper is organized as follows: Section 2 details how the new status of the relative velocity and the presence of a nonequilibrium microstructure have an influence on the thermodynamics and the hydrodynamics of a suspension. Section 3 gives a general presentation of the two main origins of thermodiffusion in suspensions with due account for previous works while Section 4 presents a new approach of thermodiffusion in the special case of rigid particles suspended in a pure liquid. 2. Thermo-hydrodynamics of a suspension 2.1. Thermo-statics of suspensions The materials from which the particles and the carrier fluid are made of have a free-enthalpy (or Gibbs potential or chemical potential) which per unit mass is given by µ0p (p, T ) and µ0f (p, T ) respectively. Upon mixing together the particles

D. Lhuillier / Physica A 390 (2011) 1221–1233

1223

and the fluid, the free-enthalpy g per unit mass of the suspension becomes g = c µ0p (p, T ) + (1 − c )µ0f (p, T ) + 1g

(7)

where c is the particle mass fraction and 1g is the free-enthalpy of mixing which witnesses to the modifications of entropy and energy resulting from the mixing process. Among the phenomena included in 1g is the entropy of mixing and the potential energy of the inter-particle forces. In a molecular mixture the free-enthalpy of mixing is a function of pressure, temperature and concentration of the solute so that 1g = 1g (p, T , c ). In a suspension of deformable particles an extra term is added to represent the extra energy associated with the particle deformation relative to its equilibrium shape. The deformation is generally represented by a symmetric second-order tensor Xij which vanishes at equilibrium and one writes [3–5] 1

1g = 1g (p, T , c ) +

2

α(p, T , c )X : X ,

(8)

where α is a positive scalar related to the drop-solvent surface tension in case of emulsions and to the entropic elasticity of the polymer coil in case of polymer solutions. The important point is that a new variable is needed to depict the deformed particle shape that is to say the deformed microstructure. For a suspension of rigid particles what is meant by microstructure is much less evident. The spatial distribution of the particles is represented (for dilute suspensions at the least) by the two-particle distribution function P (R1 , R2 ) = P (r − R/2, r + R/2) ≈ n2p (r )g (R), and the microstructure is represented by the probability g (R) for finding two particles with centers separated by R. That probability has an equilibrium value g eq (R) and the microstructure deformation of rigid particles is bound to g (R) − g eq (R). In a molecular mixture the distribution of the solute molecules is always quite close to equilibrium and the above difference vanishes. This is no longer the case for suspensions in which non-uniform flows are able to distort the particle microstructure. The physics is then represented by a free-enthalpy of mixing written as [6,7]

ρ 1g = ρ 1g (p, T , c ) +

n2p



2

[Φ (R) + kB T ln g (R)]g (R)dR

(9)

where all phenomena connected to the number density np and the equilibrium microstructure g eq (R) are gathered into 1g (p, T , c ) as in molecular mixtures. The last term on the right-hand side represents the extra energy linked to the interparticle potential Φ (R) and the extra entropy linked to the distribution g (R) with kB the Boltzmann constant. Note that the equilibrium microstructure is such that g eq (R) = exp(−Φ (R)/kB T ) so that the last term exists out of equilibrium only. The above expression being a bit cumbersome we will simplify it into the symbolic form

1g = 1g (p, T , c , X )

(10)

where X is some internal variable supposed to depict the non-equilibrium microstructure g (R) − g (R). In suspensions the relative velocity is also an important variable. If vp and vf stand for the mean velocity of the particles and the fluid respectively, the kinetic energy per unit volume of the suspension appears as eq

c

vp2 2

+ (1 − c )

vf2 2

=

v2 2

+ c (1 − c )

w2 2

,

w = vp − vf ,

(11)

where v is the mass-weighted mean velocity defined in (3) and w is the relative velocity. If one wants to give a description of suspensions that can be easily compared to that of molecular mixtures we must keep v 2 /2 as the definition of the kinetic energy per unit mass of the mixture and consider the kinetic energy linked to the relative motion as a part of the freeenthalpy. As a consequence the free-enthalpy per unit mass of the suspension is finally written as g = c µ0p (p, T ) + (1 − c )µ0f (p, T ) + 1g (p, T , c , X ) + c (1 − c )

w2 2

.

(12)

Important quantities deduced from (12) are the chemical potentials µk of phase k in the mixture

µp = µ0p + 1g + (1 − c ) µf = µ0f + 1g − c

∂ 1g ∂c

(13)

∂ 1g , ∂c

(14)

and the extended chemical potentials which include the kinetic energy in relative motion 1

µ⋆k = µk + (vk − v)2 .

(15)

2

From (12) one also deduces dg = (µ⋆p − µ⋆f )dc − sdT +

1

ρ

dp +

∂ 1g dX + c (1 − c )w · dw, ∂X

1224

D. Lhuillier / Physica A 390 (2011) 1221–1233

together with g = c µ⋆p + (1 − c )µ⋆f . From the last two results follows the Gibbs–Duhem relation 1

cdµp + (1 − c )dµf = −sdT +

dp +

ρ

∂ 1g dX , ∂X

(16)

which implies that the entropy per unit mass s and the specific volume 1/ρ are defined as

 ∂µk s = csp + (1 − c )sf with sk = − ∂ T c ,p,X   1 c 1−c 1 ∂µk = + with = . ρ ρp ρf ρk ∂ p c ,T ,X 

(17)

(18)

The thermodynamic relations of a suspension are conveniently rewritten in terms of the internal energy e = g + Ts − p/ρ . Of particular interest is the Gibbs relation which appears as ⋆



de = (µp − µf )dc + T ds − pd

  1

ρ

+

∂ 1g dX + c (1 − c )w · dw. ∂X

(19)

It is clear that upon neglecting the microstructural variable X and all terms involving the relative velocity squared one recovers the thermo-static relations of a molecular mixture. 2.2. Local equilibrium The assumption of local equilibrium [8] states that when following a material point with the local velocity v(x, t ) its entropy evolves like in equilibrium. Taking the Gibbs relation (19) into account that assumption is conveniently expressed as T

ds

=

dt

de dt

− (µ⋆p − µ⋆f )

dc dt



p dρ ∂ 1g dX dw − − c (1 − c )w · ρ 2 dt ∂ X dt dt

where d/dt = ∂/∂ t + v · ∇ is the material time-derivative. One eliminates some of the time-derivatives appearing on the right-hand side with the help of the balance equations for mass and internal energy [9] dρ

= −ρ∇ · v,

dt de

ρ

dt

ρ

dc dt

= −∇ · J,

(20)

= −Π : (∇ v)s − ∇ · Q,

(21)

where Q and Π are the Galilean-invariant (non-convective) parts of the total energy flux and momentum flux respectively while J = ρ c (1 − c )w is the particle flux in the material frame moving with v and (∇ v)s is the symmetrized velocity gradient tensor. Local equilibrium thus amounts to an entropy transport written as

[ ] ∂ 1g dX dw ⋆ ⋆ ρT = ∇ · [(µp − µf )J − Q] − (Π − pI ) : (∇ v) − ρ −J· + ∇(µp − µf ) . dt ∂ X dt dt ds





s

(22)

After noticing the identity J·

dw dt

[ ≡J·

dp vp dt



df vf dt

+∇

(vf − v)2 2

−∇

(vp − v)2 2

]

− (w · ∇)v

the entropy balance is transformed into

ρT

ds dt

= ∇ · [(µ⋆p − µ⋆f )J − Q] − (Π − pI − w ⊗ J) : (∇ v)s − ρ

[ ] ∂ 1g dX dp vp df vf −J· − + ∇(µp − µf ) . (23) ∂ X dt dt dt

Except for the term involving X , this result is equivalent to the one obtained in ch.3 Section 4 of Ref. [1]. Any further progress requires the elimination of the last time-derivatives appearing on the right-hand side of the above entropy balance. In other words, one must write down the momentum balances of the particulate phase and the carrier fluid together with the transport equation of the microstructure variable.

D. Lhuillier / Physica A 390 (2011) 1221–1233

1225

2.3. Momentum balances and momentum flux In a suspension the relative motion between the particles and the carrier fluid can be driven by many kinds of forces, much more than those acting in a molecular mixture and which result in Fick’s law. For example the gravity forces enter into play (neglected in molecular mixtures) and the forces exerted by the carrier fluid do not reduce to a drag force. In fact one must take all these different forces into account with two coupled momentum balances for the particulate phase and the fluid phase. Per unit volume of the suspension these coupled equations are [10]

φρp

dp vp dt

(1 − φ)ρf

= ∇ · τ m − φ(∇ p − ∇ · τ ) + F + φρp g, df vf dt

= −(1 − φ)(∇ p − ∇ · τ ) − F + (1 − φ)ρf g.

(24) (25)

The two accelerations are expressed in terms of the convected time-derivative dk /dt = ∂/∂ t + vk · ∇ . The role of these two inertial terms is usually negligible for colloidal suspensions [11] but not for suspensions of larger particles. The external force was written with the gravitational acceleration g because it is by far the most usual case. The force F that appears with opposite signs in the two equations is the total force exerted by the fluid on the particles and pI − τ is a stress due to the fluid but which takes hydrodynamic (fluid-mediated) interactions between particles into account so that τ is a dissipative viscous stress. And τ m is a non-dissipative stress associated with the non-equilibrium microstructure, i.e. with the particle shape in case of deformable particles or with the anisotropic spatial configuration in case of rigid particles. Note that the microstructure stress τ m acts on the particles only while pI + τ is shared between the two phases in proportion to their volume fraction. When summing the two above momentum balances, one obtains the conservation of momentum for the whole suspension in the usual form

ρ

dv dt

+ ∇ · Π = ρ g,

(26)

with the suspension stress, the stress already appearing in the internal energy balance (21), defined as

Π = pI + w ⊗ J − τ − τ m .

(27)

With the above momentum balances the transport of entropy is now transformed into

ρT

ds dt

= ∇ · [(v − u) · τ + (v − vp ) · τ m + (µ⋆p − µ⋆f )J − Q] + τ : (∇ u)s + τ m : (∇ vp )s [   ] ∂ 1g dX 1 1 F −ρ − J · ∇(µp − µf ) − − ∇p + . ∂ X dt ρp ρf ρ c (1 − c )

(28)

The next step is to devise an equation for the time evolution of the microstructure so as to eliminate the last time-derivative appearing in the right-hand side of the above transport of entropy. 2.4. The microstructure stress The non-equilibrium microstructure is represented by the internal variable X . Since Hand [12] it is accepted that the best variable X to represent a microstructure distortion is a symmetric second-order tensor noted below as Xij and some recently published models [13–15] are based on that assumption. As an example, in the model by Stickell et al. [15] the second-order tensor is related to the shape of the free volume surrounding a particle and where no other particle is likely to be present. The point is that, whatever the interpretation given to Xij , its evolution is bound to the particle mean motion in the general form [4]

∂ Xij ∂ 1g + (vp · ∇)Xij = (∇ vp )ik Xkj + (∇ vp )jk Xki − Mijkl , ∂t ∂ Xkl

(29)

and which means that the tensor representing the non-equilibrium microstructure is convected with vp , affinely deforms with the gradients of vp and relaxes towards its equilibrium value Xij = 0 (the simplest type of relaxation involves a relaxation time θ with Mijkl ∂ 1g /∂ Xkl = Xij /θ as can be deduced from the particular expression (8)). As a consequence of the above evolution equation for the microstructure variable

[ ] ∂ 1g dXij ∂ 1g m ρ = τ : ∇ vp − ρ(vp − v) · ∇ Xij − T ∆m , ∂ Xij dt ∂ Xij

(30)

where T ∆m is the dissipation bound to the relaxation towards equilibrium T ∆m = ρ

∂ 1g ∂ 1g Mijkl ∂ Xij ∂ Xkl

(31)

1226

D. Lhuillier / Physica A 390 (2011) 1221–1233

while τ m is a symmetric non-dissipative stress related to the microstructure deformation and defined as

τijm = 2ρ Xik

∂ 1g . ∂ Xjk

(32)

It is possible to refine the evolution equation (29) with the consequence of obtaining slightly modified expressions for the microstructure stress and the dissipation rate but the main features of (30) are conserved. The statistical description of suspensions ends with a quite similar result where instead of an internal variable X one introduces the vector R separating the centers of two particles and considers its pair-distribution function g (R). As explained in Refs. [6,7] instead of (30) one finds from (9)

ρ

d dt

(1g )p,T ,c = τ m : ∇ vp − ρ(vp − v) · ∇(1g )p,T ,c − T ∆m ,

(33)

with T∆ = m

n2p 2



∂ D ∂ (Φ + kB T ln g ) · · (Φ + kB T ln g )g (R)dR, ∂R kB T ∂ R

(34)

and

τm =

n2p 2



(R ⊗ I + 2C) ·

∂ (Φ + kB T ln g )g (R)dR ∂R

(35)

where D(R) is the diffusion tensor of a pair of particles, I the unit tensor and C(R) a third-rank tensor depicting the hydrodynamic interactions between a pair of particles [16]. It is worthy to note that both ∆m and τ m vanish at equilibrium. Based on a pair-distribution function the above results must be modified for concentrated suspensions. Brady [17] derived general results for concentrated suspensions starting from the N-particle Smoluchowski equation and the microstructure tensor one can deduce from his results is

  ∂ 1 ( Φ + k T ln P ) τ m = np (R ⊗ I + RSU · R− ) · B N FU ∂R

(36)

where the bracket ⟨· · ·⟩ means a statistical average over all configurations and involves the N-particle distribution function PN . Moreover the role of hydrodynamic interactions is now written in terms of the resistance matrices RSU and RFU which depend on N-particle configurations and for the definitions of which the reader is referred to Ref. [17]. With result (30) or (33) the entropy balance is now transformed into

ρT

 [  1 1 = ∇ · [(v − u) · τ + (v − vp ) · τ m + (µ⋆p − µ⋆f )J − Q] − J · ∇(µp − µf ) − − ∇p dt ρp ρf ] 1 F − ∇(1g )p,T ,c + + τ : (∇ u)s + T ∆m . c ρ c (1 − c ) ds

(37)

2.5. Entropy and energy fluxes The next issue is to present the entropy balance in conformity with the second principle of thermodynamics

ρ

ds dt

+ ∇ · H = ∆ ≥ 0,

(38)

where ∆ is the entropy production rate and H is the non-convective (Galilean-invariant) entropy flux. In a suspension the entropy of the particles is transported with vp and that of the carrier fluid with vf . Writing the total entropy flux as the sum ρ sv + H implies that H includes the contribution (sp − sf )J. It is important to distinguish the role of that contribution and to write H = h + (sp − sf )J

(39)

where h is some reduced entropy flux. It is to be emphasized that the reduced entropy flux is itself possibly dependent on the relative velocity and when writing (39) we just wanted to extract a trivial dependence of H on the relative motion. The entropy balance is now transformed into

 T

ρ

ds dt

 +∇ ·H

= ∇ · [T H + (v − u) · τ + (v − vp ) · τ m + (µ⋆p − µ⋆f )J − Q] − h · ∇ T + τ : (∇ u)s + T ∆m [   ] 1 1 F 1 − ∇ p − ∇(1g )p,T ,c + . − J · ∇(µp − µf ) + (sp − sf )∇ T − ρp ρf c ρ c (1 − c )

(40)

D. Lhuillier / Physica A 390 (2011) 1221–1233

1227

For that entropy balance to merge into the desired form (38) the energy flux Q must be related to the other fluxes so as to cancel the first term on the right-hand side of (40) Q = T H + (v − u) · τ + (v − vp ) · τ m + (µ⋆p − µ⋆f )J.

(41)

The right-hand side of (40) has now the expected form of a sum of terms written as the product of a flux by a force and we now focus on the last term. 2.6. The osmotic pressure The interphase force F is now split into a thermodynamic force Fth and other forces represented by f F = (1 − φ)(Fth + f).

(42)

The prefactor 1 − φ is for convenience only but implies that the thermodynamic force is defined as Fth = −c ρf

[   ] 1 1 1 ∇(µp − µf ) + (sp − sf )∇ T − − ∇ p − ∇(1g )p,T ,c . ρp ρf c

(43)

The role of the gradient of chemical potential in the thermodynamic force (hence in diffusion phenomena) is well known in molecular mixtures [1] and for suspensions [16]. Here this role is completed by a temperature gradient force, a pressure gradient force and by a microstructural gradient force. For those interested in thermodiffusion the presence of a temperature gradient force is welcome but it happens that the force (sp − sf )∇ T is cancelled by the force coming from the temperature dependence of the chemical potentials. And that cancelling is also true for the pressure gradient force because according to (17) and (18)

∇(µp − µf ) + (sp − sf )∇ T −



1

ρp



1



∇ p = ∇(µp − µf )p,T .

ρf

(44)

This result is of utmost importance: it means that there is no ∇ p contribution to the thermodynamic force provided the role of pressure is correctly taken into account in the momentum balances, and that there is no ∇ T contribution to the thermodynamic force provided the entropy flux is written as in (39). We conclude that the role of the thermodynamic force is limited to concentration-diffusion and to microstructural-diffusion, and that the origin of thermodiffusion must be found elsewhere. We can give an even more compact form to the thermodynamic force with the Gibbs–Duhem relation (16) which results in c ∇(µp − µf )p,T + ∇(µf )p,T = ∇(1g )p,T

(45)

and consequently

[

Fth = ρf ∇(µf )p,T = −ρf ∇ c 2

∂ ∂c



1g

]

c

p,T

.

(46)

And when ρf is a constant (independent of pressure and temperature) it is possible to write th

F

= −∇(p

osm

)p,T

∂ = −∇ φ ∂φ [

2



ρ 1g φ

] (47) p,T

where posm is the so-called osmotic pressure related to the volume fraction dependence of the free-enthalpy of mixing per unit volume ρ 1g. For example, hard spheres with volume V and volume fraction less than about 0.4 are depicted by [18]

ρ 1g =

[ ] 4 − 3φ φ ln φ + φ 2 . V (1 − φ)2

kB T

(48)

We acknowledge some pedagogical value to the concept of osmotic pressure but this concept should be used with care in so far as the pressure or temperature dependence of the osmotic pressure has no influence at all on the thermodynamic force as emphasized in (47). With the above writing of the interphase force F the entropy balance of suspensions is finally expressed in a form compatible with (38)

ρ

ds dt

+∇ ·H=

1 T

[f · (u − vp ) − h · ∇ T + τ : (∇ u)s ] + ∆m .

(49)

The entropy production rate is different from the one that holds for molecular mixtures, with the specific roles played by the hydrodynamic force f and the volume-weighted mean velocity u, together with the dissipation associated with the microstructure relaxation towards equilibrium.

1228

D. Lhuillier / Physica A 390 (2011) 1221–1233

3. Thermodiffusion in suspensions A suspension of particles with a single pressure and a single temperature is described by a set of equations for mass (20)a and (20)b, for internal energy (21) and for momentum (26). These equations have the same form as those describing molecular mixtures [1,9]. The main difference between suspensions and molecular mixtures lies in the (27) and (41) for the overall stress and the energy flux, in the entropy balance (49) and, last but not least, in the momentum balance (24) for the particulate phase. All phenomena related to the relative motion between the particles and the carrier fluid are embodied in the momentum balance for the particles. Combining (24) and (57) with definition (42) for the thermodynamic force, one arrives at a convenient form for that momentum balance

φρp

dp vp dt

− φρf

df vf dt

= ∇ · τ m − (∇ posm )p,T + f + φ(ρp − ρf )g,

(50)

where the left-hand side can be neglected in case of colloidal suspensions. Note that one can combine the gradient of osmotic pressure with the microstructure stress into a single quantity that can be coined the osmotic stress and which writes −posm I + τ m . However, while this combination is permitted by the above special writing of the momentum balance, it no longer exists when the same momentum balance is written like in (24) with (42) and (47) taken into account. Hence, the concept of an osmotic stress is coincidental only while, considered separately, the (equilibrium) osmotic pressure and the (out of equilibrium) microstructure stress are two fundamental quantities related to two distinct and well-defined physical phenomena. Moreover one must keep in mind that the osmotic pressure acts in opposite directions for the fluid phase and the particulate phase while the microstructure stress acts on the particulate phase only. Hence the osmotic pressure must not appear in the overall suspension stress, at variance with the microstructure tensor. Considering the momentum balance (50), one sees that besides the reduced gravity force (responsible for sedimentation) all the phenomena inducing a relative motion must be present either in the microstructure stress τ m or in the hydrodynamic force f, or in the thermodynamic force −∇(posm )p,T . Since the thermodynamic force does not depend on the temperature gradient we conclude that thermodiffusion in suspensions is either due to the temperature dependence of τ m or to a part of f depending on the temperature gradient. These two sources of thermodiffusion will be considered separately. 3.1. Stress-induced thermodiffusion Dhont [19] convincingly suggested that the inter-particle potential Φ is temperature-dependent and he deduced the complete temperature dependence of the osmotic pressure. Unfortunately we have seen above that the temperature dependence of the osmotic pressure does not play any role in thermodiffusion. However, we can take Dhont’s suggestion for granted and apply it to the microstructure stress (35) which becomes temperature-dependent. Hence ∇ ·τ m (T ) will give rise to a temperature gradient force in the particle momentum balance. Since the microstructure stress is an out-of-equilibrium quantity that exists only in presence of non-uniform flows it is possible to propose for it a closure like

τ m = 2ηp (φ, T )(∇ vp )s + · · ·

(51)

where ηp is some effective viscosity of the particulate phase, while the dots stand for possible non-Newtonian contributions [11]. There is no doubt that a part of thermodiffusion comes from the microstructure stress but it is far from easy to calculate its explicit temperature dependence and, being of order n2p or φ 2 at low concentration (see (35)) it does not explain thermodiffusion in dilute suspensions. Moreover because τ m exists only in a flowing suspension it leads to a very special kind of ‘‘convection-induced’’ thermodiffusion which disappears when the flow stops. 3.2. Force-induced thermodiffusion In what follows we focus on the role of the hydrodynamic force f which is the main source of thermodiffusion. An important step is to acknowledge a possible dependence of the reduced entropy flux h on the relative velocity. Such a contribution to the reduced entropy flux cannot take part in the entropy production rate because it is odd in a time reversal. The simplest solution is to suppose that a part of the hydrodynamic force f depends on the temperature gradient and that the two extra terms cancel each other in the entropy production. Hence (49) is compatible with the expressions T h = −λ∇ T + np e⋆ (vp − u), ⋆

f = −np ζ (vp − u) − np e ∇ T /T ,

(52) (53)

where the thermal conductivity λ and the friction coefficient ζ are two positive transport coefficients while there is no constraint on the sign or magnitude of the transport coefficient e⋆ which has the dimension of an energy per particle. It is clear that, upon neglecting the acceleration forces and the microstructure stress in (50) one arrives at np e ⋆ 1 np (vp − u) = − (∇ posm )p,T − ∇T , ζ Tζ

(54)

D. Lhuillier / Physica A 390 (2011) 1221–1233

1229

meaning that the collective diffusion coefficient D, the thermodiffusion coefficient DT and the non-dimensional Soret coefficient TST are given by D=

1



ζ

∂ posm ∂ np



,

TST = T

DT D

=

e⋆

(∂

posm

/∂ np )

.

(55)

Moreover from (39), (41) and (52) can be deduced the expression of the energy flux Q = −λ∇ T + (u − v) · τ + (vp − v) · τ

pp

[ + hp − hf +

np e⋆ c ρf

]

+ (vp − v) /2 − (vf − v) /2 J 2

2

(56)

where hk is the enthalpy per unit mass of phase k. As a consequence hp − hf + np e⋆ /c ρf is nothing but the ‘‘energy of transfer’’ introduced by De Groot and Mazur [1] for molecular mixtures. Or conversely, only the difference between the energy of transfer and hp − hf is involved in thermal diffusion [20]. It is clear from (52) and (53) that the energy e⋆ can be deduced either from a contribution of the temperature gradient to f or a contribution of the relative velocity to T h. The former route was used by most investigators, but the latter one was sometimes preferred [21]. 3.2.1. Deformable particles A few exact results exist concerning fluid particles. Subramanian [22] calculated the force on a very dilute suspension of spherical fluid droplets (with radius R) and found a particular case of (53) f = −np (4π R ) 2

[

1 + 3p/2 ηf 1+p

R

(vp − vf ) +

dγ /dT

(1 + p)(2 + q)

] ∇T

(57)

where γ is the surface tension, p = ηp /ηf is the viscosity ratio while q = λp /λf is the thermal conductivity ratio. Hence, for fluid-like particles thermophoresis is entirely due to the temperature dependence of the surface tension, i.e. a Marangoni effect. Together with the dilute result ∇ posm = kB T ∇ np , the above result for f means e⋆ =

4π R2

∂γ T . (1 + p)(2 + q) ∂ T

(58)

From the momentum balance (50) one can deduce the temperature gradient that compensates for buoyancy effects [23]

∇T =

(1 + p)(2 + q) (ρp − ρf )gR gˆ 3 dγ /dT

(59)

where gˆ = g/g is the unit vector in the direction of gravity. We are not aware of any generalization of (57) to concentrated suspensions. 3.2.2. Rigid particles With a gas as the carrier fluid Despite the interest is on particles surrounded by a liquid, it is worthy to remember that when a spherical particle with radius R is surrounded by a gas it will be acted upon by a force similar to (53) with e⋆ = 4π R2 l0 pH (q, l0 /R)

(60)

where l0 is the mean-free path and p is the gas pressure. When the Knudsen number l0 /R is very large H is a constant independent of the thermal conductivity ratio q [24] while for small Knudsen numbers H ∝ (l0 /R)/(2 + q) [25], and the intermediary regime is now known in detail [26]. A singular case is that of a perfectly heat conducting particle (q → ∞) for which a negative thermophoresis exists for small Knudsen numbers [27,28]. With a liquid as the carrier fluid Referring to (57), the case of solid particles is obtained with p → ∞ and it is clear that if surface tension still plays a role, it involves a different form of the energy e⋆ . Semenov [29] proposed to take into account the London–van der Waals forces between a molecule of the solvent and a molecule of the solid particle material. These liquid–solid forces involve the Hamaker constant A (an energy per couple of particles) and they are effective in a thin boundary layer of thickness l ≪ R around the particles. They provide an excess pressure Π ≈ A/l3 in the boundary layer and when for some reason a gradient of the excess pressure exists along the particle surface it will induce a tangential Stokes flow of the solvent which, when integrated over the particle surface amounts to a slip velocity [30] of order (l2 /ηf )∇ Π between the particle and the bulk of the solvent. Since the friction coefficient is of order ηf R, when the excess pressure depends on temperature one gets e⋆ ∝

R l

[ A

T ∂Π

Π ∂T

]

.

(61)

Semenov neglected the temperature dependence of A but supposed that l behaves like the solvent molecular volume so that (T /Π )∂ Π /∂ T = T αT where αT = −(1/ρf )(∂ρf /∂ T )p is the thermal expansion coefficient of the solvent. Instead of

1230

D. Lhuillier / Physica A 390 (2011) 1221–1233

Semenov’s excess pressure, Parola and Piazza [31] insisted on the anisotropy of the pressure tensor close to the solid particle and the sign of the anisotropy led them to introduce a surface tension γ and a thermophoresis depicted by e⋆ ∝ Rlγ

T ∂(γ l)

[

]

γ l ∂T

.

(62)

Assuming that both γ and l depend on the temperature like the solvent mass density, Parola and Piazza [31] also concluded that the thermal expansion coefficient αT was one of the main factors governing thermophoresis. Note that the Hamaker constant is linked to surface tension as A ∝ γ r02 where r0 is a molecular size so that the two above results are close to each other. But the existence of an anisotropic pressure tensor in the vicinity of the interfaces appears to be a general and useful concept [30] which we want to develop hereafter in a continuum-mechanical framework, offering a different perspective to the approaches used by Semenov and by Parola and Piazza. 4. A Cahn–Hilliard approach to thermophoresis in pure liquids The presence of particle–liquid (i.e. colloid–solvent) molecular forces (with the London–van der Waals forces as an example) concentrated in a thin boundary layer close to the particles can be taken for granted. These forces will induce local modifications of the molecular number density, the magnitude of the modifications depending on the compressibilities of the carrier fluid and the particles. In case of rigid particles the modifications will concern the liquid density only. The range l of the intermolecular forces is usually much larger than the molecular size r0 allowing to represent these forces by a continuum-mechanical quantity, the spatially varying mass density ρ(r , t ) of the carrier fluid in the boundary layer, with value ρ0 at the particle surface and value ρf (the same ρf as in (25) for example) outside the boundary layer. To give some explicit expression for the varying mass density we will use a slightly modified form of the Cahn–Hilliard model [32] supposing that while the bulk has a free-energy per unit volume F (ρf , T ) the boundary layer has a free-energy per unit volume F ⋆ = F (ρ, T ) +

λ(ρ, T ) 2

(∇ρ)2 .

(63)

The boundary layer free-energy not only depends on the mass density but also on the mass density gradients and the magnitude of the positive coefficient λ is related to the thickness of the boundary layer, itself related to the correlation length in the solvent. Since that thickness is usually much smaller than the particle size we simplify the description by considering a one-dimensional configuration in which the interface is located at z = 0 and the bulk of the fluid at z = ∞. The equilibrium density profile ρ(z ) is the solution of the equation [32] F (ρ, T ) − F (ρf , T ) =

λ



2



2 (64)

dz

with the boundary conditions ρ(0) = ρ0 and ρ(∞) = ρf . The sign of the difference ρf − ρ0 depends on the details of the intermolecular forces but the magnitude of the mass density variation is small enough to justify the simple quadratic expression [33] F (ρ, T ) − F (ρf , T ) =

(ρ − ρf )2 2κf ρf2

(65)

where κf = (1/ρf )(∂ρf /∂ p)T is the isothermal compressibility coefficient of the carrier fluid. The equilibrium profile is exponential-like

 z ρ = ρf + (ρ0 − ρf ) exp − ,

l=

l



λκf ρf2 .

(66)

The equilibrium pressure tensor of the Cahn–Hilliard model is anisotropic and gives rise to a surface tension

γ =





λ



0

∂ρ ∂z

2 dz =

l(ρ0 − ρf )2 2κf ρf2

.

(67)

Out of equilibrium the local pressure anisotropy varies along the interface and a gradient of shear stress develops in the boundary layer which is described by [30]

∂ ∂z



∂vx η ∂z



∂ + ∂x



(ρ − ρf )2 κf ρf2

 =0

(68)

D. Lhuillier / Physica A 390 (2011) 1221–1233

1231

where the x-axis is along the interface. The fluid viscosity η depends on the local mass density but the difference ρf − ρ0 is generally too small for this dependence to matter, hence η ≈ ηf in the whole boundary layer. And a stick boundary condition vx (0) = 0 is supposed to hold on the surface of the particle. Two immediate consequences concern the strain rate

∂vx ∂vx 1 d (z ) = (∞) + [γ exp(−2z /l)] ∂z ∂z ηf dx

(69)

and the velocity slip

v = lim



s

z →∞

∂vx vx − z ∂z

 =

1 d(γ l/2)

ηf

dx

.

(70)

When the velocity slip resulting from the temperature dependence of γ l is integrated over the surface of a sphere of radius R one obtains the force acting on a particle and in case of a dilute suspension one gets the force per unit volume of the mixture f = −np 6π ηf R(vp − vf ) − np 2π R(∇ T )0

∂(γ l/2) ∂T

(71)

where (∇ T )0 , the mean temperature gradient in the vicinity of the particle surface, is related to the temperature gradient in the suspension by

(∇ T )0 =

3 2+q

∇T ,

(72)

where q is the thermal conductivity ratio. The total force exerted by the fluid can be written as in (53) with e⋆ =

6π R

∂ 2 + q ∂T T



λ 4

 (ρ0 − ρf )2 .

(73)

This result for rigid particles is to be compared with (58) for deformable ones. The three quantities λ, ρ0 and ρf depend on the temperature and provide three different origins for thermophoresis. The temperature dependence of ρf and the role of the thermal expansion of the liquid has already been stressed [29,31,34]. Less usual is the temperature dependence of ρ0 which witnesses to the temperature dependence of the affinity (or phobia) of the solvent for the surface of the particle, a factor already suggested by experimental results [35] and which possibly accounts for the change of direction of thermophoresis in protein solutions [36]. The capillary coefficient λ at a liquid–vapour interface was studied by Fisk and Widom [37] who concluded it was almost independent of temperature except for a very weak divergence close to the critical temperature. If that result can be transposed to a liquid–solid interface it means that, far from the solvent critical point, the sign and temperature dependence of the density jump ρf − ρ0 is likely to play the main role in the direction and magnitude of thermodiffusion of solid particles. Unfortunately we are not yet ready to calculate ρf (T ) − ρ0 (T ) from the knowledge of the solvent–particle intermolecular forces. To conclude that section we note that Derjaguin et al. [38] have proposed a model of thermophoresis based on a local excess enthalpy in the interfacial layer. Their local excess enthalpy is here represented by the local derivative T ∂(λ(∇ρ)2 )/∂ T . 5. Conclusions We strived to describe thermodiffusion in suspensions with thermo-mechanical tools only. Our main result is the momentum balance (50) for the particles completed by expression (53) for the interaction force. We proved with (44) and (43) that the osmotic pressure is not involved in thermodiffusion and we were led to distinguish between stress-induced and force-induced thermodiffusion, the latter prevailing for dilute suspensions. While the force-induced thermodiffusion is rather well understood for fluid-like particles, it is still mysterious for solid ones. We suggested after Anderson [30] that a non-isotropic pressure tensor in a thin liquid layer close to the particles is the key factor. And we proposed an explicit expression for that pressure tensor based on a Cahn–Hilliard description of the colloid–solvent boundary layer, focussing on the liquid density profile in the vicinity of the solid particles. We concluded that the liquid’s coefficient of thermal expansion should not be considered alone but only with the quantity ∂(ρf − ρ0 )2 /∂ T representing the temperature dependence of the small jump in liquid density between the bulk of the solvent and the interface. We proposed that the sign and magnitude of (ρf − ρ0 )∂(ρf − ρ0 )/∂ T gives the direction and the intensity of thermodiffusion for solid particles suspended in a pure fluid. Unfortunately we have not enough experimental data on the density jump to substantiate that proposal. Acknowledgement I am indebted to Eric Bringuier for drawing my attention to some issues concerning thermodiffusion in suspensions and to him and Michel Martin for stimulating discussions.

1232

D. Lhuillier / Physica A 390 (2011) 1221–1233

Appendix. Per mole vs. per unit mass Instead of chemical potentials µk per unit mass, one can be tempted to use chemical potentials µ ˜ k per molecule and to use molar fractions instead of mass fractions. This Appendix presents the main features of that second description of mixtures and suspensions of particles. Concerning thermo-statics the main difference is the use of quantities per unit volume instead of quantities per unit mass d(ρ e) = µ ˜ ⋆p dnp + µ ˜ ⋆f dnf + Td(ρ s) + ρ c (1 − c )w · dw + ρ(∂ 1g /∂ X )dX

(74)

ρ e + p − T (ρ s) = np µ ˜ ⋆p + nf µ ˜ ⋆f ,

(75)

µ ˜ ⋆k = µ ˜k +

mk 2

(vk − v)2

(76)

where nk is the number of molecules of phase k per unit volume of the mixture, mk is the molecular mass and the relation between the molecular chemical potential used here and the chemical potentials per unit mass used in the main text is obviously µ ˜ k = mk µk . The local equilibrium assumption is now expressed as T

∂(ρ e) ∂ np ∂ nf ∂w ∂ 1g ∂ X ∂(ρ s) = −µ ˜ ⋆p −µ ˜ ⋆f −J· −ρ . ∂t ∂t ∂t ∂t ∂t ∂X ∂t

(77)

Eliminating the time-derivatives with the hydrodynamic equations

∂ nk + ∇ · (nk v + jk ) = 0, (k = p, f ) ∂t ∂(ρ e) + ∇ · (ρ ev) + Π : (∇ v)s + ∇ · Q = 0 ∂t

(78) (79)

and taking thermodynamic relations (74) and (75) into account one obtains

[ T

] [ ] ∂(ρ s) ∂w + ∇ · (ρ sv) = −∇ · Q − (Π − pI ) : (∇ v)s + µ ˜ ⋆p ∇ · jp + µ ˜ ⋆f ∇ · jf − J · + (v · ∇)w ∂t ∂t [ ] ∂ 1g ∂ X −ρ + (v · ∇)X . ∂X ∂t

(80)

Noticing that mp jp = −mf jf = J

(81)

it is clear that the above entropy balance is completely equivalent to (22), hence to all what is deduced after (22). One can start with either molecular quantities or quantities per unit mass, the final results are the same. The part of the entropy production rate which depends on the gradient of chemical potentials is sometimes written as j · ∇(µ ˜p − µ ˜ f ) where j = jp = −jf . According to (81) this is possible only when mp = mf but, except for that special case, one must work with J · ∇(µp − µf ). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

S.R. de Groot, P. Mazur, Non-Equilibrium Thermodynamics, North-Holland, Amsterdam, 1962 (Chapter XI). B.U. Felderhof, Diffusion in binary mixtures and osmotic pressure gradient, J. Chem. Phys. 118 (2003) 11326–11334. D. Lhuillier, A. Ouibrahim, A thermodynamic model for dilute polymer solutions, J. Méc. 19 (1980) 725–741. A. Onuki, Dynamic equations of polymers with deformations in the semidilute regions, J. Phys. Soc. Japan 59 (1990) 3423–3426. H.C. Ottinger, Modeling complex fluids with a tensor and a scalar as structural variables, Rev. Mexicana Fís. 48 (2002) 220–229. D. Lhuillier, Internal variables and the non-equilibrium thermodynamics of colloidal suspensions, J. Non-Newton. Fluid Mech. 96 (2001) 19–30. N.J. Wagner, The Smoluchowski equation for colloidal suspensions developed and analyzed through the GENERIC formalism, J. Non-Newton. Fluid Mech. 96 (2001) 177–201. See Ref. [1] (Chapter III). L.D. Landau, E.M. Lifchitz, Fluid Mechanics, Pergamon Press, 1959 (Chapter 6). R.I. Nigmatulin, Dynamics of Multiphase Media, vol. 1, Hemisphere Publishing Corp., 1991 (Chapter 1). D. Lhuillier, Migration of rigid particles in non-Brownian viscous suspensions, Phys. Fluids 21 (2009) 023302. G.L. Hand, A theory of anisotropic fluids, J. Fluid Mech. 13 (1962) 33–46. N. Phan-Tien, Constitutive equations for concentrated suspensions in Newtonian liquids, J. Rheol. 39 (1995) 679–695. J. Goddard, A dissipative anisotropic fluid model for non-colloidal particle dispersions, J. Fluid Mech. 568 (2006) 1–17. J.J. Stickel, R.J. Phillips, R.L. Powell, A constitutive equation for microstructure and total stress in concentrated suspensions, J. Rheol. 50 (2006) 379–413. G.K. Batchelor, Brownian diffusion of particles with hydrodynamic interaction, J. Fluid Mech. 74 (1976) 1–29. J.F. Brady, The rheological behavior of concentrated colloidal dispersions, J. Chem. Phys. 99 (1993) 567–581. N.F. Carnahan, K.E. Starling, Equation of state for nonattracting rigid spheres, J. Chem. Phys. 51 (1969) 635–636. J.K.G. Dhont, Thermodiffusion of interacting colloids. I. A statistical thermodynamics approach, J. Chem. Phys. 120 (2004) 1632–1641. E. Bringuier, What Soret can learn from Seebeck, in: S. Wiegand, W. Köhler, J.K.G. Dhont (Eds.), Thermal Nonequilibrium, Jülich, Germany 2008, Godesberg, 2008, pp. 261–268. F. Brochard, P.-G. de Gennes, Effet Soret des macromolecules flexibles, C. R. Acad. Sci., Paris 293 (1981) 1025–1027. R.S. Subramanian, The Stokes force on a droplet in an unbounded fluid medium due to capillary effects, J. Fluid Mech. 153 (1985) 389–400.

D. Lhuillier / Physica A 390 (2011) 1221–1233

1233

[23] N.O. Young, J.S. Goldstein, M.J. Block, The motion of bubbles in a vertical temperature gradient, J. Fluid Mech. 6 (1959) 350–356. [24] L. Waldmann, On the motion of sherical particles in non-homogeneous gases, in: L. Talbot (Ed.), Rarefied Gas Dynamics, Academic Press, New York, 1961, pp. 323–344. [25] P.S. Epstein, Zur theorie des radiometers, Z. Phys. 54 (1929) 537–563. [26] S. Takata, Y. Sone, Flow induced around a sphere with a non-uniform surface temperature in a rarefied gas, with application to the drag and thermal force problems of a spherical particle with an arbitrary thermal conductivity, Eur. J. Mech. B Fluids 14 (1995) 487–518. [27] S. Takata, K. Aoki, Y. Sone, Thermophoresis of a sphere with a uniform temperature: numerical analysis of the Boltzmann equation for hard-sphere molecules, in: B.L. Shizgal, D.P. Weaver (Eds.), Rarefied Gas Dynamics: Theory and Simulations, AIAA, Washington DC, 1994, pp. 626–639. [28] Y. Sone, Molecular Dynamics, Theory, Techniques and Applications, Birkhauser, 2007, Secs. 4.5, 4.6 and 5.3. [29] S.N. Semenov, Mechanism of particle thermophoresis in pure solvents, Phil. Mag. 83 (2003) 2199–2208. [30] J.L. Anderson, Colloid transport by interfacial forces, Annu. Rev. Fluid Mech. 21 (1989) 61–99. [31] A. Parola, R. Piazza, Particle thermophoresis in liquids, Eur. Phys. J. E 15 (2004) 255–263. [32] J.W. Cahn, J.E. Hilliard, Free-energy of a non-uniform system. I. Interfacial free energy, J. Chem. Phys. 28 (1958) 258–267. [33] B. Widom, Surface tension and molecular correlations near the critical point, J. Chem. Phys. 43 (1965) 3892–3897. [34] H. Brenner, Elementary kinematical model of thermal diffusion in liquid and gases, Phys. Rev. E 74 (2006) 036306. [35] A. Regazzetti, M. Hoyos, M. Martin, Experimental evidence of thermophoresis of non-Brownian particles in pure liquids and estimation of their thermophoretic mobility, J. Phys. Chem. B 108 (2004) 15285–15292. [36] S. Iacopini, R. Piazza, Thermophoresis in protein solutions, Europhys. Lett. 63 (2003) 247–253. [37] S. Fisk, B. Widom, Structure and free energy of the interface between fluid phases in equilibrium near the critical point, J. Chem. Phys. 50 (1969) 3219–3227. [38] B.V. Derjaguin, N.V. Churaev, V.M. Muller, Surface Forces, Plenum, New York, 1987.