Two-dimensional study of drop deformation under simple shear for Oldroyd-B liquids T. Chinyoka a Y.Y. Renardy a,1 M. Renardy a D.B. Khismatullin b a Department of Mathematics, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061-0123, USA b Department of Biomedical Engineering, Duke University, Durham, NC 27708, USA Abstract The effect of viscoelasticity on the deformation of a circular drop suspended in a second liquid in shear is investigated with direct numerical simulations. A numerical algorithm based on the volume-of-fluid method for interface tracking is implemented in two dimensions with the Oldroyd-B constitutive model for viscoelastic liquids. The code is verified against a normal mode analysis for the stability of two-layer flow in a channel; theoretical growth rates are reproduced for the interface height, velocity and stress components. Drop simulations are performed for drop and matrix liquids of different viscosities and elasticities. A new feature is found for the case of equal viscosity, when the matrix liquid is highly elastic and surface tension is low; hook-like structures form at the drop tips. This is due to the growth of first normal stress differences that occur slightly above the front tip and below the back tip as the matrix elasticity increases above a threshold value. Key words: Droplet deformation; Viscoelastic fluids; Direct numerical simulation 1 Introduction The investigation of a single drop deforming under simple shear is a model problem in emulsion technology [1, 2]. In shear mixing, when the shear rate is large enough to overcome the interfacial tension between the fluids, initially large drops of the dispersed fluid deform and rupture into smaller drops [3]. 1 Corresponding author. E-mail address: renardyy@math.vt.edu (Y. Renardy) Preprint submitted to Elsevier Science 12 May 2005 The ultimate sizes and distribution of the droplets play crucial roles in the physical and rheological properties of an emulsion. New products manufactured from polymer blends combine favorable characteristics of the constituent viscoelastic fluids. In this paper, we consider a horizontal array of dispersed drops at low concentration. When the constituent fluids are Newtonian, numerical simulations have been performed in three dimensions with a boundary integral method for Stokes flow [4, 5], a finite element method applied to Stokes flow [6], and with the volume-of-fluid method for the inclusion of inertia [7]. The latter approach was also used for yield-stress fluids [8]. Viscoelastic effects have mostly been simulated in 2D or axisymmetric situations; e.g., a time-periodic extensional flow at non-zero Reynolds number was treated with a front-tracking finitedifference method in [9], a level-set method was used for the deformation of a circular ‘drop’ under shear [10], and a diffuse-interface model was used in [11]. Sec. 2 is an outline of the dimensionless continuum equations governing the flow system for Oldroyd-B liquids. Sec. 3 is a description of the viscoelastic algorithm based on the volume of fluid (VOF) approach. Results are compared with the Newtonian boundary-integral results of [12], and theoretical stability analyses of two-layer channel flow of Oldroyd-B liquids. Sec. 4.1 is a study of 2D drop deformation for parameters similar to [10]. 2 Mathematical Modeling Fig. 1 is a sketch of the initial condition of the model problem. The boundary conditions are that the walls at z = 0, Lz move at speeds UB = −UT and UT , respectively. Periodic boundary conditions are posed at x = 0, Lx . The boundaries are taken far enough away so that the drops do not interact with neighbors or walls. We use an initial drop radius of a = Lz /8 with Lz = 1 throughout, while Lx is usually 16a, except in Figure 10 below, where Lx = 32a. The initial conditions we use for velocities and the viscoelastic stresses are a uniform shear flow, and either (i) the stresses which would prevail in a steady flow of this shear rate, or (ii) zero stresses. Comparisons showed that the choices of the initial stresses made little difference for those cases resulting in a stationary drop shape; at higher rates of stretching, the difference can become substantial, as we shall see below. The conditions at the interface are that the velocity and shear stress are continuous, and the jump in the normal stress is balanced by surface tension. The initial condition is not required to satisfy these conditions but they are satisfied as the computation progresses. In the volume-of-fluid approach, rather than explicitly tracking the interface, a color function C is used to track the 2 Fig. 1. Initial condition for the drop-in-a-flow system. position of the liquids: C(t, x, z) = ½ 1 in fluid 1 (drop), 0 in fluid 2 (matrix). (1) This function is advected by the flow. The position of the interface is defined as the surface across which C jumps in value. Densities are denoted ρi , solvent viscosities ηsi , polymeric viscosities ηpi , total viscosities ηi = ηsi + ηpi , and relaxation times λi , where i refers to fluid 1 or 2. The horizontal channel length Lx and vertical plate separation Lz are taken large enough compared with the initial drop radius a, so that the drop is effectively in an infinite medium at constant shear rate γ̇ = (UT − UB )/Lz far enough away. The interfacial tension is denoted σ. The dimensionless parameters are the capillary number, Reynolds number, viscosity ratio, and Deborah numbers, denoted respectively by: Ca = aγ̇η2 /σ, R = ρ2 γ̇a2 /η2 , m = η1 /η2 , Dei = λi γ̇. (2) The droplet-matrix system is governed by the incompressible momentum equations: ρ( ∇ · u = 0, (3) ∂u + u · ∇u) = ∇ · ~τ + F, ∂t (4) where the total stress tensor is ~τ = −pI + T + ηS, I is the unit tensor, T is the extra stress tensor, and S = [∇u + (∇u)T ] is the deformation rate tensor. In the continuous surface force (CSF) model of [21], the body force F includes 3 the interfacial tension force Fs = σκns δs . (5) Here, κ is the total interface curvature, ns is the normal to the interface, and δs is the delta function at the interface. The Oldroyd-B constitutive equation for the extra stress tensor T is ∇ λ T +T = ηp S, (6) ∇ where the upper convected time derivative T is defined as, ∇ T= 3 ∂T + (u · ∇)T − (∇u)T − T(∇u)T . ∂t (7) Numerical algorithm The Newtonian counterpart of our algorithm was originally developed in [13, 14, 15] for high Reynolds numbers, and modified in [16, 7] to handle low Reynolds numbers. In this section, we sketch this framework and focus on the viscoelastic component [17]. The discretization of the governing equations is based on a rectangular Cartesian mesh and staggered grid, on which finite-differences are taken. The grid places the pressure p, volume fraction function C and diagonal stress components T11 and T22 at the center of the rectangular cell; the velocity components u and v are on its edges, and non-diagonal stress component T12 = T21 is at a corner. In the discretized problem, C represents the volume fraction of fluid 1 in a grid cell; thus, when a cell is cut by the interface we have 0 < C < 1. The density ρ and viscosity η of each fluid is constant. In cells cut by the interface, the physical properties are interpolated: ρ = ρ1 C +ρ2 (1−C), η = η1 C +η2 (1−C), and λ = λ1 C + λ2 (1 − C). At each timestep, the interface is reconstructed by a piecewise linear interface calculation (PLIC) technique [18, 19, 20]. The interface moves according to a Lagrangian method: xn+1 = xn + u(∆t). The VOF-CSF formulation of [21] uses a mollified color function for the calculation of (5). The mollified color function c is obtained by taking the con4 volution of the volume fraction function C with a kernel K: c(t, x) = Z C(t, x′ )K(x′ − x, ε)dx′ . Ω c is a smooth function that transitions from 0 to 1 in a narrow region localized around the interface cells. This facilitates the calculation of ns = ∇c/|∇c| and δs = |∇c| . Suppose the solution at the nth timestep is known. The numerical solution for the next timestep begins with the projection method for the momentum equation (4) which decouples the pressure computation [22]. Let u∗ denote the intermediate velocity field which needs to be calculated in the projection method. The semi-implicit time-integration scheme of [16] is used to decouple the computation of the components of u∗ = (u∗ , v ∗ ): u∗ − u(n) (n) (n) = −(u(n) u(n) uz ) x +v ∆t à ! ∂ 1 ∂ (n) ∗ ∗ (n) , (2ηux + T11 ) + (T12 + η(uz + vx )) + F1 + ρ ∂x ∂z v ∗ − v (n) = −(u(n) vx(n) + v (n) vz(n) ) ∆t ! à ∂ 1 ∂ (n) (n) ∗ ∗ . (η(uz + vx ) + T12 ) + (T22 + 2ηvz ) + F2 + ρ ∂x ∂z This scheme allows for the decoupling of variables, which is an advantage over a fully implicit scheme, and has an associated error term of O((η∆t/ρ)2 ). The implicit terms in the u∗ -equation are then à " à ! à à !) ∆t ∂ ∂ ∂ ∂ I− 2η + η ρ ∂x ∂x ∂z ∂z !#! u∗ . The operator on u∗ factorizes; ( à ∂ ∆t ∂ 2η I− ρ ∂x ∂x !)( ∆t ∂ ∂ I− η ρ ∂z ∂z u∗ = explicit terms. (8) An analogous system is solved for v ∗ . Hence, the solution procedure for u∗ reduces to inversions of tridiagonal matrices. The incompressibility condition ∇ · u(n) = ∇ · u(n+1) = 0, and ∇p un+1 − u∗ =− , ∆t ρ (9) 5 lead to the pressure equation: à ∇p ∇· ρ ! = ∇ · u∗ . ∆t (10) This is solved with a multigrid method. The updated solution u(n+1) is then found with (9). The constitutive equation (6) is treated with an analogous semi-implicit time integration scheme. The terms which involve spatial derivatives of the extra stress tensor are taken implicitly, as is T: ∂ ∂ T(n+1) − T(n) + (u(n) + v (n) )T(n+1) + T(n+1) . λ ∆t ∂x ∂z ! à The equations for T11 , T12 , and T22 decouple, and a factorization as in (8) leads to à ∂ λ u(n) ∆t (λ + ∆t) 1 + λ + ∆t ∂x !à ! ∂ λ v (n) ∆t 1+ T(n+1) λ + ∆t ∂z (11) = explicit terms, with an associated error O ((λkukmax ∆t)2 /(λ + ∆t)). The solution procedure reduces to inversions of tridiagonal matrices. The interpolation of fluid properties in the grid cells surrounding the interface can lead to a local violation of the positive definiteness of the configuration tensor. We have observed instabilities resulting from this [23], and we suppress these by adding a multiple of the identity matrix to the stress as needed, to ensure positive definiteness of the tensor T + (ηp /λ)I. 3.1 Code validation 3.1.1 Newtonian drop deformation We investigate the parameter range of [12]; the height of their computational domain is denoted 2H, the top wall moves to the right with velocity U and the lower wall moves with velocity −U . Lengths are non-dimensionalized using half the height; hence their square computational domain is [0, 2]×[0, 2]. Their dimensionless drop radius is denoted a/H = 0.5. The drop viscosity is mη where η = ηs is the matrix liquid viscosity. The shear rate is U/H; hence, the Reynolds number is R = ρU 2a/η and the capillary number is Ca = U η/σ. The cases Ca = 0.2 and Ca = 0.4, m = 1 and m = 10 and R = 1, 10, 50 and R = 100 are addressed in Figs. 3 and 5 of [12]. Drop shapes and transient data are reproduced [17]: Fig. 2 illustrates our results for the evolution of drop deformation in time for the same parameters as those of [12]. The deformation is, in the usual way, represented by a deformation parameter D defined as 6 D = (Rmax − Rmin )/(Rmax + Rmin ), where Rmax and Rmin are the largest and shortest distances of the interface from the drop center, respectively. In particular, for the case m = 1, R = 1 and both Ca = 0.2 & Ca = 0.4, Fig. 2d is a composite of Fig. 3a of Sheth & Pozrikidis [12] and ours, showing that our results accurately reproduce them. We also have agreement at other parameter values: Fig. 3 shows the drop shapes for all the relevant parameter values in [12], the two capillary numbers, the two viscosity ratios and the four Reynolds numbers. As mentioned before on the evolution of drop deformation, the steady state shapes in Fig. 3 (cases R = 1, 10) accurately depict the corresponding shapes in Fig. 5 of [12]. The cases R = 50, 100 continue to deform and the corresponding qualitative features shown in Fig. 3, such as the dumbbell shapes, accurately portray those in Fig. 5 of [12]. 0.7 0.8 Re=1 Re=10 Re=50 Re=100 0.6 Re=1 Re=10 Re=50 Re=100 0.7 0.6 0.5 0.5 D D 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 0.5 1 1.5 2 2.5 3 time 0 0.5 1 1.5 2 time (a) (b) (c) (d) Fig. 2. Evolution of drop deformation for m = 1; (a) Ca = 0.2, (b) Ca = 0.4, (c) Fig. 3a of [12], and (d) is a superposition of (c) and our work. 7 2.5 3 Re=1 Re=10 Re=50 Re=100 1.8 1.4 1.4 1.2 1.2 1 1 z 1.6 z 1.6 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 (a) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 1.2 1.4 1.6 1.8 x Re=1 Re=10 Re=50 Re=100 1.8 1.6 1.6 Re=1 Re=10 Re=50 Re=100 1.4 1.4 1.2 1.2 1 z z 0 (b) 1.8 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 (c) Re=1 Re=10 Re=50 Re=100 1.8 0 0 0.2 0.4 0.6 0.8 1 x 1.2 1.4 1.6 1.8 (d) 0 0.2 0.4 0.6 0.8 1 x Fig. 3. Drop shapes at t = 3γ̇ −1 (a) Ca = 0.2, m = 1, (b) Ca = 0.2, m = 10, (c) Ca = 0.4, m = 1 & (d) Ca = 0.4, m = 10. 3.1.2 Viscoelastic two-layer channel flow The interfacial stability of a channel flow consisting of a highly elastic OldroydB fluid superposed on silicone oil is investigated experimentally and theoretically in [24, 25]. The fluid properties in [24] are used to compare the growth rates predicted from the linear stability analysis of two-layer Couette flow with our direct simulations. The dimensionless parameters are: undisturbed interface height Li = L∗i /Lz = 0.73, m = 0.95, ρ1 /ρ2 = 1.003, Reynolds number R1 = Ui L∗i ρ1 /η1 = 0.85, Weissenberg number W1 = 0, W2 = Ui λ2 /L∗int = 0.14, interfacial tension parameter σ/(η2 Ui ) = 88, gL∗ /Ui2 = 3, retardation parameter β1 = 0, β2 = η2s /η2 = 0.84. * denotes dimensional quantities, 1 refers to the lower liquid, 2 refers to the upper liquid, and i refers to the interface. In the initial condition for the numerical simulation, the interfacial eigenmode from the stability analysis is superposed on the undisturbed two-layer solution. The computational domain is [0, 2π/α]×[0,1], where α denotes the wavenumber of the disturbance. Fig. 4 shows agreement between the computed (dash) and theoretical decay rates at α = 6. Similar evolutions for the maxima of |u|, |T11 |, |T12 |, and |T22 | are obtained and details are given in [17]. 8 −6.85 −6.9 −6.95 −7 logh −7.05 −7.1 −7.15 −7.2 −7.25 −7.3 −7.35 0 0.01 0.02 0.03 0.04 0.05 0.06 time Fig. 4. Evolution of ln h(t), where h(t) denotes the difference of the interface height at time t with the initial undisturbed interface height. Numerical (- - -) vs theoretical (—). 9 4 Numerical results on drop deformation for viscoelastic liquids A Newtonian drop deforming in a viscoelastic fluid and a viscoelastic drop in a Newtonian liquid are both addressed in [10]. We begin with their parameter range. With R = 0.0003, Ca = 0.24, De = 0.4, ηs = ηp and equal viscosities (η1 /η2 = 1), [10] finds that a Newtonian drop in a viscoelastic liquid reaches a stationary state with D ≈ 0.48. We find D ≈ 0.3. Fig. 5 shows our results for R = 0.3 for various drop-matrix liquid combinations including the case addressed by [10]. For the two-letter notation in Fig. 5, the first letter gives the droplet type and the second represents the matrix fluid; e.g., NV stands for a Newtonian drop in a viscoelastic matrix. Comparison with results at smaller R assure that inertial effects are negligible; we use the larger R to increase ∆t and reduce computational effort. Spatial and temporal refinements on the transient evolution for the VV case give the same results: ∆x = a/8n, ∆z = a/16n, n = 1, 2, 4, and ∆t = 0.001γ̇ −1 , 0.0005γ̇ −1 . An alternative initial condition with zero stress gives the same results. 0.35 0.3 0.25 0.2 D 0.65 0.6 0.15 z 0.55 0.5 0.1 NN VN NV VV 0.05 0.45 0.4 0.35 0 0 0.5 1 1.5 2 2.5 Time 3 0.3 0.7 0.8 0.9 1 1.1 1.2 1.3 x (a) (b) o System D φmax NN 0.29 32.3 VN 0.28 31.2 NV 0.26 28.2 VV 0.26 28.2 Fig. 5. Ca = 0.24, R = 0.3, De = 0.4. (a) Temporal evolution of D for the NN, VN, NV and VV systems, (b) Stationary shapes. The table gives data corresponding to (b) at t = 3γ̇ −1 . The table in Fig. 5 shows that the NN case settles to a higher angle than the NV. In fact, ‘the orientation of the drop towards the direction of the flow is significantly enhanced in the viscoelastic case than in the Newtonian case’ for both the small deformation theory of [26], and the larger deformations 10 addressed in the experimental work of [27]. We find that D is higher in VN than in NV or VV, in line with the conclusion in [11] that for small De (at correspondingly low Ca) ‘viscoelasticity in the matrix suppresses deformation more than if it is in the drop’. Fig. 6 shows our results for Ca = 0.6 with the other parameters as before. If the matrix fluid is Newtonian then the drop continues to deform (and in a 3D investigation may break), but settles to a stationary state when the matrix fluid is viscoelastic. In particular, Fig. 6 (b) shows that the viscoelastic drop in the VN system (dashed) has the characteristic dumbbell shape showing continued deformation whereas the Newtonian drop in the NV system (dotdashed) has an almost elliptical shape. [10] reports the opposite scenario, that the NV case continues to deform and the VN case settles. 0.8 0.7 0.6 0.5 D 0.65 0.4 0.6 0.55 z 0.3 0.5 0.2 0.45 NN VN NV VV 0.1 0.4 0.35 0 0 1 2 3 4 5 6 7 8 9 time 10 0.3 0.7 0.8 0.9 1 1.1 1.2 1.3 x (a) (b) System D φmax NN 0.77 14.8 VN 0.72 15.7 NV 0.60 14.7 VV 0.55 15.8 Fig. 6. Ca = 0.6, De = 0.4, R = 0.3. (a) Temporal evolution of D for NN, VN, NV and VV systems; (b) Drop shapes at t = 10γ̇ −1 . The table gives data on (b). We next explain our results in terms of the first normal stress difference N1 = T11 − T22 , the deformation angles and shear strengths. From Fig. 7, we notice that the high N1 in the NV system acts to force the Newtonian drop to align more with the direction of the flow, leading to the lower deformation angle than that in the VN system. This means that the (lower angle) Newtonian drop in the NV case is subject to much less shear than the (more inclined) viscoelastic drop in the VN case explaining why it (the Newtonian drop) would be less deformed than the viscoelastic drop. This argument would however not explain the differences in D for the (NN, VN) or the (VV, NV) systems, i.e. 11 0.65 0.65 0.018 0 0.6 0.6 0 0.016 0.55 z z 0.045 0.009 0.55 0.018 −0.008 0.5 0.018 0.5 −0.008 0.45 0.45 0.016 0.018 0.4 0.35 0.3 0.7 0 0.009 0.045 0 0.4 0.35 0.018 0.3 0.8 0.9 1 1.1 1.2 1.3 0.7 0.8 0.9 1 x x (a) (b) 1.1 1.2 1.3 Fig. 7. Contours of N1 for two of the systems shown in Fig. 6 (b): (a) VN, (b) NV. systems in which the matrix is kept constant but the drop is changed between the two possible choices. In this case we mention the usual expectation that the viscoelasticity inside the drop tends to reduce drop deformation. These expectations would dictate that the VN system should be less deformed than the NN and similarly that the VV system should be less deformed than the NV, as evident in Fig. 6. The case of a viscoelastic drop in a Newtonian fluid at R = 0.0003, Ca = 60, De = 8.0, ηs = ηp , m = 1, is investigated in [10]; at t = 10γ̇ −1 , they find D = 0.8152, and conclude that the drop would eventually break up. We find that the drop deforms into a long and ever-thinning ‘sheet’ or finger shown in Fig. 8 (a). The droplets at the ends of the drop are sub-grid and mesh refinements confirm that the drop is a single sheet. Similar sheets result from the NN and VV systems. The NV case shown in Fig. 8 (b) has a unique structure and is addressed in the next section. (a) 0.8 0.7 z 0.6 0.5 0.4 0.3 0.2 −2 −1 0 1 2 3 4 x (b) 0.65 0.6 z 0.55 0.5 0.45 0.4 0.35 0.3 −2 −1 0 1 2 3 4 x Fig. 8. Ca = 60, De = 8, R = 0.3, ηs = ηp , m = 1. Drop shapes at t = 10γ̇ −1 . (a) VN system, (b) NV system. 4.1 Development of hooks on a thinning drop We consider the case of a Newtonian drop in a viscoelastic liquid with R = 0.3, Ca = 60, ηs = ηp , m = 1, and examine the effect of increasing the elasticity 12 of the matrix liquid. For the following results, the initial stresses in the matrix fluid are those of a uniform shear flow. Hence, our results correspond to a situation where a drop is injected into a pre-existing shear flow. If, on the other hand, the flow is started with the drop already in place, then zero initial stresses would be appropriate for the initial condition. In the situation considered here, this makes a significant difference, because the drop is undergoing significant stretching during the time the elastic stresses take to build up. For zero initial stresses, the drop therefore stretches more slowly. Moreover, the vortices and hooks shown below were not found to develop in this case. Fig. 9 shows velocity vectors at t = 3γ̇ −1 as the Deborah number of the matrix liquid increases. We also tabulate D and φmax . The plots have been cropped in order to show details in the drop region; the computational domain is much larger. Note that the flow becomes simple shear at some distance away from the drop, as is evident in Figure 10. Evidently, viscoelasticity in the matrix fluid enhances deformation in this case. In particular, the vortices around the drop tips gradually appear and the hooks develop. Table 1 D and φmax corresponding to Fig. 9. De Rmax Rmin D φmax 0.01 0.41 0.04 0.83 16.9 0.1 0.41 0.04 0.83 16.9 0.5 0.41 0.04 0.84 16.2 2 0.47 0.03 0.88 15.0 4 0.56 0.02 0.92 12.4 8 0.72 0.02 0.95 8.9 An understanding of these results is achieved by looking at the precise positioning of the viscoelastic stresses. In particular, we calculate the first normal stress difference N1 , and how this evolves with De. Fig. 11 shows contours of N1 , at six values of the Deborah number. In (a), N1 is extremely small compared with (b), and the maximal areas are at the drop tips. At the higher De, the regions of maximal N1 move away from the tip, to above the front tip and below the back tip. This forces the drop to align more with the direction of flow, thus reducing φmax as shown in Fig. 9 f. Fig. 12 shows contours of the pressure field. At low matrix elasticity, Fig. 12 a shows that the pressure is maximal inside the drop at the tips. This location of high pressure regions and the corresponding positioning of the high N1 regions just at the drop tips result in mostly extension. This severely diminishes the mechanisms for the creation of vortices around the tips, and hence the hooks do not appear. On the other hand, Fig. 12 b shows that the high (and positive) 13 0.65 0.6 0.5 0.4 (a) 0.35 0.3 1 1.7 X−Axis 0.65 0.6 0.5 0.4 (b) 0.35 0.3 1 1.7 X−Axis 0.65 0.6 0.5 0.4 (c) 0.35 0.3 1 1.7 X−Axis 0.65 0.6 0.5 0.4 (d) 0.35 0.3 1 1.7 X−Axis 0.65 0.6 0.5 0.4 (e) 0.35 0.3 1 1.75 X−Axis 0.65 0.6 0.5 0.4 (f) 0.35 0.3 1 1.7 X−Axis Fig. 9. Ca = 60, R = 0.3, m = 1. Newtonian drop in viscoelastic matrix with η s = ηp at t=3γ̇ −1 . Velocity vector plots for (a) De = 0.01, (b) De = 0.1, (c) De = 0.5, (d) De = 2, (e) De = 4, (f) De = 8. pressure regions lie outside the drop, below the tip of the front hook and above the tip of the lower-end hook. The lowest (and negative) pressure regions lie inside the ‘neck’ regions of the hooks. Comparing Fig. 12 to the velocity vectors for the same parameters, say Fig. 9 f, we notice that the high pressure regions force the fluid back towards the low pressure zones at the same time as the correspondingly high N1 is rotating the drop towards the direction of flow. These effects combine to create the vortices around the drop tips and hence 14 1 0 0 1 2 3 4 X−Axis Fig. 10. The results for Fig. 9 f are obtained in the larger computational domain, for which the boundaries are taken sufficiently far away from the drop. also the hooks shown in Fig. 9 f. For comparison, Figure 13 shows a drop for De = 4, t = 10γ̇ −1 , and otherwise the same parameters as Figure 9. We see that no vortices have formed. 15 0.5 0.65 0.5 0 1 0 0.5 0.6 0.5 0.55 0 0 0 0.5 z 0 0.5 0.5 0.5 0.5 0.5 0 0 0 0.45 0 0.5 0.5 0.5 0.4 0 1 0.35 0 values * e−03 0.5 0.5 0.5 (a) 0.3 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.4 1.5 x 0.65 3.3 0.6 4.4 1.1 0.55 4.4 z 0.5 4.4 1.1 3.3 3.3 1.1 4.4 4.4 3.3 0.45 1.1 3.3 0.4 4.4 1.1 0.35 (b) values * e−03 0.3 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 x 2.7 0.65 2.4 2.7 0.6 2.4 0.3 2.4 0.55 0.3 z 2.7 0.5 2.4 2.7 0.3 0.45 2.4 0.3 2.7 0.4 2.4 2.7 0.35 values * e−02 (c) 0.3 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 x 0.1 0.65 0.18 0.1 0.6 0.1 0.1 0.02 0.55 z 0.02 0.1 0.5 0.1 0.45 0.1 0.18 0.02 0.4 0.35 0.1 (d) 0.3 0.2 0.4 0.6 0.8 1 0.1 1.2 1.4 1.6 1.8 x 0.65 0.3 0.3 0.54 0.12 0.6 0.12 0 0.55 z 0.12 0.5 0.12 0.12 0.45 0 0.54 0.3 0.4 0.3 0.12 0.35 (e) 0.3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 0.48 0.65 0.12 0.48 0.6 z 0 0.48 0.12 0.55 0.12 0.5 0.12 0.12 0.12 0.12 0.45 0.48 1.2 0 0.4 0.48 0.12 0.35 0.48 (f) 0.3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x Fig. 11. Contour plots of N1 at Ca = 60 t = 3γ̇ −1 , (a) De = 0.01, (b) De = 0.1, (c)De = 0.5, (d) De = 2, (e) De = 4, (f) De = 8. Contour values should be multiplied by 10−3 for Figures a and b, and by 10−2 for Figure c. 16 0 0 0.65 5 0.6 15 5 5 0.55 −5 5 z 5 0.5 0 5 5 0.45 0 5 5 −5 0.4 15 5 0.35 (a) values * e−04 0.3 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 x 0.65 5 0 0.6 0 5 0.55 0 z 0 −5 0 5 0 0.5 0 5 −5 0.45 0 0 5 0.4 0 15 5 0 0.35 values * e−04 0 (b) 0.3 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 x values * e−03 0.3 0 0.65 1.2 −1.8 3 0.6 −3 0 −1.8 0.3 0.55 −1.8 z 0.3 1.2 1.2 −3 −3 0 0.5 −1.8 0.45 −1.8 −3 0.4 3 −1.8 1.2 0.35 (c) 0 0.3 0.3 0 0.3 0.4 0.6 0.8 1 1.2 1.4 1.6 x 5 0.8 −1 −5 −1 0.7 −5 5 −15 −1 −15−5 5 0.6 z 5 −1 0.5 15 −15 −5 5 −15 15 −15 0.4 −5 5 −1 −5 5 −1 0.3 values * e−03 0.2 (d) 5 0.4 0.2 0.6 0.8 1 1.2 1.4 1.6 1.8 x values * e−02 0 0.8 0 2 −2 −2 0.7 −3 0.6 −3 2 −2 −3 −4 0 26 −20 −3 z −4 −4 2 2 4 −2 0.5 4 −2 2 2 4 0.4 2 4 6 0.3 −2 0 −4 −3 −2 0 −3 −3 −2 2 −2 0 2 0.2 (e) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 0.8 0 5 0 0.7 −5 5 −5 −5−5 −5 0.6 z 5 0 15 10 5 10 −10 0 0.5 −5 0.4 10 10 15 5 −10 −5 −5 0 5 −5 0 0.3 values * e−02 (f) −5 5 0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 01.6 1.8 x Fig. 12. Contour plots of p corresponding to Fig. 11; (a) De = 0.01, (b) De = 0.1, (c) De = 0.5, (d) De = 2, (e) De = 4, (f) De = 8. Contour values should be multiplied by 10−4 for Figures a and b, by 10−3 for Figures c and d, and by 10−2 for Figures e and f. 17 1 0 0 1 2 X−Axis Fig. 13. Newtonian drop in viscoelastic matrix for the parameters corresponding to Figure 9, De = 4, t = 10γ̇ −1 , with elastic stresses initially equal to zero. 18 5 Conclusion A two-dimensional viscoelastic volume-of-fluid continuous-surface-force technique has been implemented. This algorithm is applicable to the numerical simulation of flows of two immiscible Oldroyd-B liquids. We investigate the stretching of a Newtonian drop in a viscoelastic matrix. At low capillary numbers, matrix elasticity rotates the drop into the flow direction and slightly decreases drop deformation. At high capillary number and high elasticity, the evolution depends on the initial stress conditions: for initially established shear flow, vortices form ahead of the ends of the stretching drop, leading to the formation of hooks. Acknowledgements This work is supported by NSF-CTS 0090381, NCSA, NSF-DMS 0456086 and NSF-DMS 0405810. We thank Dr Jie Li (Cambridge) for discussions, and the reviewers for insightful suggestions. References [1] D. I. Bigio, C. R. Marks, R. V. Calabrese, Predicting drop breakup in complex flows from model flow experiments, International Polymer Processing, XIII2 (1998) 192-198. [2] S. Guido, F. Greco, Dynamics of a liquid drop in a flowing immiscible liquid, Rheology Reviews 2004, eds D. M. Binding, K. Walters, British Society of Rheology (2004) 99-142. [3] T.G. Mason, J. Bibette, Shear rupturing of droplets in complex fluids, Langmuir 13 (1997) 4600 - 4613. [4] V. Cristini, J. Blawzdziewicz, M. Loewenberg, An adaptive mesh algorithm for evolving surfaces: simulations of drop breakup and coalescence, J. Comp. Phys. 168 (2001) 445-463. [5] V. Cristini, R. Hooper, C. W. Macosko, M. Simeone, S. Guido, A numerical and experimental investigation of lamellar blend morphologies, Ind. Eng. Chem. Res. 41 (2002) 6305-6311. [6] R. Hooper, V. Cristini, S. Shakya, J. Lowengrub, C. W. Macosko, J. J. Derby, Modeling multiphase flows using a novel 3D adaptive remeshing algorithm, in ‘Computational Methods in Multiphase Flow’, eds C. A. Brebbia, H. Power, Advances in Fluid Mechanics 29 (2001) Wessex Institute of Technology Press, UK. [7] J. Li, Y. Renardy, M. Renardy, Numerical simulation of breakup of a 19 [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] viscous drop in simple shear flow with a volume-of-fluid method, Phys. Fluids. Vol. 12 (2) (2000) 269-282. J. Li, Y. Renardy, M. Renardy, Shear induced rupturing of a viscous drop in a Bingham liquid, J. Non-Newtonian Fluid Mech. 95 (2000) 235-251. Sarkar, K. and W. R. Schowalter. Deformation of a two-dimensional viscoelastic drop at non-zero Reynolds number in time-periodic extensional flows. Journal of Non-Newtonian Fluid Mechanics 95 (2000) 315-342. S.B. Pillapakkam, P. Singh, A level set method for computing solutions to viscoelastic two-phase flow, J. Comp. Phys. 174 (2001) 552-578. P. Yue, J.J Feng, C. Liu, J. Shen, Viscoelastic effects on drop deformation in steady shear, J. Fluid Mech. submitted (2004). K.S. Sheth, C. Pozrikidis, Effects of inertia on the deformation of liquid drops in simple shear flow, Computers & Fluids 24 (2) (1995) 101-119. B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging and fragmentation in multiphase flows with SURFER, J. Comp. Phys. 113 (1994) 134-147. R. Scardovelli, S. Zaleski, Direct numerical simulation of free surface and interfacial flow, Annu. Rev. Fluid Mech. 31 (1999) 567-604. D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, S. Zaleski, Volume-offluid interface tracking and smoothed surface stress methods for threedimensional flows, J. Comp. Phys. 152 (1999) 423-456. J. Li, Y.Y. Renardy, M. Renardy, A numerical study of periodic disturbances in two-layer Couette flow, Phys. Fluids, 10 (12) (1998) 3056-3071. T. Chinyoka, Numerical Simulation Of Stratified Flows And Droplet Deformation In 2D Shear Flow Of Newtonian And Viscoelastic Fluids, Ph.D. thesis, Virginia Polytechnic Institute and State University, 2004. Web address http://scholar.lib.vt.edu/theses/available/etd-11292004-130236/ J. Li, Calcul d’Interface Affine par Morceaux (Piecewise Linear Interface Calculation), C. R. Acad. Sci. Paris”, t.320 série IIb (1995) 391-396. J. Li, Résolution numérique de l’équation de Navier-Stokes avec reconnexion d’interfaces. Méthode de suivi de volume et application à l’atomisation, Université Pierre et Marie Curie PhD thesis (1996) J.E. Pilliod, E.G. Puckett, Second-order accurate volume-of-fluid algorithm for tracking material interfaces, Tech. Rep. LBNL-40744, Lawrence Berkeley National Laboratory. J.U. Brackbill, D.B. Kothe and C. Zemach, A continuum method for modeling surface tension, J. Comp. Phys. 100 (1992) 335-354. J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comp. Phys. 2 (1967) 12 - 26. I. M. Rutkevich, The propagation of small perturbations in a viscoelastic fluid, J. Appl. Math. Mech. (PMM) 34 (1970) 35-50. B. Khomami, Y. Renardy, K. C. Su, M. A. Clarke, An experimental/theoretical investigation of interfacial stability in superposed pressuredriven channel flow of well-characterized viscoelastic fluids. Part II: Nonlinear stability, J. Non-Newt. Fluid Mech. 91 (2000) 85-104. 20 [25] H. Wilson, M. Renardy, Y. Renardy, Structure of the spectrum in zero Reynolds number shear flow of the Upper Convected Maxwell and Oldroyd-B liquids, J. Non-Newtonian Fluid Mech. 80 (1999) 251-268. [26] F. Greco, Drop deformation for non-Newtonian fluids in slow flows, J. Non-Newtonian Fluid Mech. 107 (2002) 111 - 131. [27] S. Guido, M. Simeone, F. Greco, Deformation of a Newtonian drop in a viscoelastic matrix under steady shear flow, Experimental validation of slow flow theory, J. Non-Newtonian Fluid Mech. 114 (2003) 65-82. 21