Modelling and asymptotic analysis of particle-interface interaction

MATHEMATICAL COMPUTER MODELLING Mathematical

PERGAMON

and Computer

Modelling

36 (2002)

147-156 www.elsevier.com/locate/mcm

Modelling and Asymptotic Analysis of Particle-Interface Interaction L. HADJI Mathematics

Department,

The University AL 35487-0350 [email protected]

of Alabama

Tuscaloosa,

(Received September 2000; accepted February 2002) Abstract-An asymptotic study of the interaction between a foreign particle and a solidifying interface is undertaken. The analysis focuses primarily on the influence of the disjoining pressure on particle pushing or engulfment. The analysis considers only thermal effects in a pure substance in which a neutrally buoyant spherical particle, whose thermal conductivity is equal to that of the melt, is positioned near the solid-liquid interface. The interface equilibrium temperature includes the undercooling effects due to both the front curvature (the Gibbs-Thomson effect) and the long-range intermolecular forces (nonretarded vanderwaals interactions) in the thin film behind the particle. A uniformly valid asymptotic representation for the front shape is derived and used to infer the gap width separating the particle from the solid front. The latter is utilized to calculate the growth rate that results from a balance of the van der Waals and Stokes forces. An attempt is also made at deriving an expression for the critical velocity for continuous pushing from the balance of forces. @ 2002 Elsevier Science Ltd. All rights reserved.

Keywords-Solidification, Materials processing flows, Particulate-reinforced posites, Solid-liquid transitions, Stokesian dynamics.

metal-based

com-

1. INTRODUCTION The interaction

between

an advancing

solid-liquid

interface

and solid particles

is a phenomenon

which occurs in several industrial and natural processes. This is the case, for instance, of the growth of Y123 superconductors [l], the freezing of water in the presence of soil particles during the formation of frost heave in the winter season [2], the cryopreservation of biological cell suspensions [3], and the solidification processing of particulate reinforced metal matrix composites [4]. It is well known that when a particle interacts with a solidifying interface, it is either pushed or engulfed in the solid. Several experimental and theoretical studies have been performed for the purpose of quantifying the pushing/engulfment phenomena [5-111. These studies have demonstrated the existence of a critical value of the growth rate which separates pushing from engulfment. The dependence of this critical velocity on the physical and processing parameters has not been fully resolved and remains an intensive area of research. This paper considers the influence of the disjoining pressure on particle’s engulfment. We proceed by allowing one neutrally buoyant and insoluble spherical particle to be near enough This research was supported by a grant from the National Science Foundation (DMS-9700380). 0895-7177/02/$ - see front matter PII: SO895-7177(02)00111-5

@ 2002 Elsevier

Science

Ltd. All rights reserved.

Typeset

by &S-‘&X

148

L. HADJI

the interface that the gap between the planar front and the nearest point on the sphere is much smaller than the radius of the particle. This situation arises when the particle is in near-contact with the interface. Our model consists of the equations governing the solidification of a pure substance

in the presence of an inert spherical particle in the melt ahead to the interface.

The

model accounts for (i) thermal diffusion in the melt, solid, and particle, (ii) conservation of energy at the solid-liquid interface, and (iii) the modification of the local melting temperature ing pressure effects.

due to the Gibbs-Thomson

and disjoin-

The flow field is uncoupled from the thermal field due to the small value of the particle’s Reynolds number [12,13].

Thus, the equations for the conservation

of momentum in the melt, as well as

the convection term in the energy equation may be ignored in the formulation (to the order at which the present calculations are carried out of course). The determination of the profile of the solid-liquid interface,

being an equilibrium situation,

is not affected by fluid motion.

The

effect of the fluid flow does appear, however, as additional terms in the boundary condition for the interface temperature (see [14]). These two terms owe their existence to the hydrodynamic pressure and the disjoining pressure in the film of melt between the particle and the solid front. The former is much smaller in magnitude than the latter and is not included in the analysis (see [14]). Because our concern here is with the effect of the disjoining pressure on the front morphology, the appropriate length scale is selected by equating the pressure variations induced by the particle’s presence to those due to the front curvature. A nondimensionalization procedure is then undertaken to reduce the large number of physical and processing variables to only a few dimensionless numbers. Two of these groupings stand out in the analysis, namely, a disjoining parameter fl and a front curvature parameter I’. Using a small-gap asymptotic analysis [15,16], we derive a uniformly valid expression for the front shape from which we derive an expression for gap thickness as function of the relevant physical and processing parameters. Then the forces acting on the displaced particle are investigated. It is these forces that will determine the pushing or capture of the particle by the moving solidifying interface.

These forces, that depend on the

gap separation, consist of the viscous drag force, which opposes the movement of the particle, and the van der Waals force which pushes the particle away from the interface (assuming a positive disjoining pressure). The application of the steady-state force equilibrium condition then yields an expression for the growth rate.

2. MATHEMATICAL

MODEL

We consider the directional solidification of a pure substance in such a way that the solidifying front is moving vertically

upward with velocity V. An inert and neutrally buoyant spherical particle is located at a distance h, from the planar solid-liquid interface. See Figure 1 for a schematic diagram. The system is completely described by the heat conduction equations in the liquid, solid, and particle. In addition, the equation for the mass conservation is implicitly accounted for. The latter is not shown in the formulation given that its effect enters the problem only through the boundary condition, i.e., it gives rise to an extra undercooling term in the condition for the local equilibrium temperature (see [10,14]). The influence of this term on the local melting point is found to be much smaller than the other undercoolings and, therefore, is dropped from the analysis. We consider an axisymmetric geometry with a moving cylindrical coordinate system which immobilizes the solid-liquid interface. With the vertical coordinate denoted by Z* and the radial coordinate, which is taken along the planar solid-liquid interface, by r*, the governing equations for the dimensional variables (denoted by *) are dT* --v&$ at*

(r*g)

-t-K-

@T*

az*2 ’

Modelling

and Asymptotic

Analysis

149

Liquid phase Particle

1”

h;

ace

aterf

Solid phase Figure 1. Sketch of the geometry of the problem. A particle of radius a is positioned at a distance h, from the planar solid-liquid interface. The deformable interface moves upward with a prescribed velocity V.

dT.9’ -at*

i a aTs* a2T, * = KS-r*+r;sr+ dr* ( dr* > a2e2 ’ aTp dT’ a2T, i a _-_V__LEK _+ QdZ*23 at* az* ‘r* dr*

where

vaTs’ a2

T’ and ts are the temperature

subscripts

S and p pertain

diffusion

boundary,

T* = T;

coefficient

in the solid phase and particle,

z* = h’ + a - dm and of the heat flux yield

At the particle-liquid temperature

and the thermal

to properties

and

in the liquid, the continuity

of the

n = 0,

k and k, are the thermal conductivities of the melt and denotes the normal to the particle-melt boundary that is pointing

where

and the

respectively.

E B*(r*),

V (kT* - kpTi)

(2)

(4)

the particle, respectively, n into the liquid and V is the

gradient vector. Furthermore, we allow Ti to have some prescribed value To* at some elevation z* = d*. If we neglect the kinetic and hydrodynamic effects, the solid-liquid interface equilibrium temperature

is given by [lo] T = Ts = T,,, i- AT,,,,

where AT,,,,

is the curvature

undercooling

interface

energy,

(5)

given by the Gibbs-Thomson

ATcur, = g [ where Q is the crystal-melt

+ AT,,,

AH,

f

formula

6

(6)

1

is the molar entropy

of fusion,

R denotes

the

molar volume of the liquid phase (taken to be the same as the molar volume in the solid phase) and K is the mean curvature of the interface (positive when the center of curvature is in the liquid phase). The undercooling term AT,, represents the change in the melting temperature due to the vander

Waals forces in the film separating

the particle

from the interface

and is given

by 1141

ATdp = [&][67rlc*)] ’

(7)

where g(r*) is the gap thickness, i.e., the depth of the melt film that is sandwiched between the particle and the solid-liquid interface and where A is the Hamaker constant. The Hamaker constant is taken here to be positive when the disjoining pressure decreases as g(r*) decreases. Far away from the interface, the temperature gradient in the solid phase is maintained at G;. When the particle enters the near-contact region, it induces pressure variations in the melt film. These pressure changes owe their existence to the disjoining pressure PDp, and the pressure

150

L. HADJI

that is associated spherical particle

with the front curvature of radius a, the pressure

while the disjoining

pressure

Since we are concerned by setting

surface

-A

DP=

with the interface

energy

]A(/67re3 = 2a/a.

For a

is given by P

of interfacial

as a result of the Gibbs-Thomson effect PGT, variation induced by the front deformation is

6rg3(r*)

’

deformations

and the disjoining

pressure,

that

are induced

the relevant

by the coupled

length

effects

scale is determined

We obtain

al-4

e =

1’3

(10)

127ra L-1

In the nondimensionalization of the equations and boundary conditions, we use e as lengthscale in the vertical direction and the radius of the particle, a, as a lengthscale in the horizontal direction, C2/K as a scale for time and the fusion temperature of the pure substance of temperature. The following equations in dimensionless form are obtained:

T, as a scale

(11) (12) (13) where E = e/a and v = N/K. The corresponding boundary conditions in dimensionless form are discussed in the following. Since we are focusing on the liquid region that is sandwiched between the particle and the solidliquid

interface,

the upper

boundary

for the liquid region is given by (see [16]) cz=?+l-&7,

(14)

where h, is the distance between the planar interface and the nearest point on the particle, i.e., the minimum clearance between the particle and the solid front prior to the deformation of the interface (see Figure 1). The particle is assumed to be in the near-contact region so that h, << a. We now assume that h, = O(e) and so h,/a = EH where H = O(1). At this upper boundary, which in dimensionless form is denoted the temperature imply V(T - uTp). n = 0 where v is the ratio of the thermal conductivity The first of equations (15) translates into 8B aT dT __-_2__=y a.2 ar ar Let q(r,t)

be the deviation

aT, -_E2_P az

of the interface

by B(r), the continuity Tp = T,

and

of the particle

dt3 6’T dr ar 1

from planarity,

of the heat flux and of

,

to that

(15) of the liquid

on z = B(r).

then K in equation

(V = k,/k).

(16) (6) is given by

(17)

Modelling and Asymptotic Analysis where the subscript perature

T denotes

at the interface

differentiation

z = n(r,t)

with respect

(equation

to r. In dimensionless

1+c2p2 1+&r3b97 717-T

where

cd

r and p describe

-AS2

P=

r = !AH,T,

(19)

6d3AHfTm.

the influence

of the front curvature

and the pressure

forth by the molecular forces on the local melting point, respectively. We proceed with an estimation of the order of magnitude of the various rameters

that

A M 10-15erg acteristic

emerged (lerg

length

form, the tem-

(5)) satisfies

T=Ts=l+

The two parameters

151

from the analysis.

We consider

= 10P7 J), (Y N”50erg/cm2

a typical

and particle

physical

radius

dimensionless

situation

pa-

[10,14] with

a M 10F3 cm. Then

scale is e M 10V7 cm which yields E M 10m4. For V z 10e4 cm/s,

put

the char-

K z 10-l cm2/s,

we obtain I’ M .13, ,8 z 1.6 x 10Y4, and v x 10-i’. Hence, and CllAHfTm z 3 x 10-gcm3/erg, appropriate scalings for the dimensionless parameters are w = ~~6, ,0 = eb, and I = O(1) where the hat symbol

denotes

O(1) quantities.

3. FRONT In this section,

we follow the analysis

of Cox and Brenner [17] to conduct an asymptotic analysis in the vicinity of the origin. For this purpose, we introduce

of equations (ll)-(13) and (15)-(18) the inner variable r = &R. This scaling system

of equations

and boundary dT -at

SHAPE

is dictated

conditions

by the form of equation

(14). The following

results:

dT 1 d ~‘6~ =E-Rc3R

(20)

(21) (22) The upper

boundary

is now given by R2 R4 z=H+~+E~+...=D(R).

On z = B(R), th e continuity

of temperature

dT dB dT ---_--~~ dR dR a.% The temperature

and of the heat flux yields

aT, dB c3T -__E_-z dR dR I [ dz

in the solid and liquid phases T=Ts=l+

(23)

at the interface

and

T = Tp.

(24)

z = q(R, t) is given by

VRR + cBg-3(R). 1 + cr]R2 I

(25)

Since the aim of this study is analyzing the effect of the disjoining pressure, we set the thermal conductance ratio Y to unity. This corresponds to a particle having the same thermal conductivity as the melt. Therefore, as far as thermal effects are concerned, the particle is indistinguishable from the melt. We solve the steady-state version of the above problem by expanding the independent variables, including the interface, in E as

L. HADJI

152

The reduced

problem

(E = 0) pertains

to the case of a planar

interface

and is described

by

72 +1, Tp = F(* - d) + To,

onH+2-
(23)

TS = Gsz + 1,

z I 0,

(29)

T(z) =

where AT = To - 1. Here AT/d liquid film and a portion profiles

because

The particle’s

the effect of the disjoining

The ratio of the temperature

pressure

H is the depth

gradient presence

difference

across a layer consisting at O(1) but enters

the temperature

AT over the depth

of the

does not affect the temperature

does not appear

i.e., O(E). Consequently,

AT -= d where

R2

is simply the temperature

of the particle.

next order in the calculations, the vertical coordinate.

(27)

profiles

d can be written

at the

are all linear as

AT, f AT, L+H

of the melt layer under

(30)

’ L = d - H, AT,

the particle,

in

and AT,

are the

temperature changes over the depths H and L, respectively (see Figure 1). In to the continuity of the heat flux at the particle-melt interface the O(1) contribution

dimensionless addition, implies

AT, _ AT, -E---H. By combining

equations

(30) and (31), we obtain d=H

Proceeding

to the O(E) problem,

we find that

[

(31)

an expression

I+%.

for the parameter

d, namely,

m 1

(32)

the temperature

perturbation

0 is given by

(33)

Application

of the heat energy

balance

at the solid-liquid

V[KT,-T].n=C

V+z [

interface

implies

1

arl

,

where K = ks/kl denotes the ratio of the thermal conductivity of the solid phase (ks) to that of the liquid phase (kL), C is the nondimensional latent heat of fusion, C = AHftcL/kLT,,, (of growth rate. The leading order term of the numerical order unity) and V is the dimensionless latter is calculated from equation (34) to get

The last term in equation (35) follows from equation (32). This expression for the growth velocity is similar to the one corresponding to the directional solidification stance without the particle’s presence with AT,/H being the thermal gradient in effect of the particle does not appear at this order of the calculations. At O(E) of we obtain an evolution equation for the interface perturbation q(R, t), namely,

dimensionless of a pure subthe melt. The equation (34),

(36)

Modelling

If we make the change of variable that

is bounded

and Asymptotic

z = Rv’m,

at x = 0 and satisfies

kind,

and Ko(x)

la(x)

the symmetry

respectively.

are the modified

Upon using

where The

the positive parameter

parameter

y depicts

condition

dy - Ko(x)

Bessel functions

2

(1 +yy2j3

y = I’/2AT,

the competition

$J

s s=

of equation

is introduced

dY

to simplify

the interfacial

1’ (37)

= AT, (37) reduces

(1 +-fY2)3

0

”

zero of the first and second

Ko(X)Ylb(Yl)

-

between

of order

(36),

= 0, is given by

YIoo(Y) o (N + (I’~/~AT)Y~)~

AT,, + AT,

cx)~oo(X)YKo(Y) dY

[S

solution

I

(32) and the fact that

ALH2AT,

q(x) =

153

then the steady-state

YKO(Y) O” (H + (I’~/~AT)Y~)~ where

Analysis

to

1

(38)

’

the expression

and thermal

effects.

for v(x). A plot

of the interface shape, represented by equation (38), is shown in Figure 2. The front profile is convex for p > 0 and concave otherwise. Far away from the origin, i.e., away from the particle’s position,

the interface

interface

shape,

dimensionless

approaches

equation

the planar

(38), takes

profile

the following

rapidly.

Note that

form when

written

the equation

in terms

for the

of the original

variables:

-aa”-------------’ 2 4

6

X

Figure 2. A plot of the scaled interface deflection, H2AZ’,q(z)/j~I, for three values of the parameter y: y = 0.05 (dash-dotted), y = 0.2 (dotted), and y = 2 (continuous line).

4. GROWTH

RATE

FOR CONTINUOUS

Now that the profile of the solid-liquid

interface

is determined,

REPULSION

we proceed with the derivation

of

an asymptotic representation for the gap thickness g(z). Upon using the expressions for the meltparticle boundary and the solid-liquid interface given by equations (23) and (38), respectively, we obtain g(x) = B(x) - q(x) If we neglect

any buoyancy

change

= H + 2

due to the thermal

(~Hx)~ ___ 8 expansion

_ rl(x)

1

(40)

of the melt then the particle,

assumed neutrally buoyant, is subject to two main forces, namely, the van der Waals force Fvdw, which is either compressive or repulsive, and the Stokes force Fvisc which opposes the particle’s movement [16]. These forces are given by

154

L. HADJI

where p is the melt’s dimensional

dynamic

gap separation.

V is the velocity

viscosity, Continuous

pushing

of the solidifying

of the particle

Fvdw = Fvisc. Upon using equation (41), the balance for the growth velocity resulting from the force balance:

provided equation

front,

by the moving

and g* is the

interface

of forces yields

occurs

the following

(42) of V,, on the order unity

where the dependence

on H. In the following,

the gap thickness

free parameter

H comes from the dependence

of

for V,, by making

use

we look for a simplified

solution

of the value of the gap thickness at the origin. This approximation will not alter the functional dependence of V’, on the physical variables which is one of the aims of this study. Moreover, given that major units),

the interface

contribution

profile decays

rapidly

to zero (the planar

to the forces comes from near the origin.

we make use of Brenner’s

analytical

solution

interface),

we expect

that

the

For g*(O) < 1 (in some appropriate

[16] for the Stokes force

6npVa2

Fviscx ~ while the disjoining

pressure

(43)

’

force is given by [18]

F Upon equating

g*(O)

Ma

(44)

vdw = 6 (g* (0))2 ’

the two forces, we obtain

I-4

Kq(H) = 36rpae[H - q(O)]

’

where we have used the fact that g*(O) = ![H - q(O)] with C being defined asymptotic representation for V,, is, thus, given by

Kq(W

where Z(y)(=

(10). An

JAI 1

- -36rpae H

n(O)) is given by Z(Y) =

Note that

by equation

the function

Z(y) satisfies

O” YKo(Y) I0

(47)

(I + YY2)3’

0 5 Z(y) 5 1, Z(y) + 1 as y -+ 0 and Z(y) -+ 0 sharply

as

y --f co. The O(E) contribution that appears in equation (46) represents the influence of the interfacial deformation on V,,. This contribution is shown to depend on the disjoining pressure parameter 6, the temperature difference across the melt film AT,, the parameter y = I’/2AT,, and the dimensionless thickness H. Thus, one may undertake a qualitative analysis of the effect of each one of these parameters on V,, directly from equation (46). Upon using the facts that efi = p and that the parameter H was related to the dimensional gap thickness by H = h,/e, we write equation (46) as (keeping in mind that it is valid only asymptotically)

(48) In our expression decreasing

for Veq, the dependence

with H because

on the gap thickness

the term (&?(y)/ATmH3)

H is such that

< 1 and H = O(1). Thus,

V,, is monotonic Veq(H) defined

Modelling

by equation

(46) would attain

the asymptotic

constraint

we may postulate

that

and Asymptotic

its maximum

that

value at the smallest

H be order unity.

the maximum

Analysis

155

value of H subject

possible

From the asymptotic

of V,, is attained

analysis

point

at H = 1, which corresponds

to

of view,

to h,

= C,

to obtain K+&[1+~]. From

a physical

defined

(49)

to the critical distance of view, h, = constant x e. This corresponds et al. [4]. This critical distance is determined and Rogge [ll] and Shanguan

point

by Pijtschke

by the range of the van der Waals forces and depends on the particular Upon using equations (10) and (19) for JJ and P, respectively, we obtain v =

(B12’3(241’3

c

where AT, depicts

* is the dimensional

36rpa4/3 temperature

difference

of V, on the physical

the dependence

2fwf)Q aAHfTAAT,*

experimental

1’

(50)

across the melt film. Equation

parameters,

situation.

is the main motivation

(50), which

of the foregoing

analysis.

5. CONCLUSION A new approach and a solidifying

was undertaken interface.

ods, is used to investigate

to investigate

The approach, the influence

the interaction

between

which is based on the application of the disjoining

pressure

an insoluble

particle

of asymptotic

on particle

engulfment.

methAn

insoluble and neutrally buoyant spherical particle is assumed to be in the near-contact region. The latter pertains to a situation when the gap separation between the particle and the planar solid-liquid interface is much smaller than the particle’s radius. In such a situation, length scale, which is determined from a balance between the pressures induced

the relevant by the front

curvature, as a result of the Gibbs-Thomson effect, and the disjoining pressure in the melt film behind the particle, is given by equation (10). The nondimensionalization of the equations of thermal conduction in the melt, solid, and particle, which incorporate the moving interface, and corresponding solidification boundary conditions give rise to several dimensionless parameters. Of those, only two stand out in the analysis, namely, the disjoining pressure parameter p and the front curvature parameter P. Using a small-gap asymptotic analysis, we have derived an expression for the front shape presence

(see equation

to the disjoining

which is uniformly (39)). The amount

pressure

parameter

valid for r > 0 and which accounts of deflection

/? but is inversely

of the solid-liquid proportional

interface

for the particle’s is proportional

to the temperature

difference

across the melt film AT,. The interface profile is utilized to derive an asymptotic expression for the gap width. However, only the value of the latter at the origin is used in evaluating the leading order terms of the hydrodynamic and disjoining pressure forces acting on the particle. The resultant of these two forces yield an asymptotic representation for the equilibrium velocity (see equation (45)). A sort of maximization procedure is then adopted to find an expression for V,. Given that V,. decreases with H, we postulated that the maximum velocity is attained at the smallest possible H such that H is order unity. dependence of V, on the various parameters.

This led to equation

(50) which shows the

REFERENCES 1. A. Endo, H.S. Chauhan, T. Egi and Y. Shiohara, Macro-segregation of YzBaiCuiOa particles in Yi BasCus 07-s crystals grown by undercooling method, J. Mater. Sci. 11, 795, (1996). 2. K.A. Jackson and B. Chalmers, Freezing liquids in porous media with special reference to frost heave in soils, J. Appl. Phys. 29, 1178, (1958). 3. C. Korber, G. Rau, M.D. Cosman and E.G. Cravalho, Interaction of particles and a moving ice-liquid interface, J. Crystal Growth 72, 649, (1985).

156

L. HADJI

4. D. Shangguan, S. Ahuja and D.M. Stefanescu, An analytical model for the interaction between an insoluble particle and an advancing solid-liquid interface, Met. nans. A 23, 669, (1992). 5. D.R. Uhlmann, B. Chalmers and K.A. Jackson, Interaction between particles and a solid-liquid interface, J. Appl. Phys. 35, 2986-2993, (1964). 6. G.F. Boiling and J. Cisse, A theory for the interaction of particles with solidifying front, J. Crystal Growth 10, 56-66, (1971). 7. J. Cisse and G.F. Bolhng, Study of trapping and rejection of insoluble particles during freezing of water, J. Crystal Grou~th 10, 67, (1971). 8. .I. Cisse and G.F. Bolling, Steady state rejection of insoluble particles by salol grown from the melt, J. Crystal Growth 11, 25,(1971). 9. S.N. Omenyi and A.W. Neumann, Thermodynamic aspects of particle engulfment by solidifying melts, J. Appl. Phys. 47, 3956-3962, (1976). 10. A.A. Chernov, D.E. Temkin and A.M. Mel’nikova, Theory of the capture of solid inclusions during the growth of crystals from the melt, Son. Phys. Crystallogr. 21, 369-373, (1976). 11. J. Potschke and V. Rogge, On the behaviour of foreign particles at an advancing solid-liquid interface, J. Crystal Growth 94, 726-736, (1989). 12. L. Hadji, Asymptotic analysis of the crystal-melt interface profile near a foreign particle, Europhysics L&t. 51, 413-420, (2000). 13. L. Hadji, State of the Art in Cast Metal Matrix Composites in the Next Millennium, (Edited by P.K. Rohatgi), pp. 101-114, Warrendale, PA, (2000). 14. P. Gasses and M.A. Azouni, Thermal effects on the shape of a solidifying interface near a foreign particle, J. Crystal Growth 130, 13-20, (1993). 15. M.E. O’Neill and K. Stewartson, On the slow motion of a sphere parallel to a nearby plane wall, J. Fluid Mech. 27, 705-724, (1967). 16. H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface, Chem. Eng. Science 16, 242-251, (1961). 17. R.G. Cox and H. Brenner, Chem. Engng. Sci. 22, 1753, (1967). 18. E.J. Verwey and J.Th.G. Overbeek, The Theory and Stability of Lyophobic Colloids, Volume 18, pp. 239-241, Elsevier, New York, (1948).