Deformation of a viscoelastic droplet passing through a microfluidic contraction

J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

Contents lists available at ScienceDirect

Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm

Deformation of a viscoelastic droplet passing through a microfluidic contraction D.J.E. Harvie a,∗ , J.J. Cooper-White b , M.R. Davidson a a b

Department of Chemical and Biomolecular Engineering, University of Melbourne, Victoria 3010, Australia Division of Chemical Engineering, University of Queensland, Queensland 4072, Australia

a r t i c l e

i n f o

Article history: Received 23 October 2007 Received in revised form 6 May 2008 Accepted 6 May 2008 Keywords: Viscoelastic Volume of fluid Contraction Microfluidics Non-Newtonian Droplet Oldroyd-B Drug delivery

a b s t r a c t A computational fluid dynamics algorithm that uses a volume of fluid method to model immiscible phase behaviour is adapted to simulate the flow of viscoelastic fluids. Two-dimensional Cartesian simulations of a viscoelastic droplet passing through a contraction are compared against experimental results generated using a three-dimensional planar contraction. Viscoelastic effects are simulated using the Oldroyd-B rheological model. The physical and numerical models capture the deformation process well, particularly the characteristic forked tail that develops on the viscoelastic droplet while it is within the contraction. Encapsulation of continuous phase fluid is observed in the experiments, but not predicted by the simulations. We suggest this is due to inadequate mesh refinement at the length scale of a fine thread entrainment process. The numerical results show the existence of an elastic stress within the droplet that could cause the entrainment and subsequent encapsulation. The process has application potential in the manufacture of micron-sized drug and chemical delivery capsules. © 2008 Elsevier B.V. All rights reserved.

1. Introduction With the increasing popularity of multiphase microfluidic devices for analysis, diagnostic and research purposes comes the need to accurately simulate fluid flow in micron-sized systems [41,31,33]. To this end we have been using volume of fluid (VOF) simulations to examine how Newtonian and generalised non-Newtonian droplets deform when moving through a simple microfluidic geometry — a cylindrical 4 : 1 contraction [13–15]. In this study we consider a viscoelastic droplet passing through a planar contraction. Many fluids, which are either already in use, or offer strong application potential in microfluidic devices, are viscoelastic. Examples include various biofluids (such as DNA and protein suspensions), polymer melts and polymer suspensions [3,40,32,11]. On the macrolength scale, the flow of a single phase viscoelastic fluid through planar, square and cylindrical contractions has emerged as a benchmark test for both experimental and numerical methods [42,5,4,2]. The contraction owes its popularity not only to the ease with which it can be manufactured and represented numerically, but also to the regions of strong extensional and shear strain which exist within it that critically test the elastic properties of the fluid. One difficulty with simulating the single phase

∗ Corresponding author. E-mail address: [email protected] (D.J.E. Harvie). 0377-0257/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2008.05.002

system is in dealing with the stress singularity that occurs at the ‘re-entrant’ corners at the entrance to the contraction, particularly for highly elastic fluids. Some researchers have avoided this difficulty by considering geometries that have corners rounded with a small radius [6,1]; simulation of highly elastic fluids in the sharpcornered geometry remains problematic however. On the micron scale, the single phase problem has not received much attention until recently: Rodd et al. [27] reviews the few relevant works. Using microfluidic contractions experimentalists have been able to probe viscoelastic behaviour at elasticity numbers not previously accessible using macro-scaled systems. The multiphase viscoelastic contraction problem has received less attention but holds great promise both as a benchmark test for researchers and because of its application potential. A viscoelastic droplet which is contained within a Newtonian fluid and flows through a contraction experiences regions of strong extensional and shear strain (as in the single phase problem), but does not directly experience the re-entrant corner stress singularity as it is separated from the device walls by Newtonian fluid (assuming that the droplet does not wet the walls). Thus, this problem has the elasticity-testing properties of the single phase problem but is easier to tackle numerically. An added advantage over the single phase problem is that the droplet shape, which is directly influenced by the elastic stress field present within the droplet, can be compared between simulations and experiments. In an application sense multiphase problems will become more relevant as devices

68

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

which utilise interfacial forces to accomplish fluid containment, reaction and mixing become more prevalent [31,33]. While numerical methods capable of simulating single phase viscoelastic systems are reaching maturity, methods capable of simulating immiscible multiphase viscoelastic systems have received less attention [35]. In a pioneering study, Keunings [19] used a finite element method on a deforming Lagrangian mesh to study the breakup of a jet of Oldroyd-B fluid. Ramaswamy and Leal [26] used a finite difference method, also on a deforming mesh, to study the behaviour of a FENE-CR droplet in uniaxial extension flows. Pillapakkam and Singh [23] appear to be the first researchers to perform immiscible phase viscoelastic simulations on a stationary Eulerian mesh. They used a finite element method and tracked the location of the boundary between phases using the Level-Set method. More recently Davidson et al. [10] adapted a VOF code to model Oldroyd-B fluids and were able to reproduce the ‘beads-ona-string’ structures observed during dilute PEO solution pendant drop experiments. Khismatullin et al. [20] used a VOF method to simulate droplets of Giesekus fluid deforming in simple shear flows, while Tome et al. [38,37] used a marker and cell (MAC) approach to simulate various channel and cylindrical free surface flows involving viscoelastic fluids. Eulerian mesh-based immiscible boundary techniques such as level-set, VOF and MAC offer great promise in simulating multiphase viscoelastic flows; however, further validation of the methods, in terms of comparing experimental and simulation results, is required. In this paper we detail a viscoelastic computational fluid dynamics (CFD) algorithm that uses a VOF method to model immiscible phase boundaries. The algorithm is used to simulate the deformation of a two-dimensional Cartesian Oldroyd-B droplet as it passes through a contraction, and the results compared against results obtained experimentally for a ‘flattened’ Boger fluid droplet passing through a three-dimensional planar contraction. The comparison serves two purposes: (a) to validate the physical and numerical models used in the simulations, and (b) to examine how viscoelastic fluids behave in this system. On the second point we are particularly interested in how the polymers move, orientate and stretch as the droplet deforms, and how the resulting polymer configuration influences this deformation.

2. Simulation method

Fig. 1. Flow geometry considered numerically. All lengths are nondimensionalised ∗ . by the inlet region half width, Win

imately replicate the rounded corners of the experimental flow device. The Oldroyd-B rheological model is used to capture viscoelastic behaviour. Under this model individual polymers are represented as ‘dumbbells’ — that is, two spheres connected via a thin spring. The orientation and length of individual dumbbells changes in response to local fluid strain rates. The Oldroyd-B model is one of the simplest dumbbell models in that the spring that connects the spheres is linear and infinitely extensible. The model is reasonably well suited to representing Boger fluids as both Boger and OldroydB fluids have a shear viscosity that is independent of the shear rate. Boger fluids typically consist of Newtonian solvents containing a dilute concentration of high molecular weight polymers. In our numerical implementation of the Oldroyd-B equations the polymer concentration is assumed to be uniform within each phase, but different between phases. We assume that the polymer relaxation time is the same in both phases. In the simulations reported here only the disperse phase contains polymers. The five equations describing the motion of the two phases, including viscoelastic effects, are

∇ ·u=0

(1)

∂ + ∇ · u = 0 ∂t

(2)

∂u 1 1 ˆ + ∇ · uu = −∇ p + ∇ · nı(|x − xs |) + We Re ∂t

(3)

∂A 1 + ∇ · Au = A · ∇ u + (∇ u)T · A − (A − I) De ∂t

(4)

 = [∇ u + (∇ u)T ] +

2.1. Physical assumptions and relevant equations The system under consideration consists of a droplet of viscoelastic liquid that is surrounded by a Newtonian liquid and is passing through a planar contraction. The terms ‘disperse phase’ and ‘continuous phase’ are used to refer to the droplet and surrounding liquids, respectively. Surface tension acts over the interface between the phases, and the walls of the flow device are assumed to be nonwetting with respect to the disperse phase, and also nonslip. Fig. 1 shows the geometry of the numerical flow domain. It is comprised of three regions: an inlet, a contraction and an outlet. Two-dimensional Cartesian coordinates are used as an approximation to the three-dimensional geometry used in the experiments. ∗ (we All lengths are nondimensionalised by the inlet half-width Win use an asterisk to indicate a dimensional quantity). To approximate the experimental conditions described in Section 3 the nondimensional contraction half-width is set to Wcon = 0.209, the initial droplet diameter to D = 0.864 and the length of the contraction region to Lcon = 5.364. The inlet region length is Lin = 3 and the total length of the flow domain is L = 12. The contraction corners are ‘bevelled’ with a nondimensional side length of 0.1 to approx-

p (A − I) De

(5)

These equations are presented and employed in a nondimensional form with the local velocity u scaled using the average fluid veloc∗ /u ¯ ∗ , density  ity entering the inlet channel u¯ ∗ , time t using Win using the continuous phase density c∗ , solvent viscosity  using the continuous phase solvent viscosity ∗c , and polymer concentration p also using the continuous phase solvent viscosity ∗c . In this formulation the polymer concentration is defined as the increase in solution viscosity due to the presence of polymers when the solution is subjected to steady shear [8]. Hence, we use the terms polymer viscosity and polymer concentration interchangeably. With these scalings, the nondimensional groups appearing in the above equations become Re =

∗ c∗ u¯ ∗ Win

∗c

,

We =

∗ c∗ u¯ ∗2 Win

∗

,

and

De =

tp ∗ u¯ ∗ ∗ , Win

where is  ∗ the surface tension (or energy per unit area) of the interface between the phases and tp∗ is the polymer relaxation time. A surface tension number S can also be defined as S=

1 We(1 + 1/Re)

(6)

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

This number measures the strength of surface tension forces relative to inertial and viscous forces combined [14]. The local disperse phase volume fraction is . As local fluid properties depend on the fluid type at a particular location, they can be defined as functions of . Specifically, the local solvent viscosity is () = d  + 1 −  where d = ∗d /∗c is the nondimensional disperse phase solvent viscosity. Similarly, the local fluid density is () = d  + 1 −  where d = d∗ /c∗ is the nondimensional disperse phase density. Finally, as the polymer concentration within the continuous phase is zero in the considered problem (∗p,c = 0), the local polymer concentration is p () = p,d  where p,c = ∗p,c /∗c is the nondimensional continuous phase polymer viscosity measured during steady shear. These functions are utilised when describing the numerical algorithm in Section 2.2. Eqs. (1)–(3) are a standard starting point for incompressible multiphase flow analysis, with the possible exception of the second term on the right of momentum Eq. (3): this term represents the surface tension force that acts across the interface between immiscible phases. In this term  is the signed curvature of the ˆ a unit vector normal to the interface and directed into interface, n the disperse phase, ı a Dirac delta function and xs the interface location. Eq. (4) describes the evolution of the dumbbell configuration tensor A [8]. This tensor is defined as A = RR where R is the endto-end vector of individual dumbbells, and the ensemble average is taken locally over many dumbbells. The length of R is normalised so that in the relaxed state, |R| = 1. Initially A = I within the droplet (the identity matrix), representing polymers that are both relaxed and randomly orientated. Two properties of A are  used when analysing the results. The tr(A) is used as a measure of average square root of the trace  √ polymer length. When the polymers are relaxed, tr(A) = 2 in two dimensions. The average orientation of the polymers is given by the principal eigenvector of A [26]. When A = I this eigenvector is indeterminant, corresponding to no preferential polymer orientation. The total stress within the fluid is given by Eq. (5). This stress is composed of a Newtonian and viscoelastic component, the latter being a linear function of A.

69

on a conservative discretisation of Eq. (2) for cell i: n |i = o |i −

1

V |i



s F |i→j(i,sk)

where F |j is the volume of disperse phase fluid passing over face j during the timestep, and the summation cycles through all of the four faces surrounding cell i. The disperse phase fluxes F |j are calculated using geometrically based interface reconstruction techniques that preserve the discrete nature of |i [28]. Further implementation details are provided in Appendix A. VOF methods have been used extensively over the past two decades to simulate a wide variety of experimentally observed immiscible fluid flows [22,21,17,29]. Modern VOF methods have an advantage over alternative immiscible phase methods such as level-set [34,30] and front-tracking [24,39] in that fluid conservation is assured. A disadvantage of VOF methods over these alternatives however is that the calculated local surface tension forces are generally less accurate, leading to larger errors in the resulting velocity field [7]. For our VOF implementation and the employed mesh size of 128 × 1536 the maximum magnitude of these velocity errors is less than 1% of the average fluid velocity occurring within the contraction [16]. Hence, these errors are of minor significance in this study. 2.2.2. Viscoelastic configuration tensor solver Just as the discretisation of Eq. (2) recognises that the VOF function  varies abruptly within the computational domain, discretisation of the configuration tensor Eq. (4) must also recognise that abrupt changes in the polymer concentration p occur within the domain. Our starting point for the method is





∂p A 1 (A − I) + ∇ · up A = p A · ∇ u + (∇ u)T · A − De ∂t

(8)

which results from combining Eqs. (2) and (4), noting that p is just a linear function of . As for Eqs. (1)–(3), Eq. (8) is averaged and discretised using a finite volume approach as np An −op Ao

t



T

+∇ · ua op Ao =np An · ∇ ua +(∇ ua ) · An −

2.2. Numerical algorithms 2.2.1. Multiphase volume of fluid solver Eqs. (1)–(3) are averaged over time and space and solved using a finite volume method. A staggered, uniform, structured mesh is used, with pressures located at cell centres and velocity components located at cell faces. Details of the algorithm as applied to Newtonian liquids are given in Rudman [28], as well as in the previous cylindrical contraction studies by the authors [13,14]. Preliminary versions of the method as applied to viscoelastic fluids are described in the pendant drop study of Davidson et al. [10] and the cylindrical contraction study of Harvie et al. [12]. A volume of fluid (VOF) method is used to track the location of the phase interface. Under the VOF method a ‘VOF function’ is defined as the disperse phase volume fraction averaged over each computational cell. The magnitude of this function defines the location of the phase interface: If |i is the VOF function for computational cell i, then the cell is assumed to contain only continuous phase fluid if |i ≤ , only disperse phase fluid if |i ≥ 1 − , or both fluids and a phase interface if < |i < 1 − . In the reported simulations we use = 10−7 . In this study an adaption of Youngs’ directionally split algorithm is used to advect the VOF function [43,28]. This algorithm is based

(7)

k=1,2 s=±1



1 (An − I) De (9)

The implicit treatment ensures that remains stable in the limit of large ∇ ua components or small De. Eq. (9) is implemented using two sequential updates; an advection update: B∗ − op Ao

t

+ ∇ · ua op Ao = o

(10)

followed by a source term update: 1 An − A∗ T = An · ∇ ua + (∇ ua ) · An − (An − I) De

t

(11)

Here B∗ and A∗ are related by A∗ =

B∗ max(np , p )

(12)

where p = min(p ( ), p (1 − )) is the minimum non-zero polymer concentration that can be represented by the viscoelastic VOF method. Noting that the configuration tensors A|i are located at cell centres, the tensor advection update can be performed using the same

70

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

dimensionally split structure as used for the VOF advection update. Discretising Eq. (10) over cell i in conservative form gives B∗i = op |i Ao |i −

1

V |i

 

sFB |i→j(i,sk)

(13)

k=1,2 s=±1

where FB |j is the amount of polymer and associated configuration tensor passing over face j during the timestep. These quantities are calculated numerically using FB |j = [Ao ]j Fp |j and Fp |j = p (F |j )

(14)

(1) If the local value of A is specified by the linear combination (1 − x)L + xR, where x is a distance variable satisfying 0 ≤ x ≤ 1 and both L and R are symmetric two-dimensional tensors that have positive diagonal components and a positive determinant, then any average of A taken within the range 0 ≤ x ≤ 1 also possesses these two properties. (2) If the two-dimensional symmetric tensors A1 and A2 both have positive diagonal components and a positive determinant, then the linear combination of these tensors ˛1 A1 + ˛1 A1 also has positive diagonal components and a positive determinant, assuming that ˛1 , ˛2 ≥ 0.

o

where [A ]j is the average A that occurs along face j during the timestep, Fp |j is the amount of polymer moving through the face during the timestep, F |j is the amount of disperse phase fluid moving through the face during the timestep, and p ( ) is the linear function relating polymer concentration and disperse phase volume fraction as defined in the previous section. By using the VOF-calculated F |j fluxes to evaluate FB |j in Eq. (14), we ensure that viscoelastic stresses cannot be transported across immiscible phase boundaries (as is required physically) as there is no disperse phase transport across such boundaries. A one-dimensional ‘characteristic type’ method is used to calculate the timestep-averaged and face-centred configuration tensors, [Ao ]j . The basis of this method is illustrated schematically in Fig. 2: a physically consistent linear profile for Ao is first constructed in each cell using values from surrounding cells. The gradient of this profile is given by the product aG. The tensor G is a notional gradient of Ao calculated using a second-order central difference. The scalar a limits this gradient so that certain physical properties of the reconstructed A profile are satisfied at both the upper and lower face locations of the cell (in the current direction): namely that the diagonal components are positive (required physically as these are defined as the average of squared R components) and that the determinant of A is positive. That |A| must be positive is shown by Joseph [18]. Note that B, being a scalar multiple of A, must also share these properties. The face-centred tensors FB are then calculated by integrating backwards in time from face j over these reconstructed profiles, finding the average tensor configuration associated with the polymer amount that passes through face j during the timestep. The variable r specifies what proportion of the polymer contained within the upstream cell ic moves through face j during the timestep. Two results concerning symmetric two-dimensional tensors are now given to aid our discussion. They can be shown via algebraic expansion:

Fig. 2. A schematic showing how [Ao ]j is calculated for each face j. The scalar limiter a is chosen so that certain physical properties of the reconstructed Ao profile are satisfied at the two end points indicated by red triangles (grey in print). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

As the tensor gradient limiter a is chosen to ensure that the reconstructed A profile has positive diagonal components and a positive determinant at both the upper and lower faces of a cell, then via the first of the above tensor results we can say that any average of the reconstructed A taken within the cell volume will also have positive diagonal components and a positive determinant. Two implications follow from this: (1) When a proportion of B is removed from a cell via a FB flux through a cell face, the average B left behind will have positive diagonal components and a positive determinant. (2) As the quantity of B that leaves a cell through a face during a timestep is calculated as the average of B taken over a volume within that cell, each FB tensor will have positive diagonal components and a positive determinant. Further, when polymer enters a cell through a face during a timestep, the updated B value for that cell is a linear combination of the previous B value and a FB value, both tensors that have positive diagonal components and positive determinants. Using the second of the above tensor results, it follows that the updated B value will also have positive diagonal components and a positive determinant. Thus, regardless of whether a cell gains or looses polymer during an advection update, the diagonal components and determinant of its B tensor will remain positive. As a result, no additional algorithms are required to keep the configuration tensors physically meaningful under this advection scheme, in contrast to some other methods [35,20]. This conclusion is supported by simulation results. Implementation details of the viscoelastic tensor advection algorithm are given in Appendix A. The source term update for the configuration tensor is performed by solving Eq. (11) in each cell using a Newton–Ralphson technique [25]. A conservative analysis of the discretisation shows that the diagonal components of An will remain positive during this step provided that the maximum Courant number within the √ domain C = max(u t/ x) < 1/ 32. We have not been able to determine a similar criterion that guarantees that |An | remains positive after this step. Instead, in the code a check is made for |An | < 0 after the source update has been performed, and if any are found the step is repeated with a reduced t. This repetition was not required in any of the simulations presented here; however in other problems the technique has always allowed the computation to proceed. 2.2.3. Stress discretisation, boundary conditions and timestep limitations Finally, the stress acting at each cell face is calculated by discretising Eq. (5) using average-timestep values for the viscoelastic tensor Ba = ap Aa : T

 = [∇ ua + (∇ ua ) ] +

1 (Ba − ap I) De

(15)

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

Only the off diagonal component of Ba needs interpolation from its stored location in the application of Eq. (15)— a simple average of values from surrounding cells is used. The alternative calculation of this component by interpolating Aa and ap separately and then taking their product led to flow instabilities developing around the interface region. Boundary conditions on  are implemented by setting values of Ba within ‘dummy’ cells located within the domain walls, and then applying Eq. (15) over the wall region. As the disperse phase does not wet the contraction walls in this problem,  = 0 in the dummy cells. Further, as volume fraction and polymer concentration are linearly related, p ( = 0) = 0 and Ba = 0 in the dummy cells. This boundary condition is equivalent to assuming that the elastic part of the total fluid stress is zero along a domain wall. This is physically justifiable on the grounds that a thin film of Newtonian continuous phase fluid must exist between any wall and disperse phase fluid, as the disperse phase is assumed to be completely nonwetting with respect to the walls. Along with the inertial, viscous and interfacial tension timestep limitations present in the Newtonian VOF code [28], the viscoelastic VOF code has an additional elastic timestep limitation. In brief, this limitation can be calculated in an order of magnitude sense by examining how an error in a beginning-of-timestep velocity propagates via configuration Eq. (9) to become an error in an endof-timestep configuration tensor, and then subsequently becomes an error in an end-of-timestep velocity via Eqs. (15) and (3). The analysis shows that for a velocity error to decrease over time the timestep must satisfy



t < c xmin

Re De min(d , 1) max[tr(Ba )max , 10−10 ]

(16)

where tr(Ba )max is the maximum trace of Ba that occurs within the computational domain, xmin is the minimum cell dimension and O(c) = 1. In the code we set c = 0.1 in Eq. (16) and check that this constraint is satisfied at the end of each timestep. If not, the timestep is reduced and the step repeated. 3. Experimental method In this study we compare various simulation results against a single droplet deformation experiment performed in a threedimensional planar contraction. For the experiment a solution of water, glycerol and 106 g/mol polyethylene oxide (PEO) polymer was made for the disperse phase using the techniques described in Tirtaatmadja et al. [36](d∗ = 1093 kg/m3 , ∗d = 3.4 × 10−3 Pa s, ∗p,d = 2.4 × 10−3 Pa s). The continuous phase was 20 cS silicon oil

(c∗ = 935 kg/m3 , ∗c = 1.87 × 10−2 Pa s). The flow device was constructed from PDMS using standard lithography techniques and had ∗ = 101 ␮ m and a uniform depth of approxan inlet half-width of Win ∗ imately h = 34 ␮m. Other device dimensions are as given for the numerical flow geometry in Section 2.1. The average fluid velocity flowing into the inlet channel was calculated to be u¯ ∗ = 0.163 m/s. A polymer relaxation time of tp∗ = 3.7 × 10−4 s was assumed based upon the dilute concentration Zimm theory analysis presented in Tirtaatmadja et al. [36]. We have chosen a surface energy (interfacial tension) of  ∗ = 40.5 mN/m for use in the simulations. This is the value measured between the silicon oil and water/glycerol solution in the absence of PEO. As shown in Cooper-White et al. [9], PEO can be surface active. However, the data also shows that for a water/glycerol/PEO and air system, when the PEO concentration is low the energy of a freshly created surface is approximately constant until the age of the surface becomes greater than around 10 ms. In a liquid/liquid system this response time is expected to be larger. In the current experi-

71

ment the droplet is within the contraction for less than 1 ms. As this time is significantly less than the surface energy response time, the characteristic surface energy for the droplet while deforming in the contraction will be closer to the value measured in the absence of polymers, rather than the equilibrium value measured with polymers present. Hence our choice for  ∗ . It would be desirable to model the absorption, desorption and diffusion of polymers to and from the phase interface during the deformation process, however this would require considerable numerical development. 4. Results and discussion 4.1. Overall deformation behaviour Selected images from the experiment are shown in Fig. 3 a. Results from an equivalent viscoelastic simulation, as well as results from a Newtonian simulation are shown Fig. 3 b and c, respectively. The nondimensional numbers used in the simulations are listed in the figure. The Newtonian simulation is identical to the viscoelastic simulation except that it was conducted with a Newtonian disperse phase having a nondimensional viscosity equal to the total steady-state shear viscosity of the viscoelastic fluid used in the experiment (d = 0.182 + 0.128 = 0.310). Both simulation results were produced using a 128 × 1536 computational mesh (including the contraction obstacles) of uniform square cells: as shown in the mesh sensitivity study (Appendix B) with this level of mesh refinement all but the smallest scale features visible in these images are quantitatively resolved. The overall deformation behaviour of each of the droplets shown in Fig. 3 is similar, and largely resembles that observed when a moderate viscosity Newtonian droplet passes through a cylindrical contraction at low Reynolds numbers [13]. As the droplet enters the contraction, extensional strain within the fluid deforms the droplet into a long filament. Within the contraction surface tension forces act to increase the radius of the filament. As these forces are moderately large here (S = 7.26), in all cases shown the droplet occupies a large proportion of the channel width while it is within the contraction. At the exit of the contraction, fluid decelerates as the channel width increases. This deceleration results in a negative extensional strain that reduces the droplet length as it exits the contraction, forming the ‘part-moon’ shape observed at t ≈ 1.32. Once away from the contraction exit and within the slower moving outlet region, the significant surface tension forces return the droplet to the nearly circular shape observed at times beyond t = 2. A closer examination of the three droplets of Fig. 3 shows that there are some small but significant differences in the way that each behaves. There are some small differences between how both of the simulated droplets behave and how the real droplet behaves. There are also differences between how the Newtonian droplet behaves and how the viscoelastic droplets behave. These differences are potentially more interesting. We discuss both in the following. 4.2. Comparison: simulated and experimental droplets The early and late timeframes of Fig. 3 show that both of the simulated droplets move through the inlet and outlet regions slightly faster than the experimental droplet. This discrepancy is most likely a symptom of the two-dimensional assumption employed in the simulations. By neglecting the third spatial dimension two physical effects are missing from the simulations; (1) surface curvature in planes normal to the plane the contraction lies in, and

72

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

Fig. 3. Viscoelastic planar experiment and two-dimensional viscoelastic and Newtonian simulations: Re = 8.25 × 10−1 , S = 7.26 (We = 6.22 × 10−2 ) and d = 1.17. In the viscoelastic results of frame (b) the shading represents the average polymer length



tr(A) and uses the scale of Fig. 4 a. Approximate nondimensional times measured

relative to the leading tip of droplet entering the contraction are given at the bottom of each image. (a) Experimental results. (b) Viscoelastic Oldroyd-B simulation: d = 0.182, p,d = 0.128 and De = 5.98 × 10−1 . (C) Newtonian simulation: d = 0.310.

(2) viscous drag resulting from the nonuniform velocity profile in planes normal to the contraction plane. In the experiment the combination of these two factors results in a retarding force on the droplet due by the presence of the upper and lower device walls. This force slows the real droplet’s progress, particularly in the inlet and outlet regions of the flow device where the axial pressure gradient within the continuous phase is quite low. Within the contraction itself the axial pressure gradient driving the droplet motion is much higher however, and the ability of this retarding force to slow the droplet’s motion is less. As a result, a difference in axial speed between the real

and simulated droplets is only observed in the wider inlet and outlet regions of the flow device, and not within the narrower contraction. Another difference between the simulation and experimental results is that the simulated droplets occupy less of the channel width while within the contraction than the real droplet does. This difference can also be attributed to the two-dimensional assumption used in the simulations. This assumption means that interface curvature in planes normal to the contraction plane is neglected. Consequently, the surface tension force, which acts to increase the radius of the simulated droplet while it is within the contraction, is underestimated.

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

Note that there are other factors that may also cause discrepancies between the experimental and simulation results: these include wall roughness, finite wettability of the disperse phase, flow device elasticity and errors in the measurement of flowrates, device dimensions and fluid properties. Also, as discussed above, at the employed mesh resolution some of the smaller scale features visible in the simulation images will not be quantitatively resolved. 4.3. Comparison: experimental and simulated viscoelastic droplets While various modelling assumptions employed in the simulations account for the differences in behaviour between the experimental and simulated droplets, the differences in behaviour between the viscoelastic droplets and the Newtonian droplet are caused by polymer behaviour. To illustrate the physical processes that occur in the viscoelastic cases, and to allow a more detailed comparison between the simulated and experimental results, enlarged images of the viscoelastic droplets of Fig. 3 are presented in Fig. 4. Unlike in Fig. 3, image times were selected so that comparable droplets are located at approximately the same positions within the flow device. Average polymer length, average polymer orientation and the recirculation patterns occurring within the simulated droplet are also shown. Fig. 4 a and c show that when the tail of the droplet has just entered the contraction the simulation and experimental forms of the viscoelastic droplet are qualitatively similar; however, the width of the real droplet is greater than that of the simulated one, and the tail of the real droplet appears ‘sharper’ than that of the simulated droplet. As discussed above, the simulated droplets are generally thinner than the real droplet while within the contraction as the twodimensional assumption used in the simulations underestimates the surface tension force acting on the droplet. Why a sharper tail develops on the real droplet than the simulated droplet is most likely due to geometrical factors: As shown in Fig. 3 a, the real droplet is positioned slightly to the right of the contraction centreline prior to entering the contraction. Also, the experimental contraction is slightly asymmetrical about the device centreline, particularly at its entrance. This asymmetry will influence the dynamics of the real droplet as it enters the contraction, potentially accounting for the difference in tail shapes observed between simulation and experiment. Average polymer configurations within the simulated droplet are displayed in Fig. 4 a. The polymers are generally orientated in an axial direction at this time, with polymers at the rear of the droplet being longer than polymers at the front. There is a thin ‘shell’ of extended polymers existing around the front and sides of the droplet. Eq. (4) shows that as well as being advected (transported) by the local velocity, polymers change their length and orientation in response to local fluid strain rates. They also experience relaxation, which continuously drives their length towards the equilibrium value without changing their orientation. With this behaviour in mind, the polymer configuration pattern of Fig. 4 a can be explained. As polymers enter the contraction, they experience a strong extensional strain which extends them in the axial direction. Once within the contraction, the polymers start relaxing; however, at the same time internal recirculation patterns that have formed within the constrained droplet start stretching the polymers in new directions. As shown in Fig. 4 b, at t ≈ 0.42 these recirculation patterns consist of two counter-rotating vortices which move droplet fluid forwards along the channel centreline towards the front of the droplet, and move droplet fluid backwards around the phase interface towards the back of the droplet

73

again. These patterns alter polymer configuration in three principal ways: (1) Next to the phase interface, but within the droplet, the tangential stress caused by the movement of continuous phase fluid past the droplet results in a thin region of high shear strain developing. Within this region polymers are extended in a direction parallel to the phase interface. The results of this are particularly evident around the leading tip and sides of the droplet where extended polymers appear as a thin ‘shell’ lining the internal surface of the droplet. (2) Along the droplet centreline and towards the front of the droplet a region of negative extensional strain exists where the droplet fluid is decelerating. This strain reduces the polymer lengths here, forming the low polymer extension region observed in the figure. Interestingly, as the polymers are generally orientated in the axial direction prior to experiencing the negative axial strain, minimum values of tr(A) near the front of the droplet √ are less than 2, suggesting that some polymers have a length that is less than their relaxed length in this region. (3) Along the droplet centreline and towards the rear of the droplet a region of positive extensional strain exists which increases the local polymer length there in a similar manner. Maximum polymer lengths in this region are around six times  the relaxed polymer length based on average values tr(A). of By t ≈ 0.7 Fig. 4 d and f shows that a wide fork has developed on the tail of both the simulated and experimental droplets. The experimental fork appears to be a little ‘deeper’ than that predicted at this time, however this conclusion is quite subjective given the resolution of the experimental image. By t ≈ 0.90 the fork on the tail of both droplets has grown and in the experimental image is now more clearly visible. Fig. 4 g and h (and also Fig. 3 a and b) indicates that its form is similar in both the experiments and simulations. To understand how this fork develops, we examine polymer behaviour. Comparing Fig. 4 a, d and g shows that the longer the droplet is in the contraction, the greater the extension of polymers within it, particularly those lining the internal surface of the droplet as described above. Towards the rear of the droplet, streamlines show (Fig. 4 b and e) that disperse phase fluid moves backwards along the outer radius of the droplet, and then inwards towards the droplet’s centreline. As polymers are advected by this flow, and as the polymers lying along the outer radius of the filament are stretched, there is a continual flux of extended polymers moving towards the rear tip of the droplet while it is within the contraction. How these extended polymers affect the droplet shape is determined largely by their orientation. The extended polymers that line the front and sides of the droplet while it is within the contraction are generally orientated tangential to the local droplet interface, as the dominant shear strains are in this direction in these regions. At the rear of the droplet this is not the case however. Polymers at the rear of the droplet have just been advected there from the sides of the droplet. As this advection has occurred quickly, and as there is only minimal shear strain present at the rear of the droplet, reorientation of the polymers has not occurred. Instead, polymers at the rear of the droplet are orientated in a mainly axial direction, which is almost perpendicular to the interface there (see Fig. 4d). As shown by Eq. (5), polymers exert an extensional stress on the fluid in the same direction as their orientation. As the polymers at the rear of the droplet are orientated normal to the droplet interface, the stress that they exert on the fluid pulls the rear interface of the droplet forward. This causes the forked tail observed in Fig. 4

74

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

Fig. 4. A comparison of simulation and experimental results for the viscoelastic droplet at various times. In frames (a), (d), (g), (i) and (k) the average polymer length,



tr(A),

is indicated by shading, and the average polymer orientation (scaled using the average length) is indicated by indigo ‘pointers’ (dark grey in print). The same shading scale for tr(A) is used in all of these frames. In frames (b) and (e) streamlines are shown in indigo (dark grey in print) over the simulated droplet shape. These streamlines are based

upon the fluid velocity taken relative to the leading tip of the droplet, thus showing the internal recirculation patterns within the droplets. Enlargements of the experimental images of Fig. 3 a are shown in frames (c), (f), (h), (j) and (l). In all images the fluid is flowing from left to right. (a) Polymer length (



(b) Relative velocity streamlines at t ≈ 0.42. (c) Droplet shape (experimental) at t = 0.49. (d) Polymer length (



streamlines at t ≈ 0.70. (f) Droplet shape (experimental) at t = 0.70. (g) Polymer length (

tr(A)) and orientation at t ≈ 0.42.

tr(A)) and orientation at t ≈ 0.70. (e) Relative velocity

tr(A)) and orientation at t ≈ 0.96. (h) Droplet shape (experimental) at t = 0.90.

(i) Polymer length and orientation at t ≈ 1.12. (j) Droplet shape (experimental) at t = 1.10. (k) Polymer length and orientation at t ≈ 1.32. (l) Droplet shape (experimental) at t = 1.32. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

d and g to develop. Further, the longer the droplet is in the contraction, the greater the number of polymers, extended by a greater amount, are advected into the rear region of the droplet. Consequently, the larger the fork that develops. The images at t ≈ 0.7 and t ≈ 0.9 support this conclusion. Fig. 4 j and l shows that as the real droplet continues through the contraction, its tail fork continues to grow, eventually entraining a thin cylinder of continuous phase fluid towards the rear of the droplet at t = 1.10. This fluid forms the continuous phase ‘packet’ encapsulated within the main droplet at t = 1.32. The simulated droplet develops in much the same way, however entrainment of continuous phase fluid does not occur at t ≈ 1.12, and in the final timeframe of Fig. 4 k only a very small amount of continuous phase fluid exists within the main droplet. The overall droplet shapes and tail forms predicted by the simulations are very similar to those observed experimentally however. As shown by Fig. 4 i, by the time the simulated droplet begins to exit the contraction, internal recirculation within the droplet has not only moved extended polymers from the sides of the droplet to its rear, but has also moved a thin stream of extended polymers from the rear of the droplet forward and along its centreline. Close

examination of the figure shows that towards the rear of the droplet and near the apex of the tail fork, these polymers are orientated in a direction which is slightly away from the axial direction. This orientation results in a stress on the droplet interface which acts to keep the two halves of the tail fork separate. The corresponding experimental image of Fig. 4 j shows that in reality this stress is indeed able to keep the two halves of the tail fork apart (a rotated fork in reality), and encapsulation of continuous phase fluid results. In the simulated case significant entrainment does not occur. Experience suggests that this is because the numerical method is not able to resolve flow details on length scales that are similar to or smaller than the computational cell size. Instead, coallescence occurs in the simulations at the rear of the droplet when the radius of the entrained stream is of the order of the computational cell size. The very small amount of continuous phase fluid that exists within the simulated droplet in the final timeframe of Fig. 4 k is formed as the forked tail of the droplet exits the contraction. Given the small size of these fluid packets relative to the computational cell size, their existence should be treated with caution however. Encouragingly, even though the entrainment and encapsulation process is not resolved by the simulations, the length of the

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

75

5. Conclusions

Fig. 5. Polymer length and orientation in a simulated viscoelastic droplet which has an increased polymer viscosity of p,d = 0.256. As in Fig. 4, the shading represents the average polymer length,



tr(A), and the average polymer orientation (scaled

using the average length)  is indicated by indigo ‘pointers’ (dark grey in print). The scale of Fig. 4 a is used for

tr(A). (a) t ≈ 0.96. (b) t ≈ 1.12. (c) t ≈ 1.32. (For inter-

pretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

entrained continuous phase fluid stream in the experimental image of Fig. 4 j is similar to the length of the central extended polymer stream observed in the simulation image of Fig. 4 j. This suggests that both the movement and extension of polymers within the droplet are captured by the physical and numerical models employed in the simulations. In order to accurately capture the entrainment and encapsulation process as well, a significantly finer computational mesh would be required; using such a mesh may be prohibitively expensive using currently available computing technology however. 4.4. Effect of increased polymer concentration To examine the sensitivity of droplet deformation to polymer concentration, the viscoelastic simulation of Fig. 4 was repeated using twice the previous polymer concentration (p,d = 0.256). Selected late time images are shown in Fig. 5. A comparison of the two simulations shows that the behaviour of the droplet is quite insensitive to polymer concentration. There does appear to be slightly more continuous phase fluid encapsulated within the droplet of Fig. 5 than the droplet of Fig. 4, however given the difficulties in resolving the entrainment process on this computational mesh, this difference is not significant. The maximum polymer extension that was experienced within tr(A) ≈ 120, and this extenthe simulated droplet of Fig. 4 was sion occurred as the leading tip of the droplet was just exiting the contraction, at t ≈ 0.9. For the higher polymer concentration  simulation of Fig. 5 the maximum extension experienced was tr(A) ≈ 100, also occurring around t ≈ 0.9. That droplet deformation is reasonably insensitive to polymer concentration (for the two non-zero concentrations considered), and that the maximum polymer extension experienced decreases as the polymer concentration is increased suggests that for a given contraction geometry, the stress field generated by the polymers is not sensitive to the polymer concentration. Therefore, to either enhance or suppress the encapsulation process other factors such as polymer relaxation time, contraction length or flow velocity should be altered. These factors will be examined in future investigations.

An existing two-dimensional finite volume immiscible fluid CFD algorithm has been extended to simulate flow with viscoelastic Oldroyd-B fluids. A polymer configuration transport equation is solved in each cell using two sequential updates per timestep: an advection update and a source update. The advection update uses fluid phase fluxes calculated using VOF techniques to calculate the flux of polymers crossing each cell face. In this way the algorithm recognises the discrete nature of the phase interface by not allowing polymer transport over immiscible phase boundaries. The source update is accomplished implicitly, partly to allow for finitely extensible dumbbells in the future, and partly to ensure solution stability in the limit of large ∇ ua components or small De. When calculating the elastic component of the total fluid stress, interpolation of the configuration tensor and polymer concentration product, rather than both separately, promotes stability around the phase interface. By comparing experimental results of a viscoelastic droplet passing through a thin three-dimensional planar contraction against simulation results calculated using a two-dimensional Cartesian contraction we have shown that the employed physical and numerical models are able to capture multiphase viscoelastic behaviour well. The forked tail that the experimental droplet adopts within the contraction is predicted by the simulations. The encapsulation of continuous phase fluid that occurs subsequently is not well resolved by the simulations, most probably due to inadequate mesh resolution. However, there is strong evidence that the polymer behaviour that causes the encapsulation is resolved, and that with the future availability of more powerful computers the encapsulation process will be accurately predicted. As well as allowing validation, the simulation results have been used to examine polymer behaviour within this microfluidic system. Shear strains near the front and sides of the droplet cause a thin shell of extended polymers to form around the droplet as it moves through the contraction. These polymers are advected by local recirculation currents towards the rear of the droplet, and also forwards from its rear along its centreline. The stress that these polymers exert cause the droplet to behave quite differently to a Newtonian droplet. In particular, the polymer stress allows a forked tail to develop on the droplet while it is within the contraction, which later results in the entrainment and eventually encapsulation of continuous phase fluid within the droplet as it leaves the contraction. Simulations suggest that the elastic stress that causes this entrainment is only weakly sensitive to polymer concentration within the droplet. Acknowledgement This research was supported by the Australian Research Council Grants Scheme. Appendix A. Numerical implementation of the advection algorithms In this section we provide pseudo-code flowcharts detailing the implementation of both the VOF and viscoelastic advection algorithms. Although the VOF advection algorithm has been described previously [28], as the two algorithms are linked in both structure and implementation we detail both for consistency and clarity. To simplify the VOF algorithm description we assume that the VOF function is stored on the standard (pressure) mesh, rather than on a mesh that is twice as fine (as used in real algorithm).

76

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

In the present context this difference is not significant: When a  value is required in a specific standard (pressure) cell in the real algorithm, it is calculated by volume averaging  values stored on the finer VOF mesh. Also, we describe only a single first-order time update: In the real algorithm two repeating first order updates are used to advance the solution in time (for all transient equations), resulting in second-order time accuracy overall [28]. Table A.1 defines the numerical notation used in our descriptions. Cell centred and face-centred variables are stored in one dimensional arrays, with the right most subscript indicating the variable’s location within the array. The ‘i’ subscripts are used to locate cell centred quantities within a total of I cells, while the ‘j’ subscripts are used to locate face-centred quantities within a total of J faces. Superscripts on the right indicate the time the variable’s value represents: ‘n’, ‘o’ and ‘a’ imply end-of-timestep, beginning-of-timestep, and timestep-averaged quantities, respectively. Coordinate directions are specified by ‘k’ indices, with the unit vectors ek defined in each coordinate direction. Here k = 1, 2 only. The orientation of each cell face is specified by a unit vector, n|j , directed normal to the face and in a positive coordinate direction. Integer functions are used to locate the neighbours of a cell or face. Fig. A.1illustrates their use. The functions ‘max’, ‘min’ and ‘sign’ are defined as in Fortran. The pseudo-code of Fig. A.2 details the VOF advection algorithm. As shown, the cell centred n |i and V |i variables are first initialised to their beginning-of-timestep values. The algorithm then loops through each coordinate direction (the order of this loop changes at each timestep), first calculating the disperse phase volumes fluxed over each cell boundary in the current direction, before updating the cell centred volume fractions using these fluxed amounts. After all coordinate directions have been considered, the calculation of n |i is complete.

Table A.1 Numerical variables used in the algorithm flowcharts Nomenclature Cell centred variable arrays |i p |i A|i B|i = p |i A|i

V |i V |i

Face-centred variable arrays u|j F|j = ua |j A|j t F |j = [o ua ]j A|j t

Fp |j = [op ua ]j A|j t FB |j = [op Ao ua ]j A|j t

A|j n|j

Definition Disperse phase volume fraction (VOF function) in cell i Polymer concentration in cell i Dumbbell configuration tensor in cell i Product of configuration tensor and polymer concentration in cell i Volume of cell i (constant) Total fluid volume contained within cell i which may not be equal to V |i during intermediate stages of the advection update Fluid velocity at face j in the direction of n|j Total fluid volume passing over face j during a timestep in the direction of n|j Total disperse phase volume passing over face j during a timestep in the direction of n|j Amount of polymer passing over face j during a timestep in the direction of n|j The tensor quantity of B passing over face j during a timestep in the direction of n|j Area of face j Unit vector perpendicular to face j

Integer functions used to locate neighbouring cells and faces j → i(j, s) Index of the cell adjacent to face j in the direction of sn|j from the centre of face j, where s = ±1 i → i(i, sk) Index of the cell next to cell i in the direction of sek from the centre of cell i, where k = 1, 2 and s = ±1 i → j(i, sk) Index of the face bordering cell i in the direction of sek from the centre of cell i, where k = 1, 2 and s = ±1 Global variables:

t

Timestep duration

Fig. A.1. A selection of computational cells and faces and some examples of the integer functions used to relate their indices.

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

77

Fig. A.2. Pseudo-code flowchart describing the VOF advection update performed at each timestep.

The pseudo-code of Fig. A.3 details the viscoelastic advection algorithm. The structure is essentially the same as that of the VOF advection algorithm: A first loop cycles through all cells initialising variables. A second loop then cycles through each of the coordinate directions. Within this loop the timestep-averaged configuration tensors [Ao ]j and fluxed amounts FB are calculated for each face orientated in the current direction using the procedure described in Section 2.2.2. The cell centred tensors B are then updated using these amounts. When both coordinate directions have been considered, the calculation of B is complete. The technique used to directionally split the advection of the viscoelastic tensor is analogous to the technique used to split the advection of the VOF function. As each coordinate direction is considered sequentially in the VOF advection procedure, intermediate  values must be bounded within the physical range of 0 ≤  ≤ 1 so that valid phase fluxed amounts F |j can be calculated during the next directional sweep. As shown in the final loop of Fig. A.2, the VOF algorithm achieves this by using the current total fluid volume in the cell V |i , rather than the total cell volume V |i (these differ between directional sweeps), to perform the update of |i [28]. The viscoelastic algorithm is similar. As shown in the final loop of Fig. A.3 it also uses the current total fluid volume in the cell V |i rather than the total cell volume V |i to perform the update of each B∗ |i . This ensures that intermediate values for |i and B∗ |i are sychronized (i.e. refer to the same phase volumes), allowing physically valid A profiles and hence fluxed polymer amounts FB |j to be calculated during the next directional sweep.

Fig. A.3. Pseudo-code flowchart describing the viscoelastic advection update performed at each timestep.

Appendix B. Mesh independence study The maximum droplet size in each coordinate direction as a function of time is shown in Fig. B.1 for simulations conducted using three computational meshes (and three corresponding timesteps). All results are qualitatively similar. The largest quantitative difference between results occurs around t ≈ 0.4 when the droplet is completely within the contraction and fully extended in the axial direction. At this time there is around a 10% difference in axial droplet size (Ymax ) between the 32 × 384 and 64 × 768 results, and around a 4% difference between the 64 × 768 and 128 × 1536

78

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79

Fig. B.4. Simulated droplet dimensions as a function of computational mesh size for the viscoelastic simulation of Fig. 4. Xmax (red, or grey in print) and Ymax (black) are the maximum nondimensional droplet dimensions in the x and y directions, respectively.

Fig. B.5. Simulated droplet shapes as a function of computational mesh size for the viscoelastic simulation shown in Fig. 4. Droplet forms at three different times are shown. Results produced using the 128 × 1536 mesh are given as solid black lines, using the 64 × 768 mesh as dotted orange lines (dotted light grey in print) and using the 32 × 384 mesh as dashed indigo lines (dashed dark grey in print). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

results. Significantly, the finest mesh results are bounded by the coarser mesh results during this period. In Fig. B.2 the droplet shapes produced using each of these meshes are compared at three different times. Consistent with the quantitative results of Fig. B.1, the droplet shapes produced on the finest 128 × 1536 mesh are closer to the shapes produced on the 64 × 768 mesh than those produced on the coarsest 32 × 384 mesh. Also, the finest mesh results are generally bounded by the two coarser mesh results. As would be expected in any CFD study, the finer the mesh, the finer the scale of the flow features predicted. For example, the forked tail which was experimentally observed around t ≈ 1.10 is most accurately predicted using the finest 128 × 1536 mesh. All results show some type of forked tail development however. References [1] M. Aboubacar, H. Matallah, H.R. Tamaddon-Jahromi, M.F. Webster, Numerical prediction of extensional flows in contraction geometries: hyrid finite volume/element method, Journal of Non-Newtonian Fluid Mechanics 104 (2002) 125–164. [2] M.A. Alves, F.T. Pinho, P.J. Oliveira, Visualisations of Boger fluid flows in a 4:1 square-square contraction, AIChE Journal 51 (11) (2005) 2908–2922. [3] P.-A. Auroux, D. Iossifidis, D.R. Reyes, A. Manz, Micro total analysis systems. 2. analytical standard operations and applications, Anal. Chem. 74 (2002) 2637–2652.

[4] J. Azaiez, R. Guénette, A. Aït-Kadi, Numerical simulation of viscoelastic flows through a planar contraction, Journal of Non-Newtonian Fluid Mechanics 62 (1996) 253–277. [5] D.V. Boger, Viscoelastic flows through contractions, Annual Review of Fluid Mechanics 19 (1987) 157–182. [6] D.V. Boger, R.J. Binnington, Experimental removal of the re-entrant corner singularity in tubular entry flows, Journal of Rheology 38 (2) (1994) 333– 349. [7] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modelling surface tension, Journal of Computational Physics 100 (1992) 335–354. [8] M.D. Chilcott, J.M. Rallison, Creeping flow of dilute polymer solutions past cylinders and spheres, Journal of Non-Newtonain Fluid Mechanics 29 (1988) 381–432. [9] J.J. Cooper-White, J.E. Fagan, V. Tirtaatmadja, D.R. Lester, D.V. Boger, Drop formation dynamics of constant low-viscosity, elastic fluids, Journal of NonNewtonian Fluid Mechanics 106 (2002) 29–59. [10] M.R. Davidson, D.J.E. Harvie, J.J. Cooper-White, Simulations of pendant drop formation of a viscoelastic liquid, Korea-Australia Rheology Journal 18 (2) (2006) 41–49, Jun. [11] K.V. Edmond, A.B. Schofield, M. Marquez, J.P. Rothstein, A.D. Dinsmore, Stable jets of viscoelastic fluids and self-assembled cylindrical capsules by hydrodynamic focusing, Langmuir 22 (2006) 9052–9056. [12] D.J.E. Harvie, M.R. Davidson, J.J. Cooper-White, Simulations of viscoelastic droplet deformation through a micro-fluidic contraction, pp. 69–78, in: Advances in Fluid Mechanics VI, WIT Transactions on Engineering Science. vol. 52, 2006. [13] D.J.E. Harvie, M.R. Davidson, J.J. Cooper-White, M.J. Rudman, A numerical parametric study of droplet deformation through a microfluidic contraction, ANZIAM J. 46 (E) (2005) C150–C166, Apr. [14] D.J.E. Harvie, M.R. Davidson, J.J. Cooper-White, M.J. Rudman, A parametric study of droplet deformation through a microfluidic contraction: Low viscosity Newtonian droplets, Chemical Engineering Science 61 (2006) 5149–5158. [15] D.J.E. Harvie, M.R. Davidson, J.J. Cooper-White, M.J. Rudman, A parametric study of droplet deformation through a microfluidic contraction: Shear thinning liquids, International Journal of Multiphase Flow 33 (5) (2007) 545–556. [16] D.J.E. Harvie, M.R. Davidson, M. Rudman, An analysis of parasitic current generation in volume of fluid simulations, Appl. Math Mod. 30 (10) (2006) 1056–1066. [17] D.J.E. Harvie, D.F. Fletcher, A hydrodynamic and thermodynamic simulation of droplet impacts on hot surfaces, part II: Validation and applications, International Journal of Heat and Mass Transfer 44 (14) (2001) 2643–2659. [18] D.D. Joseph, Fluid dynamics of viscoelastic liquids, Springer-Verlag, 1990. [19] R. Keunings, An algorithm for the simulation of transient viscoelastic flows with free surfaces, Journal of Computational Physics 62 (1986) 199–220. [20] D. Khismatullin, Y. Renardy, M. Renardy, Development and implementation of VOF-PROST for 3D viscoelastic liquid-liquid simulations, Journal of NonNewtonian Fluid Mechanics 140 (2006) 120–131. [21] J. Li, Y.Y. Renardy, M. Renardy, Numerical simulation of breakup of a viscous drop in simple shear flow through a volume-of-fluid method, Physics of Fluids 12 (2) (2000) 269–282. [22] M. Pasandideh-Fard, R. Bhola, S. Chandra, J. Mostaghimi, Deposition of tin droplets on a steel plate: Simulations and experiments, International Journal of Heat and Mass Transfer 41 (1998) 2929–2945. [23] S.B. Pillapakkam, P. Singh, A level-set method for computing solutions to viscoelastic two-phase flow, Journal of Computational Physics 174 (2001) 552–578. [24] S. Popinet, S. Zaleski, A front-tracking algorithm for accurate representation of surface tension, International Journal for Numerical Methods in Fluids 30 (1999) 775–793. [25] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes in FORTRAN: The art of scientific computing, Second Edition, Cambridge University Press, 1992. [26] S. Ramaswamy, L.G. Leal, The deformation of a viscoelastic drop subjected to steady uniaxial extensional flow of a newtonian fluid, Journal of NonNewtonain Fluid Mechanics 85 (1999) 127–163. [27] L.E. Rodd, J.J. Cooper-White, D.V. Boger, G.H. McKinley, Role of the elasticity number in the entry flow of dilute polymer solutions in micro-fabricated contraction geometries, Journal of Non-Newtonian Fluid Mechanics 143 (2007) 170–191. [28] M. Rudman, A volume-tracking method for incompressible multifluid flows with large density variations, International Journal for Numerical Methods in Fluids 28 (1998) 357–378. [29] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annual Review of Fluid Mechanics 31 (1999) 567–603. [30] J.A. Sethian, P. Smereka, Level set methods for fluid interfaces, Annual Review of Fluid Mechanics 35 (2003) 341–372. [31] H. Song, R.F. Ismagilov, Millisecond kinetics on a microfluidic chip using nanoliters of reagents, Journal of the Americal Chemical Society 125 (47) (2003) 14613–14619. [32] T.M. Squires, S.R. Quake, Microfluidics: fluid physics at the nanoliter scale, Reviews of Modern Physics 77 (3) (2005) 977–1026. [33] H.A. Stone, A.D. Stroock, A. Ajdari, Engineering flows in small devices: Microfluidics toward a lab-on-a-chip, Annual Review of Fluid Mechanics 36 (2004) 381–411.

D.J.E. Harvie et al. / J. Non-Newtonian Fluid Mech. 155 (2008) 67–79 [34] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, Journal of Computational Physics 114 (1994) 146–159. [35] R.I. Tanner, S. Xue, Computing transient flows with high elasticity, KoreaAustralia Rheology Journal 14 (4) (2002) 143–159, Dec. [36] V. Tirtaatmadja, G.H. McKinley, J.J. Cooper-White, Drop formation and breakup of low viscosity elastic fluids: Effects of molecular weight and concentration, Physics of Fluids 18 (2006) 043101. [37] M.F. Tome, J.L. Doricio, A. Castelo, J.A. Cuminato, S. McKee, Solving viscoelastic free surface flows of a second-order fluid using a marker-and-cell approach, International Journal for Numerical Methods in Fluids 53 (2007) 599–627. [38] M.F. Tome, L. Grossi, A. Castelo, J.A. Cuminato, S. McKee, K. Walters, Die-swell, splashing drop and a numerical technique for solving the Oldroyd-B model for axisymmetric free surface flows, Journal of Non-Newtonian Fluid Mechanics 141 (2007) 148–166.

79

[39] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of multiphase flow, Journal of Computational Physics 169 (2001) 708–759. [40] V.M. Ugaz, R.D. Elms, R.C. Lo, F.A. Shaikh, M.A. Burns, Microfabricated electrophoresis systems for DNA sequencing and genotyping applications: current technology and future directions, Phil. Trans. R. Soc. Lond. A 362 (2004) 1105–1129. [41] E. Verpoorte, Microfluidic chips for clinical and forensic analysis, Electrophoresis 23 (2002) 677–712. [42] S.A. White, A.D. Gotsis, D.G. Baird, Review of the entry flow problem: Experimental and numerical, Journal of Non-Newtonian Fluid Mechanics 24 (1987) 121–160. [43] D.L. Youngs, Time-dependent multimaterial flow with large fluid distortion, in: K. Morton, M. Baines (Eds.), Numerical Methods for Fluid Dynamics, Academic Press, 1982, pp. 273–285.