J . Chem. SOC., Faraday Trans. I , 1985, 81, 919-937 Capillary Phenomena Part 24.-Properties of Fluid Bridges between Solid Rods in a Gravitational Field BY ERNEST A. BOUCHER* School of Chemistry and Molecular Sciences, University of Sussex, Brighton BNl 9QJ AND MICHAEL J. B. EVANS Department of Chemistry, Royal Military College of Canada, Kingston, Ontario, Canada Received 22nd May, 1984 Properties are predicted in dimensionless terms for fluid bridges between vertically aligned rods in a gravitational field. Meridian curves, which on rotation give the bridges, have been obtained for reduced rod radii R = 0.5-10.0. Actual rod radii are given by r = aR, where the capillary constant a is ca. 3.8 mm for water bridges in air and 8 mm for silicon at its freezing point in an inert atmosphere subjected to terrestrial gravity.Properties such as the pressure difference across the bridge interface and the force on the upper rod have been obtained for cases of varying fluid volume at several fixed rod-separation distances and of varying rod separation at fixed volumes. Stability limits as well as equilibrium paths are discussed. The results also provide a basis for analysing ideal float-zone refining, e.g. a silicon zone formed from its own solid: in these cases there are limitations on the fluid volume and on the angle of inclination of the fluid interface needed for solidification of a rod of purified material of uniform radius. Taking a growth angle of 11" for silicon, the maximum rod radius that can be grown from a zone of melt is ca.40 mm. The purpose of this paper is to extend the study of the properties of fluid bridges between equal solid rods with planar ends subjected to a gravitational field. An account of the basic features has already been given:' further work has been carried out to cover a wider range of bridges up to limits of stability. The approach is mainly based on the accurate computation of meridian curves meeting the solid rods at their edges as indicated in fig. 1. In contrast, fluid bridges meeting planar horizontal solids of large extent at a definite contact angle have been discussed separately.2 While the present study deals with the fundamental properties of bridges between rigid rods, the systems act as idealized models for the behaviour of float-zone processes used for the purification of silicon to semiconductor q ~ a l i t y .~ In much the same way, the holm/rod systems already studied43 are models for Czochralski growth from the melt, whereby crystals are pulled from the melt contained in a crucible. The rupture of fluid bridges has been discussed by Wolfram.6 Coriell et al.7 considered the stability of a fluid zone in a gravitational field, paying particular attention to the case where the volume of the molten fluid zone equals that of the solid cylinder from which it formed. The preliminary molten-zone studies of Keck et aL8 are not completely reliable because they did not use the full range of possible meridian curves (bridge shapes). Furthermore the often-cited study of Heywang9 is limited to an approximate account of special cases. ForteslO has considered fluid bridges in the absence of a gravitational field and reached conclusions compatible with our earlier 919920 FLUID BRIDGES Fig.1. Generating R curves which on rotation give a basic fluid bridge of phase a between the solid discs of radius R and surrounded by phase p. The vertical direction is z and the horizontal x axis coincides with the top of the lower solid. The forcefapplied to the upper solid is indicated on the small bridge. acc0unt.l The differences and similarities in properties of systems in the presence and absence of a gravitational field have also been discussed for a fluid ring round a solid sphere contacting a horizontal plane.ll COMPUTED QUANTITIES BASIC EQUATIONS The results given in this paper are based on meridian curves meeting the edges of the two solids of equal radius as shown in fig.1. [See also fig. 1 of ref. (l).] Reduced or dimensionless quantities, usually signified by upper-case letters, are used such that linear dimension I is scaled by the capillary constant a = (2y/Apgp, e.g. in reduced terms L = I/a, where y is the fluid/fluid interfacial tension, Ap is the positive fluid density difference and g is the acceleration due to gravity. The equations to be solved for the meridian coordinates ( X , 2) when the arc length S is the independent variable are d@/dS = 2(H - 2) - sin @ / X (1) dX/dS = cos @, dZ/dS = sin Q,. ( 2 ) The meridian angle Q, is the angle of inclination of the meridian at ( X , Z ) to the horizontal, i.e. tan@ = dZ/dX.The shape factor is H = A P / 2 , i.e. one half of the pressure difference across the interface at the lower end of the bridge. Computation is started by choosing values of X = R, the solid rod radius, where Z = 0 corresponds to the lower rod face, a value of Q, (= &) and a value of H. Then by trial and error meridian curves are found which also satisfy the upper boundary condition, i.e. which meet the upper solid of radius R at its edge for separation distance Z* between the rods.E. A. BOUCHER AND M. J. B. EVANS 92 1 The numerical solution of eqn (1) and (2) using the Runge-Kutta or Adams methodl29l3 can now be carried out on microcomputers. Details of the extensive computation are given in the Appendix. The volume of a fluid bridge is given by1 Va = n[R2Ze+ R (sin 42 - sin dl)] (3) in terms of the boundary meridian angles b1 and 42 (fig.1). The reduced force exerted on the upper solid owing to the presence of the bridge is feXt = nR sin 4h2 - nR2(H- 20). (4) Some bridges can be regarded as being partly pendent with respect to the upper solid and partly sessile with respect to the lower solid. Under such circumstances there will be a plane of zero force at some 2 where the azimuthal and meridional radii of curvature are equal. The pressure difference across the fluid/fluid interface where the bridge meets the upper solid is The and A P = 2(H-Ze). ( 5 ) area of the fluid/fluid interface is given by numerical integration of AaB = 2nsXdS the centre of mass 2. of the a phase constituting the fluid bridge is P Z * = n P Z d Z (7) s which can likewise be integrated numerically.The Helmholtz free energy is given in reduced terms by129 l4 F = Aap+2PZ.. (8) BRIDGE QUANTITIES For a given rod radius, R, the results of finding the meridian curves satisfying the boundary conditions can be expressed as the dependence of bridge volume Va on bridge height 20 for each value, d1, of the angle made by the meridian with the lower solid (fig. 1). This representation, illustrated in fig. 2 for R = 1, is a convenient intermediate step in the interpretation of the results. The region of existence of bridges on plots of Va against Z* can be interpreted in terms of experimentally measurable quantities. The limits of stability of bridges are sought since equilibrium configurations need not necessarily represent stable bridges.There are two ideal ways in which bridges can be manipu1ated:l (i) by varying the bridge height or rod separation Z* at constant bridge fluid volume Cm ; (ii) by varying the bridge volume Cm while keeping constant the rod separation 20. When making these changes, contact-angle conditions might arise leading to movement of the three-phase confluence from the solid edges. Boundary instabilities of this kind must be anticipated, and they are discussed below, but for the present their possible occurrence will be ignored. The results for R = 0.5, 1.0, 5.0 and 10.0 show most of the features which can be expected for fluid bridges. The sets of bridges between rods of reduced radius R = 1 will be seen to show a wide range of properties, and these 31 FAR 1922 FLUID BRIDGES 0 1 2 Z" Fig.2. Dependence of bridge volume Ya on height Z* for rods of radius R = 1 and several values of dl, the lower meridian angle at the solid contact. The locus of the maxima in Z* on the right-hand side of the loops defines the envelope, beyond which there are no bridges. The broken-line portions of the loops are for d2 c O", signifying that the meridians cannot meet the upper rod at its edge (see fig. 1). are discussed in some detail. The bridges where the rods are both of radius 0.5 or 10.0 are mainly discussed only where it is predicted that properties differ in kind from those for the other radii. To give an idea of actual dimensions, for air and water subjected to terrestrial gravity the capillary constant is a/mm = 3.8, so the rods would be of approximate radius 2-40 mm.The capillary constant for liquid silicon (surrounded by gas) at its freezing point (ca. 1690 K) is 8 mm, and R = 10 therefore corresponds to a rod of radius 80 mm. PREDICTED BRIDGE PROPERTIES RODS OF RADII R = 1 Bridges of constant volume P correspond to taking a horizontal line at some value of Va in fig. 2. No bridges exist to the right of what is referred to as the envelope. Also the range of q51 values (see fig. 1) cannot exceed 180" because the lower solid has a horizontal planar face, and, although one can find meridians with dl < 0", instability can be expected in this range: the vertical walls of the lower solid limit the meridian to q51 2 -90". Clearly one cannot have d2 < 0" because of the chosen geometry of the upper solid.The forcefext applied to the upper solid owing to the existence of a fluid bridge of fixed volume can be obtained from eqn (4). Fig. 3 showsfext varying with bridge height 20 for Va = 1, 2 and 4, corresponding to ca. 64, 128 and 256 mm3, were the bridge to consist of water surrounded by air. On elongating a bridge of reduced volume, Va = 1, it is seen that a sharp maximum is encountered in the applied force. By raising the upper rod with a device such as a piano-wire tensiometer, i.e. whereE. A. BOUCHER AND M. J. B. EVANS fex' -10 0 923 I - JJ 1 p 7 \ : 1 1 1 I 10 I 1 I 1 I 1 a twisted wire provides a torque raising the balance arm and vertical link to the solid, the applied force could be made to overshoot at this maximum.This has been referred to as an example of instability by stress-controlled manipulation, as opposed to strain-controlled manipulation where displacement is controlled and the force can pass through the maximum. The maximum is analogous to the condition for rupture of an adhesive joint made by placing liquid between two plates and pulling them apart by hand, and to the yield point in the stress-strain curve for a plastic solid in tension. In both of these analogues, gravity is not usually important. In the present case apart from buoyancy, the existence of gravity gives a body force. The bridge can exert a force on the rods owing to interface curvature (Laplace pressure) and the resolved component of the interfacial tension, i.e. that on the upper solid, is given by eqn (4).We defer discussion of the minima in the curves of fig. 3, but note the following additional features. By controlling the displacement of the upper solid it is predicted that the paths to the right of the maxima for all three curves can be followed reversibly. In the case of Va = 1 it is found that the curve ends abruptly when 41 = 171 k0.5" and 42 = O", but for Va = 2 and Va = 4 there are maxima in bridge height Z e which are expected to lead to bridge rupture. It seems that the under-cut portions after these maxima cannot be reached, although they represent bridges in hydrostatic equilibrium. We can illustrate how contact-angle conditions can intervene and lead to movement of the three-phase confluence away from an edge by examining again the case of = 1 near the force maximum.Using the concept of a so-called macroscopic contact angle, 8, which for purposes of argument is regarded as unique, i.e. there is no contact-angle hysteresis, fig. 4 shows how the ' hingeing' alp interface can reduce to 8 = 180-4, and the interface recedes. The principle can be summarized: if a fluid/fluid interface meeting a solid at its edge tends to make an angle with a neighbouring solid surface which is less than the intrinsic contact angle of that surface, then the interface will move from the edge rather than infringe the contact-angle condition. Problems of dynamic contact angles, microscopic contact angles and contact-angle hysteresis are well appreciated, and the previous statement refers to an ideal model. The treatment in this paper does not cover cases once boundary instability has been encountered.For the bridges under consideration the greatest range of behaviour without the intervention of instability due to contact angles will be possible for solid rod ends readily wetted by the a phase, but with sides not wetted 31-2924 0 fext 7 6 FLUID BRIDGES 0.3 0.4 Z" 10.7 0.8 0.9 X 0.4 0.2 z 0 1 Fig. 4. (A) Expanded portion of the maximum in fext from fig. 3 for Va = 1, with (a), (b) and (c) corresponding to the meridians in (B). (B) Three of the meridian curves near this maximum : (a) 41 = 135", (b) q51 = 150" and (c) = 155". The inset shows schematically how lower edge contact would be lost with 8 = 28" before meridian (c), for which d1 = 155", is reached. 1c 5 A P C -E 0 I I I I I 1 1 2 Z* Fig. 5.Dependence of the pressure difference A P e across the fluid/fluid interface where it meets the upper solid on bridge height for R = 1 and P = 1, 2 and 4 as indicated. A P contributes t o y x t shown in fig. 3 for these systems. by the a phase. The formal incorporation of contact-angle hysteresis into these studies of capillarity has already been de~cribed.~ Fig. 5 shows the dependence of the pressure difference A P e across the a/P interface, given by eqn ( 5 ) , at the upper three-phase confluence on bridge height. Comparison of fig. 3 and 5 shows that A P gives a major contribution, namely n R 2 A P , to the forcefext. The locations A P e = 0 correspond to a bridge whose interface has zero mean curvature at Z = Z*. The variation of the area, A@, of the fluid/fluid interface and the position, Z@, relative to the lower solid surface of the centre of mass of the constant-volume bridges are shown in fig.6(a) and (b), respectively. The bridges elongate with increase in Z@, but then they tend to squat giving a Z@ maximum before the height, Z*, maximum. When combined according to eqn (8) to give the Helmholtz free energy F as shown10 Aap 6 1 E. A. BOUCHER AND M. J. B. EVANS I I I 0.8 0.6 Z . 0.4 0.2 0 15 F 10 5 I I I I 1 Z" 2 I 1 I I , 2 t 1 z* b, 925 Fig. 6. (a), (b) and (c) show, respectively, the dependence of the fluid surface area A@, the position 2. of the centre of mass and the Helmholtz free energy Fon bridge height for R = 1 with volumes Va = 1, 2 and 4.926 FLUID BRIDGES I I I 1 I I 0 L 8 12 V" Fig.7. Dependence of the forcefext applied to the upper solid on bridge volume for rods of radius R = 1 and fixed bridge heights of Z e = 1, 1.5 and 2 (see text). 0 A P - 1 0 L 8 V" 12 Fig. 8. Variation of A P with Va for the systems of fig. 7. in fig. 6 ( c ) for R = 1, the area minimum dominates to give a minimum in F, then as 20 increases there is an inflexion in each curve with the abrupt termination of the curve for Va = 1. Cusps exist at the limit of bridge height for Va = 2 and 4 (the cusp does not show on the scale of the diagram for Va = 4). A full analysis of these thermodynamic quantities must await a more complete general treatment, extending previous accounts.12* 1 4 9 l5 By fixing the gap between the rods and varying the bridge volume one can follow the equilibrium path for the corresponding quantities.Starting on the d1 = -90" loop in fig. 2 the quantitiesfext, A P e , A@, Z* and F as a function of Za for bridge heights Z* = 1, 1.5 and 2 have been obtained: their termination corresponds to the curve through the ends of the loops in fig. 2. For R = 1 and Z e = 1 fig. 7 shows that a lower limit to Va can be expected. There is a maximum in the downwards pull on the upper solid, and then as the volume increases further the bridge bulges and pushes on the upper solid, as evidenced by the negativefext values. When Z e = 1.5 there is always a downwards pull on the upper solid. It is also noted thatE. A. BOUCHER AND M. J. B. EVANS 927 0 L VQ 8 12 I I I 1 I 1 0 4 8 12 VQ I I I 0 2 6 10 I VQ Fig.9. (a), (6) and (c) show, respectively, the dependence of Aab, 2. and F on bridge volume for the three bridge heights, Z e = 1, 1.5 and 2, for each of which the rod radius is R = 1 .928 2 1 0 f e x t -1 - 2 FLUID BRIDGES I a I I -4 Fig. 10. (A) and (B) show, respectively, the dependence of the applied forceFxt and the pressure difference A P on the bridge height Z e for rods of radius R = 0.5 and bridge volumes Va: of (a) 0.25, (b) 0.5, (c) 1 and (6) 2. there is a gap in the equilibrium path between vb! x 5.3 and 7.0 owing to the fact that the curve through the ends of the Va against Z* loops in fig. 2 crosses the vertical line 20 = 1.5 twice. Thefext against P curve for 20 = 2 shows a distinct cusp at the limit of the volume. The broken-line portion of the same curve represents the dependence offext on Va when 4h2 < 0".This region has no physical significance for rods having plane horizontal ends. A similar region exists for 2e = 1 . The contribution t o p x t from A P is shown in fig. 8 for R = 1, where the portions identified with respect t o p x t are evident, but there is not the similarity in shape and features which exist for the constant-volume cases. Fig. 9 (a) and (6) show, respectively, how A@ and 2. vary with P for the three bridge gaps. There is considerable similarity in the general shape of the curves for each quantity. Their combination to give the Helmholtz free energy in fig. 9(c) shows three similar curves which are themselves different from those in fig. 6(c) for the Helmholtz free energy at constant volume.With increase in bridge volume the angle 9, decreases and becomes zero (the fluid/fluid interface is therefore horizontal at the edge contact with the solid): at this location, within the slight uncertainty on the scale of the plots, there is an inflexion where d2F/d( P)2 = 0, and this probably denotes the limit for stable bridges (with the condition 8 > 90" for liquid on the sides of the lower rod). RODS OF RADII R = 0.5 A representative set of loops of Za as a function of 20 for various 4h1 values has already been published1 for the case of rods of radius R = 0.5. The main difference between them and the set in fig. 2 is that the properties of the system have now been predicted to conditions where the boundary becomes unstable. For cases of d1 d 90" the meridian curves have inflexion points.As the limit of boundary stability is reachedE. A. BOUCHER AND M. J. B. EVANS 929 I I 1 8 F 4 0 0.5 1 .o 1.5 2" Fig. 11. Variation of the Helmholtz free energy Fwith bridge height Z e for R = 0.5 and Va = (a) 0.25, (b) 0.5, ( c ) 1 and (6) 2. the S-shaped bridges develop a very thin neck. Beyond the stability limit meridians no longer have an inflexion, and eventually Q, along the meridian now increases to 180". For d1 > 90" the loops of Va against bridge height terminate at a value of H , approached froni below. Values of H in excess of the limit give meridians which share a position where Q, = 180", whereas the systems require bridges for which the meridians have a maximum in Q, at the inflexion, and these double back to meet the upper boundary condition (solid edge). The dependence offext on 20 in fig.10(A) is dominated by the existence of maxima infext and in Z*. The dependence of A P e on 20 in fig. 10(B) shows similarity with thefext plots for Va = 0.25, but for the Va = 2 curve it is now obvious that the resolved component of the interfacial tension is also contributing. The dependence of the Helmholtz free energy F on 20 in fig. 11 shows that the curves are very similar to those when R = 1. RODS OF RADII R = 5 With rods of reduced radii R = 5 one is severely limited by the range of possible meridian curves which do not have d2 < 0" where they meet the upper solid (see discussion regarding the implications of having a particular geometry of the solid). The Va against 20 curves can no longer be called loops since they terminate when d2 = 0", as is shown in fig.12. Fig. 13(a) and (b) showingfext and F, respectively, plotted against Z* for three fixed bridge volumes, Va = 25, 50 and 100, reveal very clearly the limited range of equilibrium paths which are possible under these circumstances;930 2 00 150 V" 100 50 0 FLUID BRIDGES 1 I 1 2 Z" Fig. 12. Dependence of bridge volume P on height Z* for R = 5 and the angles dl shown. The curves corresponding to the bridges terminate at the broken portion where d2 = 0". 200 100 f C X t 0 -100 -20 0 25 0 0 5 1 0 1 5 2" I I ( b ) 50 25 L I I 0 0 5 1 0 1 5 2" 200 I50 100 F 50 0 Fig. 13. (a) and (b) show, respectively, the dependence of applied forcefext and Helmholtz free energy F on bridge height for R = 5 for volumes P indicated.Comparison with fig. 3 and 5 indicates the restrictions on the ranges of the curves for R = 5 as opposed to those for R = 1.E. A. BOUCHER AND M. J. B. EVANS 93 1 800 Va 600 LOO 1.2 1 . 5 2 .o 2" 2.5 Fig. 14. Dependence of Va on 20 in the region 20 2 1.2 for rods of radius R = 10. The curves all terminate when $z = 0". The curves for < 0" are unlikely to represent physically realistic bridges, and so the possible range of their existence is severely limited. The broken line is for Va = nR220. RODS OF RADII R = 10 For rods of reduced radii R = 10 there are even greater restrictions on the existence of bridges owing to the limited ranges of Z* at constant P, and of Va at constant Z e .The majm factor in this limited behaviour is that tjb2 falls to zero. Fig. 14 indicates the possible systems. If the upper solid had a scooped-out end some extension of range would be possible for undercut equilibrium shapes, but these are not likely to be physical 1 y realis tic. SUMMARY OF RESULTS The range of solid radii examined illustrates the behaviour to be expected for fluid bridges, including equilibrium meridians giving stable bridges, bridges which might become unstable depending on the physical system being used to manipulate them, and bridges which might become unstable owing to boundary instability, where the system can lose edge contact. The plots of bridge volume at constant lower meridian angle, tjbl, against bridge height, 20, also show regions of Va and of Z e for which no bridges are at all possible.The following discussion deals with a model for zone refining where there are limitations on bridge volume and on the meridian angle. MODEL FOR FLOAT-ZONE REFINING Technological processes of purification of semiconductor materials, e.g. silicon, by using floating zones rely on capillarity. A vertical rod of material, say 5 cm diameter, is melted using a radio-frequency ring to give a so-called zone between rods which rotate. Purification is effected by moving the ring heater relative to the rods so that melting occurs from the remaining upper solid portion and solidification, leaving impurities in the liquid, takes place at the lower solid/liquid boundary. A fluid bridge of liquid between inert solids with horizontal ends is but a crude representation of the actual floating zone.The chief simplifications are those of ignoring the rotational motion, the temperature gradients, and therefore convection and Marangoni932 FLUID BRIDGES solid Fig. 15. Idealized representation of the freezing lower solid/liquid interface showing the liquid/vapour interface orientation for the floating-zone process. effects, and the shape of the solid/liquid interfaces.12 Two features of the real process can, however, be modelled: ( a ) the zone volume at an assumed uniform temperature will be proportional to the length L of solid which has been melted (in the simplest case Va = nR2L), but silicon shrinks by ca. 9% on melting; ( h ) it is found that the lower solid will only grow at constant radius equal to that of the upper solid if the growth angle, ly, shown in fig.15, is fixed, e.g. ca. 1 1 " for silicon. The bridges deal with Laplace or mechanical equilibrium and stability, and do not specifically account for physicochemical effects such as diffusion.12 In the floating-zone process physicochemical phenomena are important, since there is transfer of substance between phases. The pressure difference AP across the liquid/vapour interface and the pull of the interfacial tensions, especially at the three-phase boundary, affect the chemical potential of the liquid and the solid. The simplest assumption is that the solid/vapour, solid/liquid and liquid/vapour interfaces adjust to satisfy the Neumann triangle of forces, usually only applied to three-fluid phases.It will be presumed that under ideal circumstances this would be one condition for physicochemical equilibrium when the solid/liquid and solid/vapour interfaces can adjust their orientation, since they are being formed infinitesimally slowly. The growth angle, ry, can be adjusted, i.e. d1 is fixed, and the solid side remains vertical. The orientation of the solid/liquid interface away from the horizontal is for the moment ignored. Fig. 16 shows Va/nR2Ze plotted against the gap Z* between the solids for the various values of R from 0.5 to 10, each with d1 = 79", i.e. the growth angle is ry = 1 lo, appropriate to silicon.16 The amount by which the quantity Va/nR2Ze deviates from 1 is a measure of how the fluid volume differs from the cylindrical space between the rod ends.For ry = 11" there are no meridians, i.e. no equilibrium bridges, for P / n R 2 Z e 2 1.054 and R 2 0.5. Also, each curve terminates abruptly, although this is only shown for R = 5 and 10. The 9.1 % volume shrinkage of silicon on melting16 is indicated by the broken horizontal line. The largest possible reduced value of R is very close to 5, with the gap between the rods 29 = 2.02. The limit in rod diameter is therefore ca. 80 mm or just over 3 in., using 720 mN m-l for the surface tension of liquid silicon at its melting point1' and a liquid density of 2.53 g cm--3, i.e. the capillary constant a is 8 mm.16 A fuller account of this model applied to silicon and germanium, including experimental uncertainties and physicochemical as well as mechanical aspects, will be given separately.There are several assumptions (above) in this model which could mean that experimental values different than these predictions could be found in practice. If, for instance, the solid/liquid interfaces are both concave then more solid will have melted compared with the model (fig. 16) and the limiting rod radius would be smaller. Nevertheless, the present state of knowledge of the zone-refining process is such that the model gives a useful guide to stability limits, and by implication the rod--holmE. A. BOUCHER AND M. J. B. EVANS 933 I V* rR2Z" 0. .O 95 0.90 Fig. 16. Dependence of Va/xR2Ze on 20 for #J = 79" (w = 11") for each reduced radius R of (a) 0.5, (6) 1, ( c ) 2, (6) 5 and ( e ) 10. The horizontal broken line represents silicon, which shrinks on melting and has the growth angle of w = 11".systems provide an analogous model for Czochralski growth. Perhaps even more importantly these models suggest that these types of process might eventually be more fully understood in terms of a combination of mechanical and physicochemical effects. In crystal growth steady-state stability should be considered, although in the present case the slow growth rates are unlikely to give conclusions differing from those given by a quasi-static model. The steady-state stability of floating zones of silicon was examined by Surek and Coriell,lG who found that with a growth angle of ll", rods of radius 25.626 mm were stable to fluctuations in zone length such as would arise from unequal rates of advance at the freezing and at the melting interface.Whether this type of stability persists up to the limiting radius of ca. 40 mm found in this paper has yet to be investigated. What is also not clear is whether this steady-state stability exists only within a certain range of growth angles. Finally it is noted that the rotational motion of one rod relative to another, used in floating-zone processes to minimize thermal asymmetry, convective and Marangoni effect^,^ is sufficiently slow that it should not greatly affect zone stability under the usual gravitational influence, but in near-zero gravity conditions the rotational effect can be important.l8 APPENDIX COMPUTATION OF MERIDIANS* The present analysis of bridges has involved the computation of a large number of meridian curves by numerical integration of eqn (1) and (2). The following summary, based on a decade of study, is relevant to the vast majority of cases in this series of papers (and not just bridges) concerning the computation of meridian curves 2 = Z(x) and related quantities for axially symmetric capillary systems in a gravitational field.It is desirable to discuss the computation in one article, especially since a study of lenses also involving a very large number of meridians * Written in conjunction with T. G. J. Jones\o w P Table 1. Summary of advantages and disadvantages of computational methods method modified Euler Taylor series Runge-Kutta versions (a) 4th order (6) Gil (c) Butcher Merson Milne Hammings self- estimated error starting per step, h advantages Yes h3 rapid, easy to use Yes hn+l straightforward n derivatives error analysis no no h5 all self-starting, h5 high accuracy, widely h6 used h5 uses step-size adjuster for accuracy h5 two calculations per step, simple predictor eqn, iteration option h5 two calculations per step, iteration and simple predictor, no stability problems disadvantages low accuracy, slow complexity of higher convergence derivatives simple versions cannot test error at each step, complicated error analysis, large number of calculations per step cannot easily meet end values of independent variable (modify standard version) possibly subject to instabilities, not self-starting not self-starting fully recommended no no yes especially ( c ) not so suitable (a) and (W,E. A.BOUCHER AND M. J. B. EVANS 935 has been undertaken, and because microcomputers can be used with care for relatively large step sizes. The three main topics discussed are: (a) recommended methods of solving the appropriate Young-Laplace equation, using the results from several methods of integration which have been compared, (b) the estimation of errors involved in these numerical integrations and (c) the foundations underlying what are usually referred to as approximate solutions of the Young- Laplace equations. In the case of zero-gravity or neutral-buoyancy conditions, solutions of the meridian shapes often exist in terms of elliptic integrals. It is, however, still possible to use the same form of the equations and methods of integration for conditions of zero gravity, microgravity and terrestrial gravity.Choices of independent variable in alternative forms of eqn (1) have already been discussed:12* l3 the independent variable must always vary in the same direction. Many techniques of integration on several computers with various step sizes, independent variables, shapes of meridian and gravitational cases have been used to form the basis of this account. Table 1 summarizes our conclusions. In some cases the software accompanying a computer set-up will influence the choice of method, e.g. the Hewlett Packard HP9821A uses Gil’s method. When there is a free choice, Runge-Kutta’s, Merson’s or Gil’s, and related methods, are recommended: they are described in many textbooks.lg 23 Incorrect meridians resulting from using a slightly incorrect form of the Young-Laplace equation should be evident by the shapes of other than small arcs of meridians.It may not be easy to detect deviations from true meridians arising from non-matching initial steps. One test of consistency is to reverse the direction of computation over a considerable meridian length and expect the initial and final values to agree within the rounding error. Rounding errors will occur with a large number of steps involved in obtaining a meridian, but with a machine precision O(lO-g) they do not appear to be serious (see later). Instability during computation can occur with very shallow (nearly horizontal) meridians, where Bessel forms of the meridians are very successful. With sharply curved meridians, changes in the step size of the independent variable will usually signify the tendency, if any, of the computation to ‘cut corners’.Two recognized error tests are (a) to vary the step size as explained in detail by Duncan,24 although this method cannot be used with algorithms automatically using variable step sizes, and (6) to compute meridians for known analytical curves, e.g. circular arcs. Both of these have been used in the current study, and a third method, not previously discussed, has also been used. In general, the second-order Young-Laplace equation appropriate to a particular fluid body, i.e. bridge or drop, can be integrated once. The resulting equation can then be numerically integrated to give the meridian. That this does not seem to have been done before is itself curious.However, it also provides the basis for approximations for meridians and it can be used to estimate the error incurred during integration, although the integration is effectively reduced to that of a first-order rather than the original second-order equation. Examples will show the principles. Eqn (1) rewritten using 2 as the independent variable for fluid bridges integrates to give When Z = 0 at the bottom of the bridge, 0 = q51 and C = -cos @ giving Using Z as the independent variable, the left-hand side can be evaluated numerically and written as (IB-cosO) for comparison with the right-hand side for chosen dl, x“ (initial) and H, to give an ‘error’ of the form 6, = ( I B -COS @) -(2HZ-Z2 -COS 41) (1 1) as I , is characteristic of that meridian, with analogous integrals for pendent and sessile drops.Approximations for meridian shapes essentially give an approximation for these integrals,936 FLUID BRIDGES although they are not usually presented in this manner. Using a fourth-order Adams-Moulton method the results for two step sizes A Z = lop2 and were 6 = 4.65661 x and 6 = -3.53903 x respectively. These and a wider range of computations lead to the conclusion that rounding errors are evident but not important when A 2 = lop3, but are less evident for A 2 = which conveniently gives much faster computation time than for the smaller step size. For pendent drops the analogous treatment gives the same integral with Z as the independent variable, or with S as the independent variable, giving 6,, = (Ip,-cos@)-(2Hz-22- 1).(13) Again it is found that rounding errors of no great importance are incurred at small step sizes, and that good agreement is found with known shapes, e.g. circular arcs for A S = The form of eqn (1) and (2) presupposes that the system is subjected to a fixed, e.g. terrestrial, gravitational field. The use of a capillary variable, [, which can be varied differentially or given a set of values, has already been and 25 The form sin @ - = 2(H-c2)-T d@ dS gives a pendent drop for each H, subjected to terrestrial gravity for [ = 1, microgravity for [ 6 1 and a spherical segment in gravity-free conditions for [ = 0. Alternatively, a factor placed before ( H - Z ) can be used for the same purpose. Finally, note that Siekmann et ~ 1 .~ ~ advocate an expanded Runge-Kutta-Fehlberg method,27 of which they state: ‘The local error . . . O(h7). . . does not only allow a step width. . .larger than that of the classical Runge-Kutta method, but it also guarantees the same accuracy as the latter. . . . the accumulation of rounding errors is less and the saving of computer time is remarkable. ’ Feldberg’s version seems very similar to Merson’s 29 PRINCIPAL SYMBOLS a AaB fext F g h H IPD 1 L LIP r R S P x, 2 S capillary constant area of a/P interface reduced externally apply force on the upper solid disc (rod) Helmholtz free energy acceleration due to gravity (not necessarily the terrestrial value) integration step size in the Appendix shape factor: one half of the pressure difference across the alp interface at the lower solid contact actual linear dimension reduced linear dimension, l / a reduced pressure difference across the a/B interface, being A P at the lower solid where 2 = 0 and A P at the upper solid where Z = 2* actual rod radius reduced rod radius, r / a denotes solid phase meridian arc length reduced volume of a phase, ua/a3 meridian coordinates, being (A?, ZO) usually with ZO = 0 at the lower solid and (Xe, Z e ) at the three-phase confluence with the upper solid position of centre of mass of a phase above plane of lower solid (X.= 0) fluid phases with a the enclosed or bridge phase characteristic integrals in the AppendixE. A. BOUCHER AND M. J. B. EVANS 937 fluid a/fluid p interfacial tension computing error defined in the Appendix capillary variable formally allowing for the acceleration to vary from g due to terrestrial gravity contact angle for three-phase confluence with plane solid, not edge contact pa-# the positive density difference between the fluids meridian angle, arctan (dZ/dX) values of the meridian angle (fig.1) at the lower and upper three-phase confluences, respectively growth angle, (90 - E. A. Boucher and M. J. B. Evans, J. Colloid Interface Sci., 1980, 75, 409. E. A. Boucher, M. J. B. Evans and S. McGarry, J. Colloid Interface Sci., 1982, 89, 154. D. T. J. Hurle, Adv. Colloid Interface Sci., 1981, 15, 101. E. A. Boucher and H. J. Kent, Proc. R. SOC. London, Ser. A , 1977, 356, 61. E. A. Boucher and T. G. J. Jones, J. Chem. SOC., Faraday Trans. I , 1980, 76, 1419. E. Wolfram, Croat. Chem. Acta, 1972, 45, 137. S. R. Coriell, S. C. Hardy and M. R. Cordes, J. Colloid Interface Sci., 1977, 60, 126. P. H. Keck, M. Green and M. L. Polk, J. Appl. Phys., 1953, 24, 1479. W. Heywang, Z. Naturforsch., Teil A , 1956, 11, 238. lo M. A. Fortes, J. Colloid Interface Sci., 1982, 88, 338. l1 E. A. Boucher and T. G. J. Jones, J. Chem. SOC., Faraday Trans. I , 1982, 78, 1491. l2 E. A. Boucher, Rep. Prog. Phys., 1980, 43, 497. l3 E. A. Boucher, M. J. B. Evans and T. G. J. Jones, Adv. Colloid Interface Sci., to be published. l4 E. A. Boucher, Proc. R. SOC. London, Ser. A , 1978,358, 519. l5 T. Surek and B. Chalmers, J. Cryst. Growth, 1975, 29, 1. l6 T. Surek and S. R. Coriell, J . Cryst. Growth, 1977, 37, 253. l' P. H. Keck and W. van Horn, Phys. Rev., 1953,91, 512. l9 W. S. Dorn and D. D. McCracken, Numerical Methods with Case Studies (Wiley, New York, 1972). 2o Numerical Solutions of Ordinary and Partial Diflerential Equations, ed. L. Fox (Pergamon, Oxford, 21 M. L. James, G. M. Smith and J. C. Wolford, Applied Numerical Methodr for Digital Computation 22 G. M. Phillips and P. J. Taylor, Theory and Application of Numerical Analysis (Academic Press, 23 H. H. Goldstine, A History of the Calculus of Variations from the 17th Through the 19th Century 24 W. J. Duncan, Philos. Mag., 1948, 39, 493. 25 E. A. Boucher and T. G. J. Jones, J. Colloid Interface Sci., 1983, 91, 301. 26 J. Siekmann, W. Scheideler and P. Tietze, Comput. Meth. Appl. Mech. Eng., 1981, 28, 103. 27 G. Jordan-Engelm and F. Reutter, B. I. Hochschultashenbucher (Biblograph. Inst., Mannheim, 1976), 28 L. F. Shampine and R. C. Allen, Numerical Computing: An Introduction (W. B. Saunders, Philadelp- 29 L. Fox and D. F. Mayers, Computing Methodr for Scientists and Engineers (Oxford University Press, W. Wuest, ESA Spec. Publ. No. 114 (European Space Agency, Pans, 1976). 1962). with Fortran and CSMP (IEP, New York, 2nd edn, 1977). London, 1973). (Springer-Verlag, New York, 1980). Bd 106, 2 Aufl. hia, 1973). Oxford, 1968). (PAPER 4/839)