Solid/fluid interfaces in non-linear steady states

摘要

Furuduy Discuss., 1993, 95, 329-346 Solid/Fluid Interfaces in Non-linear Steady States Stephen Warr and Leslie V. Woodcock* Department of Chemical Engineering, University of Bradford,Bradford,UK BD7 1DP A molecular dynamics (MD) method for investigating the steady-state coexistence of solidifying systems is described. The model system may be applicable to the formation of non-equilibrium crystals or glasses under highly non-linear conditions such as splat quenching in external fields, or to the formation of colloidal crystals or colloidal glasses in uniaxial compaction processes. As a first step towards determining constitutive relations for these processes, computer simulations have been undertaken for the continuous steady-state uniaxial compaction of soft-sphere model particles in one, two and three dimensions.The system is characterised by the usual two intensive thermodynamic state variables, density and temperature (or a Stokes friction constant for colloidal systems), and two additional steady-state variables, a uniform driving force and an interface velocity. At a sufficient driving force these conditions lead spontaneously to, and completely characterise, coexist- ing low- and high-density uniform phases, the lsquo;feedrsquo; and the lsquo;bedrsquo;, and an interface between them. For rigid spheres, a reduced interface velocity (q,), relative to the driving force (F,)determines the properties of the non-equilibrium defect crystal or glassy solid phase. Preliminary results suggest this simple idealised steady-state system exhibits a rich phase diagram, showing a kind of phase behaviour, analogous to two-phase coexistence and critical behaviour, with interfacial characteristics reminiscent of thermodynamic equilibria.In attempting to understand and describe the process of crystallisation, although there is an extensive, and long-standing literature of experimental and theoretical little progress was made until the advent of fast computers. The early computer simulations of crystal-melt coexistence have been able to yield detailed pictures of the structure of the interface between the two phases in thermodynamic coexi~tence.~The information provided by simulation is essential for a complete microscopic description of two-phase coexistence and it is largely unobtainable by conventional experimental methods.As a result of previous simulation studies, the crystal/liquid interface of simple models such as the Lennard-Jones and soft-sphere models are reasonably well under~tood.~ This implies thaeuro; all linear crystal-growth processes are likewise understood because, in these conditions, the mechanism of crystal growth is merely an imperceptible bias on the thermal fluctuations of the equilibrium coexistence interface. More recently NEMD (non-equilibrium molecular dynamics) computer simulations have been reported for the investigation of crystal growth induced by temperature gradient^.^,^ The present simulation of steady states is a novel alternative approach for highly non-linear crystal-growth processes, or even vitreous-state growth processes.This approach offers a potential means towards the characterisation of highly non-equilibrium solid states using a continuous and well defined production process for a minimum number of fixed steady-state variables. In this preliminary report, a general model is investigated in both two and three dimensions, which represents probably the simplest imaginable non-linear solidification 329 Non-linear Steady States process. A system of particles interacting through a simple pairwise additive potential, initially in thermal equilibrium, is subjected to a uniform external force which acts to accelerate all the particles unidirectionally, and simultaneously, towards an impedance in the form of a semi-permeable periodic boundary which results in a steady-state uniaxial solidification process.The present programme of computational research has been initiated in the belief that this simple model has a general fundamental significance, especially in the context of this Discussion on crystal-growth processes. The description and determination of the constitutive rheological relations for computational fluid dynamics (CFD) simulations of a diversity of solidification processes for colloidal systems, such as filtration, sedimentation, centrifugation, tabletting, slip casting, or the reverse (desolidification) process of fluidisation of powders, requires computer simulations of the appropriate steady states in the first instance.The fundamental quasi-thermodynamic description of highly defective crystalline states of matter, including the glassy state, requires characterisation, in the same sense that two thermodynamic variables characterise an equilibrium crystal.rsquo; This is achievable for crystals produced under non-linear crystal-growth processes, or glasses, by the use of steady-state variables. A topic of much discussion in the literature of the glassy state concerns the number of parameters, additional to the requisite two intensive state thermodynamic variables, that are required to uniquely characterise the state of a glass8 From the present approach, it can be seen that, in general, two additional steady-state variables are required, but this reduces to one for simple scalable models such as hard spheres.Although a great deal is now known about crystal/liquid interfaces at equilibrium, at present, virtually nothing at all is known about the interface between a highly non- equilibrium melt and its substrate, this is a new and potentially fruitful area of research, since the overwhelming majority of interfaces in laboratory and engineering processes, are dynamical and highly non-linear, just as almost all solids that we encounter in everyday life are non-equilibrium solids.8 In the following sections a simple computational model for the continuous uniaxial solidification of an infinite periodic system of interacting particles is investigated. Despite the fundamental significance of such a steady-state system, there appear to be no previous simulation studies or theoretical studies of this or similar models in the scientific literature.The following sections summarise the salient results to date of some exploratory investigations into one such steady-state model for a system of repulsive spheres. Many of the computational details are omitted here, but further details and more extensive results are available in a the~is.~ Method: Steady-state Uniaxial Compression The present simulation of solidification processes by means of a steady-state uniaxial compression represents a departure from previous simulations which have been generally either time-dependent approaches to an equilibrium statelo or involved the formation of static solid structures by accretion of single partic1es.lrsquo; In the present method, an otherwise equilibrium fluid system of interacting particles at a prescribed mean density is allowed to solidify continuously under the influence of a uniform driving force.In addition to the interparticle forces, an external uniaxial solidification driving force (F,)is imposed in the x-direction. A sufficient force would merely cause the particles to lsquo;solidifyrsquo; at a rigid boundary, or, in a periodic boundary, to move continuously across the system with an ever-increasing velocity. A steady-state solidification is effected by the use of a hybrid boundary condition, which permits the particles to recirculate at a prescribed lsquo;feedrsquo; velocity. In practice, this is equivalent to the imposition of a pseudo-rigid boundary in the x-direction which moves at a given relative velocity to all the particles.This second steady-state system variable is called the boundary velocity (u,,). S. Warr and L. V. Woodcock 331 Standard equilibrium MD programs have been modified to simulate the continuous process of steady-state uniaxial compression under both isokinetic (isothermal) conditions using modified Newtonian equations of motion and Stokesian dynamics using a simple mean-field approximation. In the isothermal steady-state compression, the total kinetic energy of the particles is maintained constant by continuous renormalisation of their peculiar velocities.12 The Stokesian conditions relate to colloidal particles using a one-body mean-field Stokesian form for the equations of motion of the particles.In this case the Stokes friction constant replaces the total kinetic energy as an independent system-state variable. The details of both these computational schemes for the equations of motion are reported elsewhere.I2J3 Computations in all dimensions have been carried out in both the isokinetic and MF-Stokesian schemes for the equations of motion. In each case where the two sets of results for the steady-state properties were compared at the same total kinetic energy, i.e. at the same steady-state point, there were no significant differences in any of the system properties cal~ulated.~Accordingly the only results reported here were obtained using the isokinetic equations of motion.The kinetic energies are fixed apriori at 0.5 kT, kTand 1.5 kTin one, two, and three dimensions, respectively, for reference to thermodynamic systems. The boundary equations-of motion are best described by reference to the hard-sphere interaction model. The particles all move under the force F,with uniform acceleration F,/m towards the semi-permeable wall in the direction of the force. u, is the relative particle velocity (the boundary wall), uais the absolute particle velocity and 1.rsquo;b is the boundary wall velocity, i.e. U,= U, -q,. On elastic collision with the wall the absolute velocity is reversed v,(new) = -ua(old) + t)b and since u,(old) = u,(old) + q, therefore u,(new) = -v,(old) -2ub If the new relative velocity is negative, the particle rebounds from the wall, whereas if the new relative velocity is positive after the wall collision, the particle is deemed to have traversed the periodic boundary and is recycled, with the new relative velocity given by eqn.(I), as part of the lsquo;feedrsquo;. In this simple scheme, the boundary velocity determines the feed velocity, and, for a lsquo;two-phasersquo; solidifying system, is the same as the interface velocity, under steady-state conditions. For soft pair potentials, particle-wall interactions can be accommodated into a finite time-difference algorithm. Following the normal upgrade of positions and velocities, for particle-particle and particle-wall interactions, every particle in the simulation cell is translated forward in the x-direction by an amount Ax = At, where At is a time increment in the integration of the equations of motion.Once all the particles have been translated, any particle with its centre outside the semi-permeable periodic boundary of the simulation cell is recycled through the opposite side to become part of the lsquo;feedrsquo; again. Full periodic boundary conditions prevail in the directions perpendicular to the direction of uniaxial flow. Under isothermal conditions, the imposition of sufficient magnitudes for F,, at a constant q,,causes a steady-state lsquo;feedrsquo;, lsquo;bedrsquo; and interface to be spontaneously generated. The interaction pair potential between particles used in the present computations is a simple generalisation of the hard- and soft-sphere models, which is also a useful representation of the effective interaction between spherical colloidal particles and which retains some simplifying scaling properties. The pair potential is defined by13 Non-linear Steady States Writing the pair potential in the form amp;j = 10-"kT(rii/a-1)-" (3) where k is Boltzmann's constant, an effective hard-sphere diameter (aeff)can be defined as the separation of particle centres for which +g = kT, then aamp;-= a( 1 + lo-"'") (4) All the computations reported here correspond to the effective pair potential for n=m= 12.One-dimensional Systems Many phase-transition and interface phenomena in two and three dimensions do not exist in analogous one-dimensional situations, but the present steady-state phenomena are to a certain extent exceptional.Indeed, the one-dimensional steady-state compaction process has a very common analogue in everyday experience in modern life in the form of traffic jams. It is well known that a partial impedance to the smooth flow of a line of traffic, such as traffic lights or a visible speed trap, can cause the coexistence of a dense logjam at the impedance, and a rather abrupt transition to the less dense, but faster, feed line of traffic arriving at the queue. The one-dimensional steady-state particle compaction is of special academic significance. It has a completely uncomplicated solid phase and it may be amenable to analytical description. It is investigated here in the first instance as the simplest of all systems to exhibit these steady-state coexistence phenomena.Isokinetic simulations have been performed for one-dimensional soft spheres interact- ing via eqn. (2) to study the effect of the external force, F,,on the equilibrium properties of beds formed when ?+,= 0. A series of runs, with the kinetic energy fixed at 0.5 kT, and a reduced density 0.5, and for F,values ranging from 0.01 to 10 gave density, kinetic energy and pressure profiles of which many of the qualitative features of these equilibrium phenomena are also exhibited in higher dimensions, at least for the static phenomena. The pressure when vb = 0, for example (Fig. I), varies uniformly from zero at the free surface of the system to the force per unit area at the wall, the gradient being simply inversely proportional to the depth of the 'bed'.The structure at vb = 0 ranges from tightly compacted beds with sharp interfaces, at high values of Fs,as seen for F,= 1 in Fig. 1, to very diffuse interfaces at low values of F,.In the limit of very low F,,approaching zero, uniform systems are obtained, without interface, with shallow density gradients across the whole system as a minor perturbation of the unperturbed thermodynamic system. It seems that the transition from one-phase to two- phase behaviour is not sharp in this case, in contrast to the higher dimensional behaviour to be discussed in the following sections. For F,values of 0.3 and above, a definite region of 'bed' exists which has an almost constant number density profile because of the low compressibility of the 'bed'.This region increases in extent from the boundary as F, increases, until when F, = 10 all the particles are in a uniform bed and its number density (No/L)= 0.9092. The corresponding effective packing fraction is, according to eqn. (3), 1.00012! This force is clearly sufficient to effect essentially maximum packing in one dimension. The kinetic energy profile in ID fluctuates around 0.5 in the 'bed', and then rises in the interfacial region before dropping to zero. The rise in interfacial kinetic energies per particle reduces as F,decreases, but it appears to be a definite feature of the equilibrium ("b = 0) configurations. The pressure profiles are all linearly varying between zero and the total force per particle acting on the wall, which is proportional to F,.The effects of non-zero and progressively increasing values of the interface velocity q,, for isokinetic steady-state solidification is one dimension, have further been determined and are also shown in Fig. 1.These profiles show that the ID system exhibits essentially the same S. Warr and L. V. Woodcock 0.5 0.0 I = l ' l ' l ' l ' l 1 a 0.2 0.0 l ' l ' l ' l ' l ' solidification phenomenology as seen, and described in the following sections, in higher dimensions. The proportion of 'bed' at steady state decreases as ub increases. In ID, the 'bed' densities are virtually the same for all values of q,,and slightly lower than the value for which ub = 0.Both the kinetic energy and pressure profiles show the same predominant features as found in higher dimensions. The 'feed' kinetic energies are approximately proportional to the 'feed' densities, and the pressure profiles are linearly proportional to the force and with the same gradient independent oft+,. Two-dimensional Steady-state Simulations A detailed study was initially performed on a system of 504 particles at the mean reduced number density p = No2/A= 0.5. The simulation cell is rectangular in shape of area, A = 55 x 18.33~~.All the simulations were started from a uniform disordered equilibrium arrangement and were equilibrated for 40 000 time steps with At = 0.01 in reduced units of wi, 0 and kT. All the systems studied reached steady-state equilibrium within this time 334 Non-linear Steady States period, and they were then run for a further 40000 time steps to enable the sampling of steady-state lsquo;ensemble averagesrsquo;.Fig. 2 shows the equilibration period and fluctuations in the mean total pressure averaged over time for both high and low F,cases for a range of boundary velocities. The magnitude of fluctuations and the total pressure increases with F, and decreases with increases in q,.Indeed, for hard spheres, there is a corresponding-states inverse relationship between these two steady-state variables; they cannot be varied independently and a single hybrid system state variable (Fs~/m~)suffices to determine the steady-state properties of the system.The stability and homogeneity of each of the two phases and the interface between them is illustrated in Fig. 3. In order to check the sensitivity of properties in the two-phase region to system size effects, detailed results have been obtained for variations in N (the total number of particles in the system) up to N = 1000in 2D. These results generally show that, provided the system size is sufficient for the interface and the atypical boundary region, which can extend a few particle diameters into the bed, not to interact, i.e.to be separated by several particle diameters of a uniform region of lsquo;feedrsquo; and lsquo;bedrsquo;, a remarkably small number of particles suffices to gain an accurate insight into the behaviour of an infinite system. In most regions of the Fs-l)b space, the essential physics, and reasonable estimates of the properties of both phases, and the nature of the solidification interface, can be obtained for a system of as few as 200 particles. To obtain an overview of the behaviour of the 2D system, five values of F, were taken (i.e.20, 10, 5, 3, 1) for four values of l)b in the range 0-1.0. Additional state points were studied for ub = 0 over a range of F,down to a value of 0.01, whereupon F,becomes ineffective and the system remains essentially homogeneous. Fig. 4 shows a computer graphics plot of instantaneous snapshots of particle positions which illustrate the general ldquo;V -40--30--20--+w+ rsquo;010; 1 ,rsquo; ,. IlIIIII ;,,,I, 0 10 20 30 40 50 60 70 80 no.of time steps x Fig. 2 Fluctuations of the bulk pressure (pA/NkT)during the equilibration period to steady-state uniaxial solidification in two dimensions for selected state points: (a) F, = 20, z+ = 0.008; (b)F, = 3, t,j, = 0; (c) F, = 3, uh = 0.008; (6)Fs = 3, ub = 0.03 S. Warr and L. V. Woodcock 335 A 15 2 C .-...Y g 10 10 a b 5 h .g 0 0 0s 0 10 20 30 40 50 B X position 0 10 20 30 40 50 60 70 80 90 100 C X position Fig. 3 Dependence of the properties of the steady-state solidification system in 2D on the number of particles (N);the computer snapshots (with their instantaneous profiles of density and kinetic energy superimposed) are at the same state point (F,= 10,211, = 0.008) and show only a weak variation of the properties of the system for as little as a few hundred particles.The time-averaged density profiles (below) indicate the insensitivity of the lsquo;feedrsquo; and lsquo;bedrsquo; properties to system size (a) 504, (b) 1008. trend in system behaviour when the boundary velocity is zero for a range of values of F,. In this case, the effect of the force is to condense the particles on to the rigid boundary substrate with a progressive increase in the density of the condensate and an attendant progressive sharpening of the free surface. In the computer graphics lsquo;snap-shotsrsquo;, particles are shown as circles of the effective hard-sphere diameter as given by eqn. (4). The profiles of number density (upper line) and kinetic energy (lower line) are also superimposed.Fig. 5 shows that the effect of increasing ?+, at fixed F, is to decrease the proportion of lsquo;bedrsquo; in the system and increase the amount of lsquo;feedrsquo; at steady-state. In all these 2D studies the solidifying phase takes the form of a hexagonal crystalline arrangement, with few defects, and the interface is rather diffuse, extending over several particle diameters. The appearance of a layered structure in the particle lsquo;feedrsquo; at some state points appears to be caused by strings of particles being recycled; it seems likely that this may be at least partly due to the smallness of the system having an increasingly significant distorting effect in 336 Non-linear Steady States A G .c( z 10I;a x 0 X position 0 10 20 30 40 50 X position Fig.4 Computer snapshots (N = 504, p = 0.5) showing the structure of the 2D system in thermodynamic equilibrium states for a range of F,when the boundary velocity (D,,) is zero; for small values of F,,below ca. 0.1, a uniform phase with a weak density gradient fills the system, at intermediate values of the order of unity there is evidence of two distinct phases, and at high values (5) there is just the solid with the odd defect and surface roughness. F,: A, 0.1, B, 0.5, C, 1 .O, D, 5.0. d .r(Y.I I; 10 a x 0 X position A c .I . aI; 10 x 0 c .AC.- I; lo a x 0 r. X Dosition c .-c .d m a 10 x 0 C .-U.4 m 0 a 10 x 0 E X position 2 G .-.-Yr; a 10 xU.-2 x 14 0 0 0 10 20 30 40 50 X position Fig.5 Computer snapshots (N = 504, p = 0.5) taken from simulations of two-dimensional steady states showing the effect of varying the boundary velocity (u,,) at a constant external driving force (F, := 3.0); there are two clearly defined regions, homogeneous at higher u,,, and two-phase for lower values of q,between ca. 0.5 and zero. c,,:A, 0.8, B, 0.3, C, 0.03, D, 0.008, E, 0. Non-linear Steady States 1.5-1.0 -0.5 -8.0 6.0 4.0 2.0 0.0 0 10 20 30 40 50 profile Fig. 6 Profiles of density, kinetic energy and pressure for steady-state simulations of the 'solidification' of the two-dimensional system at a constant external driving force (F, = 1.Oj; the interface velocity (u,,)varies from zero to 0.08 as shown, and the profile dimension is the particle core diameter (0).vb: (a) 0, (6) 0.002, (cj 0.005, (d) 0.008, (e)0.01, (f) 0.03, 0.05, 0.08.highly non-linear conditions. Qualitative observations of these snapshots, over the whole range of system-state variables examined, shows that for all the values of F,,in the limit that l)b becomes very large, the two-phase situation reverts to a single phase again. This apparent homogeneous phase resembles the unperturbed equilibrium state in structure but, in fact, it is anisotropic in both its energy and pressure profiles. Fig. 6 shows the profiles of density, energy and pressure over the whole range of q,when F, = 10.The density profiles confirm that two regions of ub-space exist; a two-phase region at low values of l)b, and a homogeneous one-phase region at higher values of ?+,. These profiles show that there is uniformity within the separate phases, and that the extent of the interfacial region appears to diverge in a manner resembling liquid/vapour interfacial profiles approaching the critical point. The corresponding kinetic energy profiles, however, show that equipartition does not prevail in these highly non-linear steady states; the feed kinetic energy always substantially exceeds the bed kinetic energy. Directionally resolved profiles, moreover, (not presented here) show that the kinetic energy profiles become increasingly anisotropic with the solidification velocity.S. Warr and L. V. Woodcock Pressure profiles, calculated using the usual Clausius virial theorem, for this series of steady states, at F, = 10, are also shown in Fig. 6. A clear trend of decreasing lsquo;feedrsquo; reduced pressure per particle, as the lsquo;feedrsquo; density increases, is observed. The decrease occurs because the lsquo;feedrsquo; kinetic energy, at steady state, drops significantly as q,increases. The interfacial region is marked by a change in slope of the pressure profile, which subsequently rises steeply over the initial surface region of the lsquo;bedrsquo;. After its initial rise, the pressure then levels off slightly to give a profile of constant gradient, which is the same for all lsquo;bedsrsquo; formed for a particular value of Fs,regardless of the thickness.It is interesting to observe that this non-uniformity in the lsquo;bedrsquo; is essentially a system-size effect; as the system size becomes very large, this pressure gradient approaches zero. The effect of varying the driving force (F,)when l)b is fixed is illustrated for a range of F, at l)b = 0.01 in Fig. 7. This indicates a rather complex behaviour which is difficult to interpret. At low values of F,, the extent of the lsquo;bedrsquo; increases from being non-existent at F, = 0, to a maximum extent around F, = 1, with a very diffuse interface. Then, as F, increases further, the interface sharpens up and the bed decreases further as the lsquo;feedrsquo; appears to become first anisotropic with vertical striations, and eventually, at the highest F, values investigated, it becomes inhomogeneous with large void regions.This complex behaviour needs to be investigated further, but it does seem that unexpected structures in both the lsquo;feedrsquo; and the lsquo;bedrsquo; can be obtained when either F, or approach large values whilst the other remains constant. Another remarkable property of this steady-state model is its resemblance to a thermodynamic system in the two-phase coexistence region. Indeed, if the effect of the external force can be regarded as a kind of perturbation which biases the unperturbed equilibrium phase diagram in favour of condensed phases at a substrate surface. It can be conjectured that the additional effect of a finite l)b will give rise to a new generalised steady- state phase diagram with its origin in the order-disorder equilibria, and/or glass transitions, of the parent thermodynamic system.To investigate this conjecture further, a range of computations have been carried out for a single state point, but over a wide range of the mean particle density. This is equivalent to changing the system volume containing a constant amount of material. In a true thermodyamic one-component system, the Gibbs phase rule dictates that properties of the coexisting phases remain invariant; only the relative amounts of each phase can vary. The results obtained in this investigation so far suggest that very similar conditions also apply to these steady-state systems in the two-phase regions. Fig. 8 illustrates that there is a range of density over which the amounts of lsquo;feedrsquo; and lsquo;bedrsquo; vary, but their properties remain approximately constant.Results for a Three-dimensional System In three dimensions steady-state properties have been obtained for the isothermal compression, over a certain range of F, and ub space. The total kinetic energy is constrained in 3D to 1.5 kT. All the results obtained are not presented in detail here; the phenomenological trends are much the same as observed in two dimensions. In 3D, a detailed study has been performed on a system of 500 particles, also at the reduced number density p = No3/V = 0.5. The snapshot results of such a study for the state points for which the constant driving force Fs = 10 are shown in Fig. 9. The corresponding profiles of density, kinetic energy and pressure are shown in Fig.10. The most significant difference between 2D and 3D is that in 3D an amorphous solid phase is obtained for the state points presently accessible; this is almost certainly partly due to the fact that the equilibration times available for the 3D steady states may be insufficient to give the true steady states which, in the long-time limit, would probably consist of crystalline lsquo;bedsrsquo; with varying degrees of defects. By analogy with the description of the metastable supercooled A 2 c.$ 10 a v)0 .- x Y.* 1 2 x 0 0 R x Dosition 2 c .e.$: a b 10 xU.-v1 laE x 0 0 C X position :: D ''4 X position 2 c x c,.C10 ;="b a 14 x 0 0 E xDosition 0 10 20 30 40 50 X position Fig.7 Computer snapshots (N = 504, p = 0.5) taken from simulations of two-dimensional steady states showing the effect of varying the external driving force (F,)at a constant boundary velocity (vb = 0.01); at very high forces the 'feed' exhibits structural anisotropy, first perpendicular to the interface, and then parallel. F,: (a)20, (h)10, (c) 5, (d)3, (e) 1. A 20 2 e A Y.-."Y.I 6" 8a 10 14 $. 0 I I I I I I I I v I I I0 , . . . . , . . . , , , . . I . , . I . . . . , 0 10 20 30 40 50 60 70 B 20 2 c Y.- .-h Y I; 10 z a 14 x 0 0 0 10 20 30 40 50 60 0 10 20 30 40 50 0 4 8 12 16 20 24 28 32 36 40 44 48 0 4 8 12 16 20 24 28 32 36 40 44 X position Fig. 8 Computer snapshots (N = 504) taken from simulations of two-dimensional steady states showing the effect of varying the mean density whilst both the boundary velocity (ub = 0.008) and the external driving force (F, = 10.0) remain constant; there appears to be a region of density for which structures of the 'feed' and 'bed' remain uniform whilst their relative amounts vary.p: A, 0.3, B, 0.4, C, 0.5, D, 0.6, E, 0.7. Non-linear Steady States 0 10 20 0 10 20 Fig. 9 Computer snapshots (N = 500, p = 0.5, F, = 10)taken from simulations of three-dimensional steady states showing the effect of varying the boundary velocity (ub) at a constant external driving force (Fs);there are two clearly defined regions, homogeneous at higher l)b, and two-phase for lower values of V, between ca.0.02 and zero. ub: A, 0.03, B, 0.01, C, 0.008, D, 0.005, E, 0.002. fluid states of matter, the solid states in Fig. 9 and 10, for high F,and/or low 2)b state points, could possibly be regarded as metastable steady states. Fig. 9 and 10 illustrate a new phenomenon, which has many common experimental counterparts in the formation of colloidal lsquo;glassesrsquo; by uniaxial compaction processes, that of the lsquo;fluid/glass interfacersquo;. The snapshots and profiles show that these interfaces are all a few particle diameters in width, i.e. they are rather diffuse, and resemble vapour/liquid interfaces at equilibrium. The properties of both the interface and the amorphous solid beds, over the whole range of i.rsquo;b, remain essentially constant, whereas the properties of the feed vary widely as in the 2D cases.To assess the effect of system size when the total system number density is constant, steady states have been generated with N = 1024, 500, 250, 128 and 54, for the state point F, = 10, vb = 0.01. When N = 1024 an initial regular cubic lattice was constructed with (8 x 8 x 80) particles, and replicated twice in the x-direction. Similar procedures were followed for the smaller systems. Irrespective of the starting conditions, in each case the same steady-state point, or possibly metastable steady-state point, was arrived at with properties independent of the starting configuration and equilibration process. The density, pressure and kinetic-energy profiles, shown in Fig.1 1, clearly demonstrate that the lsquo;feed/ bedrsquo; interface properties are remarkably invariant with system size. The absolute amount of each phase increases with system size, but both the relative proportions and the properties of each phase are essentially the same to within the statistical uncertainties of the computations. A system size of 128 is about the minimum that can give a reasonable estimate of large-system feed/bed properties. S. Warr and L. V. Woodcock c.-x 3 a 0.1 0.1 0.0 0.0 0 10 profile 20 Fig. 10 Profiles of density, kinetic energy and pressure for steady-state simulations of the lsquo;solidificationrsquo; of the three-dimensional system at a constant external driving force (F,= 10.0); the interface velocity (q,)varies from zero to 0.03 as shown, and the profile dimension is the particle core diameter (a). u,,: (a) 0, (b)0.002, (c)0.005, (d) 0.008, (e)0.01, v 0.03.The effect of varying density across a two-phase region in 3D can be observed from the profiles in Fig. 12. As in 2D, the pictures confirm that there are two regions of phase behaviour. For densities below ca. 0.3, the system appears to be homogeneous, but at the density of ca. 0.3 and higher, the system comprises two phases and an interface between them. In the two-phase region the properties of the lsquo;feedrsquo; and the lsquo;bedrsquo; remain fairly constant over a wide range of density, but the relative amount of lsquo;bedrsquo; to lsquo;feedrsquo; increases with density. In three dimensions an amorphous solidification process is invariably obtained in the present series of simulations.The circumstances that would give rise to crystalline solidification steady states in 3D remains a matter of some conjecture; clearly, however, much longer computational equilibration times would be required. Conclusions In this preliminary report a simple steady-state model for the continuous rheological deformation of a system of particles corresponding to a general uniaxial compression, or, alternatively, a continuous solidification process, has been introduced. In addition to the Non-linear Steady States 0.9 U.-h m -o 0.6 0.3 0.0 6.0 h 4.08 v._U .5! 2.024 0.0 1 0.1 0.0o.2~ 0 10 20 profile Fig.11 Profiles of density, kinetic energy and pressure for steady-state simulations of the lsquo;solidificationrsquo; of the three-dimensional system at a constant external driving force (F, = 10.0) and interface velocity (v,, = 0.01) for a range of system sizes for N = 54-1024; although the relative amounts of lsquo;feedrsquo; and lsquo;bedrsquo; vary somewhat for small systems, their structural properties remain approximately constant. N: (a)54, (b)128, (c) 250, (6) 500, (e)1024. two intensive thermodynamic system-state parameters of mean density and mean temperature, two additional state variables will generally characterise this solidification process, uniquely defining the properties of the low-density lsquo;feedrsquo; and higher density solid lsquo;bedrsquo; phase and the interface between them.The immediate objective of this approach is to gain insights into the underlying physics of solidification processes in which crystalline and glassy states of matter are produced in highly non-linear condition in external fields or driving forces. It seems likely that the same physics underpins similar phenomena with colloidal suspensions in common processes such as sedimentation, filtration, centrifugation, fluidisation (the reverse of solidification), slip casting and tabletting. In these processes, the particles move in hydrodynamic flow fields within a fluid media and the particles follow Stokesian dynamical equations of motion. These results, however, suggest that provided these colloidal systems are compared with the present isothermal systems at the same lsquo;corresponding steady-states, i.e.at the same granular lsquo;temperaturesrsquo; and osmotic pressures, then there is no new physics in the S. Warr and L. V. Woodcock 1.2 I I 0.9 0.6 0.3 0.0 9.0 6.0 3.0 I0.0 1 I I 0.20 1 0.15 1 0 5 10 15 profile Fig. 12 Profiles of density, kinetic energy and pressure for steady-state simulations of the lsquo;solidificationrsquo; of the three-dimensional system at a constant external driving force (Fs= 10.0) and interface velocity (v, = 0.08) for a range of mean system densities for N = 500; although the relative amounts of lsquo;feedrsquo; and lsquo;bedrsquo; vary, the density of each phase remains approximately constant in the lsquo;two-phasersquo; region. p: (a) 0.2, (b)0.3, (c) 0.4, (d)0.5, (e) 0.7.solidification processes. Isothermal systems may be related to Stokesian colloidal suspension systems by simple laws of corresponding states. l3 An alternative way of regarding these isokinetic simulations is that they allow many steady-state points to be reached which cannot be investigated by any other means, either computationally or experimentally. When ub = 0, for example, sharp lsquo;bedrsquo;/void interfaces cannot be generated without the use of an isokinetic constraint. In the isokinetic system, the particles may be compacted by increasing F,at constant total kinetic energy. When F,is increased under non-isokinetic conditions a rise in the system kinetic energy occurs. Since 1he ratio of external force to kinetic energy is fixed at steady state, the corresponding steady state may not be reached.Thus, relaxation times for well defined beds to form may approach the infinite in non-isokinetic systems, including Stokesian systems. Isokinetic conditions may be essential, therefore, to efficiently explore the whole of the steady-state space. Thus, although there may be no direct experimental counterpart, the isokinetic models 346 Non-linear Steady States do have a far-reaching importance for the underlying physics of solidifying systems. One of the greatest assets of computer simulations is that one is not restricted to the experimental conditions of the real laboratory. Not only are truly steady states virtually impossible to achieve in the laboratory, but they cannot in any real circumstances exist at prescribed granular temperatures, as we examine atomic or molecular systems under isothermal conditions, for example.Likewise, under these highly non-linear solidification conditions it would be impossible to devise experiments to maintain the isokinetic conditions that characterise the non-equilibrium solid states, as we see is now possible in these simulations. The longer-term objective from a computational science point of view is to determine the necessary constitutive relations (material properties) that are prerequisite for a fluid mechanics CFD-type simulation of any of these solidification processes. In a real experiment the solidification is generally time-dependent with the conditions prevailing on both sides of the interface changing with time.In a typical finite element CFD simulation of such a solidification process, call would have to be made upon all the non-linear transport coefficients of the solidification process as a function of both the system properties and system-state variables over the entire relevant range. In practice, for example, this means obtaining the tensors of pressure and energy, and the corresponding diffusivity and thermal conductivities, all of which may be highly anisotropic tensors, as a function of the values of density, temperature, boundary velocity and driving force, for both the lsquo;feedrsquo; and lsquo;bedrsquo;. This still constitutes an awesome task for even the very simplest of model materials such as hard or soft spheres.We thank the SERC for the award of a Research Studentship (to S.W.).We are also grateful to Dr. Geoff Maitland and his colleagues at Schlumberger Cambridge Research for sponsorship of this project and helpful discussions. References 1 See, e.g. N. F. Mott, Discuss. Faraday SOC.,1949, 5, 11. 2 D. P. Woodruff, The Solid-Liquid Interface, Cambridge University Press, Cambridge, 1973, Ch. 3. 3 A. J. C. Ladd and L. V. Woodcock, J. Phys. C, 1978, 11, 3565. 4 J. Q. Broughton and G. H. Gilmer, J. Chem. Phys., 1983,79, 5095. 5 U. Landeman, C. L. Cleveland and C. S. Brown, Phys. Rev. Lett., 1980,452032; see also V. Landeman, W. D. Luedtke, M. W. Ribarsky, R. N. Barrett and C. L. Cleveland, Phys. Rev. B, 1988, 37, 4637. 6 M. Schneider, I. K. Schuller and A. Rahman, Phys. Rev. B, 1987,36, 1340. 7 S. R. Elliot, Physics qf Amorphous Materials, Longman, London, 1983. 8 A. J. Kovacs, Ann. NY Acad. Sci., 1981, 371, 38. 9 S. Warr, Ph.D. Thesis, University of Bradford, 1992. 10 G. C. Ansell and E. Dickinson, J. Chem. Phys., 1986,854079. 11 P. Meakin and R. Jullien, in Physics of Granular Media, ed. D. Bideau and J. Dodds, Nova Science, New York, 1992, p. 323. 12 L. V. Woodcock, CCP5 Quarterly Newsletter, 1987, No. 24, 29 (SERC Daresbury Laboratory, UK). 13 L. V. Woodcock, Mol. Simul., 1989,2, 253. Paper 3/001445; Received 8th January, 1993