Microstructural evolution during complex laminar flow of liquid–liquid dispersions

J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

Microstructural evolution during complex laminar flow of liquid–liquid dispersions Eric D. Wetzel, Charles L. Tucker III∗ Department of Mechanical and Industrial Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Received 4 April 2001; received in revised form 27 August 2001

Abstract This paper presents a general framework for modeling microstructural evolution in liquid–liquid dispersions undergoing complex mixing flows. The method is implemented using the theory of Wetzel and Tucker [J. Fluid Mech. 426 (2001) 199] for dilute dispersions with zero interfacial tension. The local microstructure is represented by a second-order tensor that describes droplet shape and orientation, and the global scheme treats this tensor as a spatial field variable. A coupled finite element solution for the global flow and microstructural field is developed. Sample solutions for an eccentric cylinder mixer reveal the effect of viscosity ratio and flow history on microstructural evolution. For the calculations presented, the coupling between microstructure and rheology has only a small effect on the global velocity and microstructural fields. Some numerical difficulties associated with high strains and specific flows are discussed. © 2001 Elsevier Science B.V. All rights reserved. Keywords: Microstructured fluids; Droplet deformation; Laminar mixing

1. Introduction Many industrial processes involve microstructured fluid–fluid dispersions in complex laminar flows. By microstructured, we mean that the homogeneous regions of each phase in the dispersion are much smaller than the macroscopic flow domain. Examples include the production of polymer blends, paints, foods, and chemical emulsions. These materials undergo complex flows in twin screw extruders, batch mixers, molds, and dies, typically resulting in droplets, threads or sheets of one fluid dispersed in a continuous matrix of another fluid. This microstructure determines many important macroscopic properties, such as appearance, permeability, mechanical behavior, and rheology. To enable control of material properties, a tool is needed to predict microstructural development during complex processing operations. Modeling these systems is particularly challenging because relevant physics must be accounted for at two vastly different length scales. On the macroscale, the flow in the mixing device must be modeled. On ∗ Corresponding author. E-mail address: [email protected] (C.L. Tucker III).

0377-0257/01/$ – see front matter © 2001 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 0 2 5 7 ( 0 1 ) 0 0 1 6 1 - 6

22

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

the microscale, the deformation, reorientation, and advection of the microstructure must be calculated. While each case can be addressed individually using existing tools, little progress has been made on the combined problem. A strategy is needed which allows simultaneous representation of microstructural detail and global variations in microstructure, while accurately modeling both macroscopic complex flows and microstructural evolution. In this paper, we present a framework for modeling the microstructural evolution of fluid–fluid dispersions in complex flows. Although, the overall approach has broad applicability, this study will focus on a subclass of the general mixing problem. We will study two-phase suspensions, in which droplets of one viscosity are immersed in a matrix of another viscosity. The suspension is assumed to be dilute and have zero interfacial tension, in order to allow utilization of an analytical theory for droplet deformation under those conditions. While restrictive, these assumptions will provide a rich variety of behaviors for investigation. These physics also have many important applications, including the initial stages of blending of immiscible polymers. A number of researchers have modeled a closely related problem, the mixing of large fluid domains in complex flows. Tjahjadi and Ottino [2] used a passive interface tracking approach to model the stretching of a large fluid domain into a thin thread, and then applied a thread extension and breakup model to predict subsequent morphological development. More direct simulation approaches can also be used, including boundary element methods [3], finite element/volume-of-filling methods [4], and phase field methods [5]. These schemes can readily incorporate interfacial tension and different phase viscosities. However, they are not useful for modeling microstructured dispersions, where the overwhelming number of domains precludes the explicit simulation of every droplet. In another approach, the microscale physics are ignored and the macroscale flow simulation is used to calculate characteristic stretching parameters at each point in the flow [6–8]. These characteristic stretching values are intended to allow qualitative judgments about the microstructure, e.g. a finer microstructure should exist in areas of higher stretching. Since these models neglect microdynamics and do not explicitly calculate microstructure, at best they provide a comparative tool for analyzing stretching efficiency in mixing flows. A different strategy for modeling microstructural evolution is to define a microstructural state variable (MSV) which is evolved as a field variable in the macroscopic domain. The governing equations for the field variable are derived based on known microstructural physics. Calculations of this type have been performed in the fields of viscoelastic flows [9], forming of polycrystalline metals [10], and processing of short-fiber composites [11]. For our particular application, significant progress has already been made in deriving MSVs and their evolution equations. Doi and Ohta [12] proposed an MSV based on the weighted orientational average of interfacial area at the microstructural level, a model that was extended by Lee and Park [13] and Wetzel and Tucker [14]. Using the area tensor model of Wetzel and Tucker [14], Galaktionov [15] modeled microstructural mixing in complex flows. Their approach, however, is limited to time-periodic flows, and to dispersions with identical viscosities and zero interfacial tension. In this paper, we use the droplet model of Wetzel and Tucker [1], which allows us to model dilute systems with zero interfacial tension and dissimilar viscosities. Because the phases possess different viscosities, we are able to implement a microstructurally-dependent rheological constitutive equation to investigate coupled flow and microstructural effects. Also, since our macroscopic flow calculations are performed using the finite element method, the model is not limited to time-periodic flows.

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

23

This paper is organized as follows. In Section 2, we present the governing microstructural evolution and macroscopic flow equations. These governing equations are then cast in finite element form in Section 3, and details are given for their numerical implementation. In Section 4, microstructural calculations are performed in an eccentric cylinder mixer for various viscosity ratios. Section 5 presents strongly coupled flow calculations, in which the rheological constitutive equation is a function of the local microstructure. Section 6 provides some concluding remarks. 2. Governing equations 2.1. Microstructural evolution Wetzel and Tucker [1] derived an analytical model for the creeping flow deformation of an ellipsoidal Newtonian droplet, suspended in another Newtonian fluid with different viscosity and zero interfacial tension. The solution can be expressed in terms of the shape tensor Gij , where all points xi on the surface of the ellipsoidal droplet satisfy the relation Gij xi xj = 1.

(1)

The shape tensor is symmetric and positive definite, with six independent components. Table 1 shows some example shape tensors for specific droplet geometries. The eigenvalues of Gij are 1/ai2 , where ai is the length of one of the ellipsoidal axes. The eigenvectors of Gij indicate the orientation of the axes of the ellipsoid relative to the global coordinates. Therefore, Gij is a complete description of the size, shape, and orientation of a droplet.

Table 1 Shape tensors for special dispersion microstructures Geometry

Gij 

Spherical

0 1 0

 0 0  1

1 1  0  r2 0

0 1 0

 0  0 0

 1 1  0  r2 0

0 0 0

 0  0 0

1 1  0 r2 0 

Cylindrical

Lamellar

24

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

In the far-field, we assume that the suspending fluid flow is linear, possessing a uniform velocity gradient. This assumption implies that the droplet is much smaller than the macroscopic flow domain. We can then define an applied rate-of-deformation tensor Dij and vorticity tensor Wij as     ∂uj ∂uj 1 ∂ui 1 ∂ui Dij ≡ , Wij ≡ , (2) + − 2 ∂xj ∂xi 2 ∂xj ∂xi where ui is the fluid velocity. Bilby and Kolbuszewski [16] and Eshelby [17] showed that, for the case of an ellipsoidal droplet with zero interfacial tension in a linear far-field velocity field, the velocity field within the droplet is also linear, although in general different from that of the suspending fluid. The flow within the droplet can then be characterized completely in terms of a droplet rate-of-deformation tensor Dij∗ and droplet vorticity tensor Wij∗ . Using the results of Bilby and Kolbuszewski [16] and Eshelby [17], the droplet flow can be expressed in terms of the far-field flow as Dij∗ = Bijkl Dkl ,

Wij∗ = Wij + Cijkl Dkl ,

(3)

where Bijkl is called the strain-rate concentration tensor and Cijkl is called the vorticity concentration tensor. These dimensionless fourth-order tensors are functions of the droplet geometry and the dispersion viscosity ratio λ ≡ µ∗ /µ, where µ∗ and µ are the viscosities of the droplet and matrix, respectively. The concentration tensors can be written as Bijkl = [I − (1 − λ) S]−1 ijkl ,

Cijkl = (1 − λ)Tijmn [I − (1 − λ) S]−1 mnkl ,

(4)

where the fourth-order identity tensor Iijkl is defined as Iijkl = 21 (δik δjl + δil δjk ),

(5)

and Sijkl and Tijkl are called the Eshelby and alternate Eshelby tensors, respectively. These Eshelby tensors are only functions of the droplet geometry Gij , and are given in full form in Wetzel and Tucker [1]. Taking the material derivative of Eq. (1), the evolution of the droplet geometry is given by ˙ ij = Wik∗ Gkj − Gik Wkj∗ − Dik∗ Gkj − Gik Dkj∗ , G

(6)

˙ ij is a material derivative. For a given Dij , Wij and initial droplet shape Gij , Eqs. (3), (4) and (6) where G allow calculation of the instantaneous rate of change of the droplet geometry. The continuous evolution of the droplet shape, which Wetzel and Tucker [1] showed remains ellipsoidal for all flows, can be calculated by integrating the droplet evolution equation over time. Note that, since Gij describes the complete geometry of the droplet, both droplet shape and orientation are calculated through this implementation. The resulting droplet history is exact for arbitrary strains and arbitrary viscosity ratio, under the conditions of an initially ellipsoidal droplet with zero interfacial tension. These results were derived for a single droplet suspended in an infinite medium. However, Eq. (6) also describes the rate of change of droplet geometry for a dilute suspension of identical droplets. This interpretation is more directly applicable to the thrust of this paper. We will treat the suspension microstructure as locally uniform, but globally varying. Therefore, Gij is calculated at each material point in the global dispersion, but still represents the contributions of many droplets at that point.

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

25

Wetzel and Tucker [1] also provide examples of droplet behaviors in simple flows, and we summarize those results here. In a purely extensional flow, an initially spherical droplet stretches indefinitely. This behavior contrasts that of droplets with finite interfacial tension, which eventually breakup or reach a limiting droplet shape. Droplets with higher viscosity ratios stretch at a slower rate than those with a lower viscosity ratio. In simple shear, initially spherical droplets with λ < 3 will stretch indefinitely. Initially spherical droplets with higher viscosity ratios, λ > 5, undergo a periodic tumbling deformation. It is impossible to define a single viscosity ratio demarcating the transition from stretching to tumbling flows, since the droplet behavior will depend on both the viscosity ratio and the full three-dimensional droplet shape and orientation. However, in general, higher viscosity droplets are more likely to tumble. These behaviors in simple shear flows were verified experimentally by Comas-Cardona and Tucker [18]. 2.2. Continuity and momentum The bulk flow of the suspension is governed by the continuity and momentum equations, which for an incompressible fluid can be written as ∂ui = 0, ∂xi   ∂τij ∂ui ∂p ∂ui + uj + , =− ρ t ∂xj ∂xi ∂xj

(7) (8)

where ρ is the fluid density and p is the pressure. If the dispersion is assumed to be a Newtonian fluid of viscosity µ, then the extra stress constitutive equation is simply τij = 2µDij . More complicated rheological behavior is possible if the influence of the microstructure is taken into consideration. Wetzel and Tucker [1] showed that the rheological constitutive equation for a dilute suspension of identical ellipsoidal droplets with arbitrary viscosity ratio and zero interfacial tension is given by τij = 2µDij + 2φµ(λ − 1)Bijkl Dkl

(9)

where Bijkl is the strain-rate concentration tensor given by Eq. (4). Since Bijkl is a function of the droplet geometry, this relation shows that the suspension rheological properties are dependent on the microstructure. Furthermore, since in general the microstructure is spatially and temporally varying, the rheological behavior of the suspension will also vary with position and time. Wetzel and Tucker [1] provide calculations of rheological evolution in simple flows. In extensional flows, these suspensions typically exhibit strain-thickening behavior. In simple shear flows, shear thinning behavior is most likely for droplets in the stretching regime, while periodic thinning behavior is noted for droplets in the tumbling regime. However, since the theory applies only to dilute suspensions, the quantitative deviation of these fluids from Newtonian fluids is typically not large. Note that this rheological constitutive equation assumes zero interfacial tension, so that it will not exhibit the elastic effects observed in suspensions with significant interfacial tension [19,20]. However, rheological models for such suspensions do exist [12] and could be incorporated into this formulation. To be consistent with such a rheological model, the microstructural evolution equation would also have to explicitly include interfacial tension effects.

26

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

3. Numerical formulation 3.1. Selection of microstructural state variable Section 2.1 provides the physical framework for solving microstructural calculations. In field problems, however, the shape tensor is not an ideal MSV. The eigenvalues of the general shape tensor are guaranteed to be positive, but possess no upper bound. Therefore, large variations in shape tensor component magnitude can exist within a computational domain. The gradients associated with such large changes in magnitude are often difficult to resolve using reasonable disˆ ij ≡ Gij /Gc , cretization densities. To overcome this difficulty, we can define a normalized shape tensor, G where Gc is the trace of Gij . Unlike the general shape tensor, the eigenvalues of the normalized shape tensor lie in the range [0, 1]. This boundedness greatly reduces the magnitudes of MSV gradients within a computational domain, which we expect will improve solution accuracy and convergence. Using Eq. (6), the evolution equation for the normalized shape tensor is [21] ˆ ij ˆ ij ∂G ∂G ˆ kj − G ˆ ik Wkj∗ − Dik∗ G ˆ kj − G ˆ ik Dkj∗ + 2Dkl∗ G ˆ kl G ˆ ij . + uk = Wik∗ G ∂t ∂xk

(10)

ˆ 22 + G ˆ 33 = 1), ˆ 11 + G Note that the normalized shape tensor has only five independent components (since G and contains complete droplet shape and orientation information. This evolution equation is closed, since ˆ ij and viscosity the concentration tensors Bijkl and Cijkl are only functions of the normalized shape tensor G ratio λ. In general, the evolution equations for the normalized shape tensor would be complemented by the ˙ c = −2Dkl G ˆ kl Gc [21]. However, for our particular physical asevolution equation for Gc , which is G sumptions, droplets can deform but will not breakup or coalesce. Under these conditions, only the axis ratios and orientations of the droplets need to be evolved. This information can be completely recovˆ ij . Therefore, Gc does ered from the eigenvalues and eigenvectors of the normalized shape tensor G not provide any meaningful microstructural information for our problem, and will not be evolved as an MSV. Incorporation of other physical effects would, however, necessitate calculation of Gc . The Gc is a measure of the absolute size scale of the droplets. If interfacial tension relaxation were added to the model, without breakup or coalesce, then an initial size scale (initial Gc ) would be needed to calculate the relaxation time scale. However, further evolution of Gc would not be necessary since deformation would produce only changes in droplet axis ratios, and not droplet size. If droplet breakup and coalescence were allowed, then Gc would need to be evolved to reflect the changing sizes of the droplets.

3.2. Galerkin form of governing equations 3.2.1. Microstructural evolution The normalized shape tensor contains five independent components. We choose these independent ˆ 11 , G ˆ 22 , G ˆ 23 , G ˆ 31 , and G ˆ 12 , and arrange them in the column vector g i such that components to be G

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

  ˆ 11   G           ˆ   G   22     i ˆ g = G23 .         ˆ  G   31        ˆ   G12

27

(11)

We define {U } as the nodal solution vector containing the velocity, pressure, and shape tensor degrees of freedom. A superscript to the left of U indicates a particular type of degree of freedom, while a superscript to the right of {U } indicates a particular component of that degree of freedom, e.g. {u U i } for the ith component of velocity, {p U } for pressure, and {g U i } for the ith component of the normalized shape tensor column vector. Introducing a matrix of basis functions N for the normalized shape tensor degrees of freedom, we can write the finite element interpolation as g i = N † {g U i },

(12)

where a dagger (†) in superscript indicates the transpose of the matrix. Wetzel [21] presents a derivation of the Galerkin finite element form of the normalized shape tensor evolution equation, Eq. (10). The final result is   g i     d{ U } † † † u j j g i i NN dV Nϕ { U }(∇N ) dV { U } = Ns dV , (13) + dt V V V where V denotes the volume of the problem domain and ϕ is the matrix of velocity basis functions, so that ui = ϕ† {u U i }. (14) Here, s i represents the right-hand side of Eqn. (10) for the term corresponding to the ith component of column vector g i . For all computations in this paper, we will use a fully implicit backward Euler scheme for the time derivative of the solution variables. Using this approach, the temporal derivative of the nodal solution vector can be eliminated in favor of a finite difference approximation, and Eq. (13) can be written as         1 1 g g i † † † † u j j i i NN + Nϕ { U }(∇N ) dV { U } = NN {r−1 U } dV , Ns + (15) !t !t V V g

where !t is the time increment, and {r−1 U i } is the known solution at the previous time step. Note that the source term s i is evaluated from {g U i }, the solution at the new time. 3.2.2. Continuity and momentum Formulating the continuity and momentum equations, Eqs. (7) and (8), in a Galerkin finite element form is a routine procedure (e.g. Engelman [22], Reddy [23]). The finite element form of the continuity equation is   ∂ ϕ† ψ i dV {u U i } = 0. (16) ∂x V

28

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

The finite element formulation of the momentum equation is      u i   † U } ∂ ϕ ∂ϕ † d{ † j u i p + ρ ϕϕ dV ρ ϕu dV { U } − ψ dV { U } + (∇ ϕτ)dV i dt ∂x j V V V ∂x V  = ϕt i dS,

(17)

S

where ψ is the matrix of pressure basis functions so that p = ψ† {p U }.

(18)

The temporal derivative in the momentum equation can be rewritten following the same backward Euler finite difference scheme utilized in Eq. (15). For strongly coupled flow problems, we can rewrite the extra stress constitutive relation Eq. (9) in terms of the fourth order viscosity tensor Hijkl as τij = 2Hijkl Dkl ,

(19)

where Hijkl = µ(Iijkl + φ(1 − λ)Bijkl ).

Using this equation to eliminate the extra stress tensor unknowns from Eq. (17) gives     u i   † d{ U } ∂ϕ † † j ∂ϕ u i + ϕϕ dV ρ ϕu U } − dV { ψ dV {p U } j i dt ∂x ∂x V V V      † † ∂ ϕ ∂ ϕ ∂ ϕ ∂ ϕ jikl u l u k jikl + H dV { U } + H dV { U } = ϕt i dS. ∂x j ∂x k ∂x j ∂x l V V S

(20)

(21)

3.3. Implementation Eqs. (15)–(17) are solved using the commercial finite element software FIDAPTM . Implementation of the momentum and continuity equations are straightforward, since FIDAP allows for direct incorporation of a fourth-order viscosity tensor through user subroutines. The shape tensor components are solved individually as species transport equations of the form ∂cn ∂cn = Rn , + ui ∂t ∂xi

(22)

where cn is the species (shape tensor component) under consideration and Rn is the component rate of change notated earlier as s i . These species transport equations already exist within FIDAP, and the particular form of Rn for each component is specified using user subroutines. All calculations use a Newton–Raphson solution method. The tangent stiffness calculations are computed using analytically derived sensitivities, with the exception of some of the concentration tensor sensitivities with respect to the shape tensor primitive variables. These particular sensitivities are calculated using finite difference approximations, which greatly reduce the necessary computational time

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

29

without sacrificing significant accuracy [21]. Details of the tangent stiffness matrix calculation can be found in Wetzel [21]. Although, the hyperbolic nature of the governing equations makes this problem a strong candidate for streamline upwinding, solution accuracy was found to be better without upwinding. Therefore, no streamline upwinding is applied in any of the calculations. Linear shape functions are used for the velocity and microstructural variables, and a discontinuous pressure approximation is used. Wetzel [21] provides several test calculations that demonstrate that the code solves the governing equations correctly, and that the results converge with mesh refinement. 4. Weakly coupled calculations 4.1. Geometry and flow The 2D eccentric cylinder mixer is a standard geometry for testing microstructural evolutions in complex flows [2,24,25]. This configuration allows a complex flow, with both extensional and rotational components, to be created in a compact, closed system. The geometry is also simple to replicate experimentally and numerically, since the flow domain is easy to mesh and possesses no geometric singularities. The eccentric cylinder mixer geometry studied here is shown in Fig. 1. The geometry of the mixer is characterized by an axis ratio Ri /R0 and an eccentricity δ/R0 . The simulations presented in this section use the geometry Ri 1 = , R0 3

δ 3 . = R0 10

(23)

The flow in the mixer is determined by the relative rotations of the inner and outer cylinders. We choose to impose a stepwise mixing protocol defined as  28 πt for 0 ≤ t ≤ 0.5 5 (24) θi (t) = 0 for 0.5 ≤ t ≤ 1  0 for 0 ≤ t ≤ 0.5 θ0 (t) = . (25) πt for 0.5 ≤ t ≤ 1 − 14 5

Fig. 1. Eccentric cylinder mixer geometry.

30

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

Fig. 2. Mesh for eccentric cylinder mixer analysis. The dot indicates the node whose time evolution at large strains is plotted in Fig. 12.

In these equations t is time, θi is the angular displacement of the inner cylinder, and θ0 is the angular displacement of the outer cylinder. First, the inner cylinder rotates while the outer cylinder is stationary, and then the outer cylinder rotates while the inner cylinder is stationary. The maximum rotational displacements have no particular significance, but were chosen largely due to computational difficulties at higher strains (see Section 4.3). Wannier [26] provides an analytical velocity solution for this eccentric cylinder flow. This solution, combined with our microstructural evolution equations, could be used within a Lagrangian framework to calculate the global microstructural evolution. However, we have chosen to implement a finite element solution in order to demonstrate the generality of our approach. The finite element method can be used to compute microstructures in any flow, and allows the implementation of microstructurally-dependent rheology. The finite element mesh used for this flow is shown in Fig. 2. The elements are bi-linear 4-node quadrilaterals. The flow was calculated separately from the microstructure using a Newtonian constitutive relation, since it was found that performing a strongly coupled calculation produced very little change in the final flow and microstructural solutions. The flow is two-dimensional, so only three microstructural ˆ 11 , G ˆ 22 , and G ˆ 12 (G ˆ 33 is recovered from the trace condition G ˆ 11 + state variables need to be evolved: G ˆ ˆ G22 + G33 = 1). The initial condition for the microstructure is uniform spherical droplets. The streamlines for the stepwise protocol are shown in Fig. 3. The two streamline plots represent the two flow stages. During the first stage, when only the inner cylinder is rotating, most of the material points lie on streamlines that travel around the inner cylinder. There is one region near the outer wall on the right-hand side of the domain where smaller closed streamlines exist, and a stagnation point exists at the center of these streamlines. During the second stage, the flow is qualitatively similar, although a larger circulatory zone exists on the right side of the flow domain and the general flow direction is reversed. 4.2. Results The figures in this section show contours of the in-plane droplet axis ratio C, defined as the ratio of the shortest in-plane axis to the longest in-plane axis. The color bar for all of these contour plots is shown in

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

31

Fig. 3. Flow streamlines for the stepwise protocol as per Eqs. (24) and (25).

Fig. 4, together with representative ellipsoid shapes. Larger values of C (darker colors) represent more nearly spherical shapes. Lighter colors represent more elongated microstructures. Also shown are vector plots of the microstructural orientation. Each vector points in the direction of the longest in-plane droplet axis. The vector length is scaled by the in-plane axis ratio, so that vectors representing spherical droplets (with C = 1) have no length. Fig. 5 shows microstructural evolution for the stepwise mixing protocol and λ = 0.5. Recall from Section 2.1 that droplets with λ = 0.5 show continuous stretching behavior in shear flows. Early in the process, when the outer cylinder is stationary and the inner cylinder is moving clockwise, nearly symmetric bands of stretched microstructure (light colors) appear. At times t = 0.25 and t = 0.375, a band of more stretched material forms near the outer wall on the right half of the domain. This band develops because of the circulatory nature of the flow in that region, as seen in Fig. 3. At t = 0.5, the

Fig. 4. (a) Contour bar and (b) corresponding droplet shapes for eccentric cylinder mixer axis ratio C contour plots. This color scale is used in Figs. (5)–(11). Larger values of C (more spherical droplets) are darker; lower values of C (stretched droplets) are lighter.

32

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

Fig. 5. Contours of axis ratio C and orientation vectors versus time for the stepwise mixing protocol at λ = 0.5. (a) t = 0.125; (b) t = 0.25; (c) t = 0.372; (d) t = 0.5; (e) t = 0.625; (f) t = 0.75; (g) t = 0.875; and (h) t = 1.0.

inner cylinder stops and the outer cylinder begins to move counterclockwise. Interestingly, at t = 0.875, the flow begins to “demix”, as droplets which were previously stretched are now reduced in length. By t = 1.0 there is a significant region of lower aspect ratio droplets. From the orientation plot, it appears that the droplets in the demixed regions are being reoriented in directions not aligned with the streamlines.

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

33

Fig. 6. Effect of droplet orientation on its stretching behavior in simple shear.

Fig. 7. Contours of axis ratio C (left) and orientation vectors (right) versus time for the stepwise mixing protocol and λ = 1 at t = 1.0.

The demixing phenomenon is caused by the dependence of microstructural evolution on the orientation of the microstructure relative to the flow. The simultaneous mixing protocol produces regions of high shear throughout the mixer. A droplet will either stretch or unstretch in a shear region, depending on its orientation (Fig. 6). Wetzel and Tucker [1] show that, in simple shear, an initially spherical droplet with λ = 0.5 will undergo continuous stretching. However, if the flow changes slightly so that the droplet becomes oriented at a slightly negative angle, the droplet will begin to un-stretch. In the simultaneous mixing protocol with λ = 0.5, some of the droplets apparently become un-stretched as shearing progresses. These results show that mixing does not always improve with increased strain. Similar localized demixing has been observed for λ = 1 by Galaktionov et al. [15]. Figs. 7 and 8 show the final microstructural fields for λ = 1 and λ = 3 for the stepwise mixing protocol. The full time histories are given by Wetzel [21]. The results are qualitatively similar to the λ = 0.5 case, although, the droplet strains are generally lower as viscosity ratio increases. Fig. 9 shows microstructural evolution for the stepwise mixing protocol and λ = 5. Recall from Section 2.1 that droplets with λ = 5 show periodic behavior in shear flows. At t ≤ 0.25, the microstructural bands

Fig. 8. Contours of axis ratio C (left) and orientation vectors (right) versus time for the stepwise mixing protocol and λ = 3 at t = 1.0.

34

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

Fig. 9. Contours of axis ratio C and orientation vectors time for the stepwise mixing protocol at λ = 5. (a) t = 0.125; (b) t = 0.25; (c) t = 0.372; (d) t = 0.5; (e) t = 0.625; (f) t = 0.75; (g) t = 0.875; and (h) t = 1.0.

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

35

Fig. 10. Contours of axis ratio C (left) and orientation vectors (right) versus time for the stepwise mixing protocol and λ = 10 at t = 1.0.

develop similarly to the λ = 0.5 case, except with less overall stretching. At t = 0.375 and t = 0.5 a new phenomenon arises. Instead of axis ratio monotonically decreasing from the outer cylinder towards the inner cylinder, a minimum in axis ratio (maximum in stretch) appears at some intermediate position. A secondary band of unstretched material forms near the inner cylinder. From the orientation plots, it is clear that this effect is due to tumbling. Tracing a streamline from the immediate left of the inner cylinder to the circulatory region, the droplet orientation can be seen to be rotating with respect to the flow direction. The droplets are decreasing in axis ratio and orienting in non-streamline directions. During stage two, these demixed, unoriented droplets are advected along streamlines towards the circulatory region of the flow domain, where they further demix. At t = 1.0, there is a large demixed region in the circulatory region. There appears to be a strong streamline from the inner cylinder counterclockwise towards the demixed region, which continually supplies tumbling droplets to the demixed region. The results for λ = 10 are qualitatively similar to the λ = 5 case, although, the droplet strains are generally lower [21]. Wetzel and Tucker [1] showed that the period of oscillation for a λ = 10 droplet is roughly one-half that of a droplet with λ = 5. These bands develop more quickly for λ = 10 because the period of oscillation is shorter. Fig. 10 shows the final microstructural field for λ = 10 for the stepwise mixing protocol. We can also define a “simultaneous” mixing flow as θi (t) =

14π t 5

θ0 (t) = −

7π t 5

for t ≥ 0 for

t ≥ 0.

(26) (27)

Like the stepwise protocol, in the simultaneous flow the inner cylinder moves clockwise 14π/5 radians and the outer cylinder moves counterclockwise 7π/5 radians over 0 ≤ t ≤ 1. However, in the simultaneous flow both cylinders rotate simultaneously, while in the stepwise protocol first the inner cylinder rotates while the outer cylinder is stationary, and then the outer cylinder rotates while the inner cylinder is stationary. Fig. 11 shows the microstructural state for λ = 1 for the simultaneous flow at t = 1.0. The full time history is given by Wetzel [21]. This microstructure is different from the final microstructure for the stepwise protocol for λ = 1, even though both flows imparted identical total strains. This comparison demonstrates that dispersion microstructure is a function of both the cumulative strain magnitude and the details of the strain history.

36

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

Fig. 11. Contours of axis ratio C (left) and orientation vectors (right) versus time for the simultaneous mixing flow and λ = 1 at t = 1.0.

4.3. Numerical stability at high strains An important limitation of the current implementation was discovered by carrying out the simultaneous mixing flow, Eq. (27), to large strains. Fig. 12 shows the evolution of the normalized shape tensor components at the node labeled in Fig. 2. Note that slightly after t = 1.0, the shape tensor becomes

ˆ 11 ; (b) G ˆ 22 ; (c) G ˆ 12 ; and (d) G ˆ c ) for the Fig. 12. Time evolution of microstructural primitive variables at large strains ((a) G simultaneous mixing flow and λ = 1. Data is from the node highlighted in Fig. 2.

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

37

ˆ 11 < 0 and G ˆ 22 > 1, and the rate of growth of these non-physical components is non-physical, e.g. G very large. Stability problems were also found to occur in other geometries, including a two-dimensional strip flow, axisymmetric expansion flow, and lid-driven cavity flow [21]. In all cases, the stability problem originated in areas of high microstructural aspect ratio, where the shape tensor becomes strongly aligned with the flow. Typically, the instability is localized, occurring at 1-node or at a few neighboring nodes. Wetzel [21] provides a full analysis of the stability of highly oriented microstructures in various flows. The analysis shows that the current MSV formulation is likely to produce stability problems in steady shearing flows at large strains. As the microstructure is sheared, it becomes more elongated and more aligned in the flow direction. During this droplet evolution, the off-axis component of the normalized shape tensor in the shear plane approaches a value of zero. If round-off errors allow this component to switch signs (so that the angle Θ between the droplet’s long axis and the flow direction changes from positive to negative), then any non-physical eigenvalues will begin to grow exponentially. Since non-physical pointwise values are sometimes necessary to construct a solution which is correct in the weighted-average, finite-element sense, this instability mechanism is fairly likely in highly sheared flows. Wetzel [21] suggests some remedies to correct this shortcoming of the current formulation, including reformulation of the MSV or its evolution equation. It is also important to note that the addition of interfacial tension would provide a physical limit to the droplet aspect ratio, which may prevent instabilities of this type from occurring in most flows. 5. Strongly coupled calculation 5.1. Geometry and flow In this section, we calculate the steady microstructural field resulting from a pressure-driven strip flow, using the microstructurally-dependent rheological constitutive relation, Eqs. (20), to calculate the extra stress in the fluid. In this flow, the microstructure varies with position in the strip. Higher droplet strain occurs towards the walls and towards the strip exit, while the inlet droplets (initial condition) and centerline droplets remain spherical. The resulting gradients in microstructure will produce gradients in local suspension viscosity, which should alter the local velocity field. A strip with an aspect ratio of L/H = 20 is used for these calculations, where 2H is the total gap height. Only the upper half of the flow domain is modeled, using a 20 × 80 mesh consisting of 4-node bi-linear quadrilateral elements. The inlet velocity profile is defined as   x 2  2 u1 (x2 ) = uc 1 − 0.9025 , (28) H which provides a small amount of wall slip. Along the length of the wall, a velocity u1 /uc = 0.0975 is imposed, while a symmetry condition is imposed at x2 /H = 0. The wall slip is imposed to allow calculation of a steady-state microstructural solution. Without wall slip, microstructures at the wall would stretch indefinitely and no steady-state solution would exist. Here, u2 = 0 is imposed over the entire fluid domain. A spherical inlet microstructural condition is imposed. The steady-state solution could not be calculated directly. Instead, a transient analysis was performed, using the uncoupled steady velocity and microstructure solutions as an initial guess. The transient analysis was advanced in time until it approached

38

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

Fig. 13. Comparison of inlet and exit velocity profiles for strongly coupled pressure-driven strip flow. Note that velocity u1 /uc is plotted on the horizontal axis.

a steady result, as evidenced by the ability of the solver to achieve converged solutions for very large time steps. This final transient result was then used as an initial guess for the steady-state calculation, and the Newton–Raphson iteration converged to the steady-state velocity and microstructural fields. 5.2. Results Fig. 13 compares the inlet and exit velocity profiles for the case of λ = 3 and φ = 0.20. Although, a volume fraction of φ = 0.20 is beyond the normal bounds of applicability for a dilute rheological theory, this choice is made to exaggerate the observed rheological effects. From the figure, it is clear that the change in microstructure from the inlet to the exit causes very little change in the velocity profile. The exit velocity profile is slightly flatter than the inlet velocity profile. This effect is analogous to the velocity profile flattening observed in other shear-thinning fluids. Higher shear regions create lower effective viscosities, which increase the shear rates further in those regions. Notice that, to maintain conservation of mass, the centerline velocity at the exit is lower than the inlet centerline velocity. Similar calculations for viscosity ratios of λ = 0.5 and λ = 5.0 also showed very little change in the velocity profile. For pressure-driven strip flow and these viscosity ratios, the effect of rheological coupling seems to be small. This behavior was also found for other flows. Strongly coupled flow calculations were performed for flow in a concentric cylinder mixer, eccentric cylinder mixer, and flow past a sphere falling in a tube. In all cases, the velocity and microstructural fields for uncoupled and strongly coupled flow calculations were very similar. From these observations, we can state that rheological coupling is not important for slowly changing, shear-dominated flows at moderate viscosity ratios. However, further calculations need to be done before more general conclusions can be reached. For example, the rheological results from Wetzel and Tucker [1] suggest that for systems with viscosity ratios of λ = 0.1, significant shear thinning occurs. We were unable to reach a converged strip flow solution for this case. Ver Weyst [11] performed uncoupled and strongly coupled calculations for suspensions of short fibers, and found that rheological coupling was

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

39

most important for flows with sudden changes in geometry or boundary conditions. Again, we were unable to reach a converged solution for such flows.

6. Discussion and conclusions The calculations performed in this paper have demonstrated the feasibility of directly predicting microstructural evolution during complex flows of liquid–liquid dispersions. The predictions are based on an existing analytical model for the deformation of a droplet suspended in a second fluid with a different viscosity, for the case of zero interfacial tension. This model is implemented globally in a microstructured fluid by recasting it in finite element form, using a normalized shape tensor as the microstructural field variable. This finite element form of the governing equation is then solved along with the continuity and momentum equations, directly predicting the flow and simultaneous microstructure development of dispersions in complex geometries. This model differs from previous efforts in several significant ways. First, the calculations are based on an analytical solution for droplet evolution. This theory rigorously accounts for the effect of viscosity ratio, so that our calculations can demonstrate the effect of viscosity ratio on microstructural evolution. Finally, flow and microstructure have been directly coupled through a microstructurally-dependent rheological constitutive equation. The eccentric cylinder mixer examples show the importance of processing geometry, mixing protocols, and dispersion properties on global microstructural development. Microstructures generally align along streamlines, stretch more in areas of high strain rates, and deform less with higher viscosity ratios. Flows with non-uniform deformation rate develop highly non-uniform microstructures. The final global microstructure is not only dependent on the total strain imposed on the dispersion, but also on the history of its application. In fact, the ‘demixing’ phenomena demonstrates that imposing more strain does not always lead to better mixing. Strongly coupled flow and microstructural field calculations have been performed by implementing a rheological constitutive equation with explicit microstructural dependence. For the particular flows and dispersion properties under consideration, the microstructure had little effect on the global flow behavior. However, microstructural effects may be more significant in flows with sharper changes in streamline direction. The calculations also revealed an inherent instability which exists in the current numerical formulation. The instability is caused by a combination of non-physical microstructural field variable values, and the intolerance of the evolution equations to these non-physical states in certain flows. Wetzel [21] has suggested some modifications to the state variables and evolution equations in order to improve the stability of this approach, but more work is needed. Although, this paper only presents two-dimensional calculations, the analytical microstructural model which is the basis of our calculations is fully three-dimensional. Therefore, our approach can be readily extended to three-dimensional flows. While the current calculations apply only to dilute droplet dispersions with zero interfacial tension, the general framework presented here is easily extended to other mixing problems by using appropriate microstructural models. For dilute Newtonian dispersions, interfacial tension can be added to the present model [27], and droplet breakup can be modeled after the droplets have formed into threads [2]. The Doi and Ohta [12] theory and its variants offer one way to treat more concentrated mixtures, which may form

40

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

co-continuous microstructures. At the present time, most of these models are limited to Newtonian phases. The development of microstructural evolution equations for mixtures of viscoelastic phases remains as an important challenge for future work. Acknowledgements This work was supported in part by the National Science Foundation, Grant No. DMI-981320. Eric Wetzel was supported as a National Science Foundation Graduate Research Fellow. References [1] E.D. Wetzel, C.L. Tucker, Droplet deformation in dispersions with unequal viscosities and zero interfacial tension, J. Fluid Mech. 426 (2001) 199–228. [2] M. Tjahjadi, J.M. Ottino, Stretching and breakup of droplets in chaotic flows, J. Fluid Mech. 232 (1991) 191–219. [3] L. Stradins, T.A. Osswald, Predicting the effect of viscosity ratios on the mixing of polymer blends using the boundary element method, Polym. Eng. Sci. 36 (1996) 979–984. [4] D.F. Zhang, D.A. Zumbrunnen, Influences of fluidic interfaces on the formation of fine scale structures by chaotic mixing, J. Fluids Eng. 118 (1996) 40–47. [5] R. Chella, J.Viñals, Mixing of a two-phase fluid by cavity flow, Phys. Rev. E 53 (1996) 3832–3840. [6] S.J. Kim, T.H. Kwon, Measures of mixing for extrusion by averaging concepts, Polym. Eng. Sci. 36 (1996) 1466–1476. [7] C. Yao, I. Manas-Zloczower, Study of mixing efficiency in roll-mills, Polym. Eng. Sci. 36 (1996) 305–310. [8] D.I. Bigio, J.H. Conner, Principal directions as a basis for the evaluation of mixing, Polym. Eng. Sci. 35 (1995) 1527–1534. [9] F.P.T. Baaijens, Mixed finite element methods for viscoelastic flow analysis: a review, J. Non-Newtonian Fluid Mech. 79 (1998) 361–385. [10] U.F. Kocks, C.N. Tome, H.R. Wenk, Texture and anisotropy: Preferred Orientations in Polycrystals and Their Effect on Materials Properties, Cambridge University Press, Cambridge, 1998. [11] B.E. VerWeyst, Numerical Predictions of Flow-Induced Fiber Orientation in Three-Dimensional Geometries, Ph.D. thesis, University of Illinois at Urbana-Champaign, 1998. [12] M. Doi, T. Ohta, Dynamics and rheology of complex interfaces, Part I, J. Chem. Phys. 95 (1991) 1242–1248. [13] H.M. Lee, O.O. Park, Rheology and dynamics of immiscible polymer blends, J. Rheol. 38 (1994) 1405–1425. [14] E.D. Wetzel, C.L. Tucker, Area tensors for modeling microstructure during laminar liquid–liquid mixing, Int. J. Mult. Flow 25 (1999) 35–61. [15] O.S. Galaktionov, P.D. Anderson, G.W.M. Peters, C.L. Tucker, A global, multi-scale simulation of laminar fluid mixing: the extended mapping method, Int. J. Mult. Flow, in press (2002). [16] B.A. Bilby, M.L. Kolbuszewski, The finite deformation of an inhomogeneity in two-dimensional slow viscous incompressible flow, Proc. Roy. Soc. A355 (1977) 335–353. [17] J.D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. Roy. Soc. A241 (1957) 376–396. [18] S. Comas-Cardona, C.L. Tucker, Measurements of droplet deformation in simple shear flow with zero interfacial tension, J. Rheol. 45 (2001) 259–273. [19] C. Lacroix, M. Aressy, P.J. Carreau, Linear viscoelastic behavior of molten polymer blends: a comparative study of the Palierne and Lee and Park models, Rheol. Acta 36 (1997) 416–428. [20] G.K. Guenther, D.G. Baird, An evaluation of the Doi–Ohta theory for an immiscible polymer blend, J. Rheol. 40 (1996) 1–20. [21] E.D. Wetzel, Modeling Flow-Induced Microstructure of Inhomogeneous liquid–liquid Mixtures, Ph.D. thesis, University of Illinois, Urbana-Champaign, 1999. [22] M. Engelman, Fidap 7.0 Theory Manual, Fluid Dynamics International, 1st Edition, 1993. [23] J.N. Reddy, An Introduction to the Finite Element Method, 2nd Edition, MacGraw Hill, UK, 1993.

E.D. Wetzel, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 101 (2001) 21–41

41

[24] J.M. Ottino, The Kinematics of Mixing: Stretching, Chaos, and Transport, Cambridge University Press, Cambridge, 1989. [25] A.N. Beris, R.C. Armstrong, R.A. Brown, Finite element calculation of viscoelastic flow in a journal bearing. II: moderate eccentricity, J. Non-Newtonian Fluid Mech. 19 (1986) 323–347. [26] G.H. Wannier, A contribution to the hydrodynamics of lubrication, Q. Appl. Math. 8 (1950) 1–32. [27] N.E. Jackson, C.L. Tucker, A compact model for droplet deformation with interfacial tension, in: Proceedings of the PPS-17 (CD-ROM), Montreal, May 2001.