Modeling the Deformation of Shear Thinning Droplets Suspended in a Newtonian Fluid


 In this work, we carried out numerical modeling of the large deformation of a shear thinning droplet suspended in a Newtonian matrix using the constrained volume model. The adopted approach was to consider making incremental corrections to the evolution of the droplet anisotropy equation in order to capture the experimental behavior of a shear thinning droplet when subjected to deformation due to imposed flow. The constrained volume model was modified by using different models to describe the viscosity of droplet phase: the Bautista et al. model, the Carreau-Yasuda model and the Power-law model. We found that by combining the constrained volume model with a simple shear thinning viscosity model we were able to describe the available experimental data for large deformation of a shear thinning droplet suspended in a Newtonian matrix. Moreover, we developed an equation approximating flow strength during droplet retraction, and we found that the model can accurately describe the experimental data of the retraction of a shear thinning droplet.


Introduction
The modeling of polymer blends composed of Newtonian fluids have been extensively investigated by researchers [1][2][3][4][5][6][7][8], and their flow behavior-to a great extent-is well understood. On the other hand, the behavior of blends composed of non-Newtonian components remains a subject that needs further development. As an example of the complexity of the topic of blends composed of non-Newtonian components, consider the work of Greco [9]. Greco [9] has derived a small deformation model based on second order fluid model, where a second order fluid is one with constant viscosity, but its stress is also affected by the presence of the first and second normal stress coefficients. As a result of the complication in the expression of the constitutive equation that describes each phase comprising the blend, the theory that was developed by Greco to describe the effect of second order fluid is mathematically tedious. And, while the behavior of a Newtonian droplet in a Newtonian matrix is controlled by two parameters, namely the capillary number (Ca) and viscosity ratio, the Non-Newtonian droplet in a non-Newtonian matrix, according to the Greco's theory, is controlled by four more additional parameters: two Deborah numbers and two ratios of the two normal stress differences for both the matrix and the droplet phases.
The previously mentioned model and other models that have tackled the subject of blends composed of viscoelastic components assume that blends are mainly composed from Boger fluids, which are fluids that have a constant viscosity and significant normal stresses [10], and those models were limited to small deformation [9,[11][12][13]. Large deformation of droplets for such case has been modeled only by a few researchers [14][15][16] and numerically simulated by [17][18][19].
For the case of one or both of the phases composed from shear thinning components, only very few experimental studies exist for the effect of shear thinning on the deformation of a single droplet. Elmendorp and Maalcke [20] investigated a variety of combinations of droplet and matrix phase fluids, including moderately shear thinning droplets in a Newtonian matrix and the opposite setup. Both systems were studied at small droplet deformation conditions. Delaby et al. [21] carried out experimental work on non-Newtonian droplets suspended in a non-Newtonian matrix. The system was subjected to a uniaxial extensional flow, and droplet lengths as they extend in extensional flow were compared to droplet lengths assuming affine deformation. It was found that when droplet deformation was predicted using an extensional rate inside the droplet that is defined using a self-consistent approach, very good agreement with experimental data was obtained. Their analytical work was carried out based on small deformation theory and it assumed a linear relationship between droplet length during extensional flow and the affine droplet length. Another study was done by Boufarguine et al. [22]. This study used a strongly shear thinning fluid as a droplet phase and subjected the shear thinning droplet to large shearing deformation. The droplet deformation was studied in addition to droplet retraction after stopping the flow. Interestingly, due to the strong shear thinning nature of the droplet fluid, the droplet, after being extended in shearing flow, stops its retraction after a short while from stopping the flow.
Some researchers have also explored the topic of a shear thinning fluid in a Newtonian matrix theoretically. Favelukis et al. [23] theoretically studied the deformation of a shear thinning single droplet suspended in a Newtonian matrix during uniaxial extensional flow. Harvie et al. [24] studied a shear thinning droplet suspended in a Newtonian liquid as it flows in a microchannel. They found that the shear thinning droplet phase can be replaced by another fluid whose viscosity is the average viscosity of the shear thinning fluid while it is within the contraction. Giraldo et al. [25] used the boundary element method to study the behavior of a shear thinning non-Newtonian droplet in a Newtonian matrix in a Couette flow. A non-Newtonian droplet deformed more than a Newtonian droplet with the same zero-shear viscosity. Peters et al. [6] have developed a constitutive equation for polymer blends based on the microstructural evolution of the anisotropy tensor, and derived a constitutive equation to describe stresses for blends with unequal viscosity between matrix and droplet phases. Later on, the Peters et al. model was compared to the rheology of blends composed of viscoelastic components by taking the viscosity to be obtained through the Carreau-Yasuda model [26]. The approach used by Peters et al. in their original work was to model the blend's microstructure by assuming passive mixing, i.e., the droplet deforms passively with the imposed flow, and hence the droplet phase deformation is independent of the viscosity ratio between the droplet and matrix phases. Therefore, in their model, the microstructure is not affected by the viscosity ratio value, although the stress is affected.
In this work, we will attempt to describe the shear thinning behavior and incorporate it in the constrained volume droplet model by implementing an approach similar to that of Tyagi et al. [26] However, the rate of deformation inside the droplet will be used to evaluate the viscosity of the droplet phase through the three different models that will be coupled with the constrained volume model. The first is the Carreau-Yasuda model, which is a steady state model that relates viscosity to applied deformation. The second is the Bautista et al. model [27], which describes shear thinning fluids through a constitutive equation. The third is a simple power law model. Then, we will compare predictions to experimental results of Boufarguine et al. [22] for shearing flow and to Delaby's results for extensional flow. To our knowledge, there are no previous studies that couple a shear thinning constitutive model to a droplet model. Although droplet elasticity can play a role in large deformation shearing flow taking place at high flow rates, this topic is outside the scope of this current work.

The Constrained Volume Model
The following equation describes the average velocity gradient tensor (κ * αβ ) inside the droplet using the theory of Wetzel and Tucker [1]: Note that we use the Einstein notation for vectors and tensors in this equation and all subsequent equations. D αβ is the applied rate of deformation tensor, Ω αβ is the applied vorticity tensor, and B αβκλ and C αβκλ are the strain rate and vorticity fourth order concentration tensors, respectively. Those fourth order concentration tensors depend only on droplet anisotropy (q αβ = 1 V d ∫︀ nα n β dS, where nα is a unit vector normal to the interface of the droplet and S is the droplet surface area) and viscosity ratio (M = η * /η) between the blend's components [28,29], where η is the viscosity of the matrix phase, η * is the viscosity of the minor phase. The procedure for obtaining B αβκλ and C αβκλ from droplet anisoropy and viscosity ratio is detailed in Appendix A. Droplet anisotropy evolves according to the assumption that changes due to its interfacial retraction and due to convective applied flow are decoupled [30].
where t is time. The droplet convection and retraction terms are described by the following two equations [31]: where Γ is the interfacial tension, Va is the approximate volume of the ellipsoidal droplet which is expressed in the model in terms of invariants of the anisotropy tensor [4], δ αβ is the identity tensor, g (︀q αβ )︀ is a function that is introduced in order to maintain constancy of droplet volume in the droplet retraction expression [4].q αβκλ is the fourth order anisotropy tensor, which is related to the second order anisotropy tensor via a closure described elsewhere [4]. f (M, ϕ) is a function of the viscosity ratio M and the volume fraction of the dispersed phase (ϕ): where c is a constant with a value of 0.695. For the evolution of the shape of a single droplet, ϕ is given a value of zero. The viscosity of the major phase is assumed to follow Newtonian behavior, whereas for the minor phase, viscosity is evaluated at the magnitude of the rate of deformation tensor inside the droplet together with the model describing the shear thinning behavior of the droplet phase, and hence η * is a value that changes during the deformation of the droplet.

Modeling Shear Thinning Behavior of Droplet Phase
The Carreau-Yasuda model is a simple model describing the steady state viscosity of a shear thinning fluid: where η * ∞ is the viscosity at infinite shear rate, η * o is the zero-shear rate viscosity and τ CY , a and n are considered to be fitting parameters. In order to have a frame-invariant equation,˙* is replaced by Constitutive models exist for describing stresses and viscosity of shear thinning fluids such as micellar solutions. One of the constitutive models that are developed for predicting stresses in micellar solutions is the Bautista et al. model [27]. Examples of other constitutive models for micellar solutions are the Johnson Segalman model [32,33], the modified Johnson Segalman model [34] and the Vasquez, McKinley and Cook model [35]. The Bautista model is capable of accurately describing shear thinning behavior of micellar solutions, in addition to the time dependent behavior of the viscosity, unlike the Carreau-Yasuda model, which provides no time-dependent behavior for the viscosity.
The Bautista et al. family of constitutive models consist of two parts: an equation describing fluidity through a Fredrickson model for fluidity, and a constitutive equation describing stress evolution. Having an equation describing stress and another describing a structure parameter (fluidity is an example of a structure parameter) is a commonly used technique in describing thixotropic fluids [36]. The original Bautista model [27] uses the Maxwell model for stress tensor evolution, although later publications used other types of constitutive models such as the Oldroyd-B model [37]. Boek et al. [38] have modified the Oldroyd-B Bautista model, since extensional viscosity predictions of the original Bautista model showed instability at large extensional flow rate values. Such model is later known as the MBM model and was used by later researchers to describe flow of micellar solutions in a circular axisymmetric capillary channel [39,40]. In this work, we will use the original Bautista et al. model since the flow type developing inside the droplet that we will investigate in this work is predominantly shearing type.
The Bautista et al. model consists of an equation for stress evolution, which is that of the upper convected Maxwell equation: coupled to a fluidity equation by Fredrickson [41], where fluidity ≡ η −1 : The last term in the above equation has been proven and validated by Stephanou et al. [42] and Stephanou [43]. σ * αβ is the stress value of the droplet phase, Go is the elastic modulus, τ B is the relaxation time of the Bautista model, and k is a structural relaxation parameter of the droplet phase. The latter parameter is taken to be a constant in this work, although other publications have assumed that k is dependent on flow strength [39]. The symbol ∇ placed on top of the stress tensor denotes the upper convected derivative. Thus, σ * αβ and η * are updated according to the value of D * αβ in the constrained volume model. It is important at this point to note that all the modeling approaches for a shear thinning fluid that we implemented in this work will alter the microstructure evolution in the CV model, because the microstructural evolution depends on viscosity ratio, and because the velocity gradient inside the droplet depends on the viscosity ratio as well. Since the viscosity of the droplet shear-thinning fluid depends on the velocity gradient inside the droplet, the picture seems to indicate a cyclical dependence between the viscosity of the droplet phase and the velocity gradient inside the droplet. The idea is illustrated in Figure 1. Such a situation is very much amenable to a numerical solution. The modified constrained volume model is composed of coupled ordinary differential equations. In our numerical implementation, we use the fourth order-Runge Kutta method with a variable step-size and relative error acceptance criterion of less than 10 −8 .

Fitting Flow Sweep Test Data of
Boufarguine et al. Figure 2 shows viscosity as a function of shear rate for a dense aqueous suspension of star-like micelles of poly(ethylene oxide) (PEO) hydrophobically end-capped with an octadecyl group that was used as droplet phase in the work of Boufarguine et al. [22]. The data are compared to predictions of the Carreau-Yasuda model [33], to predictions of the steady state Bautista et al. model and to the power law model. The Carreau-Yasuda model shows excellent agreement with the data using the parameters provided in Table 1. Note that the parameter η * o was given a very high value (4·10 8 Pa·s) because there is no evidence in the experimental data of Boufarguine et al. of a plateau viscosity value at very low shear rate values. Moreover, the η * ∞ value was set to be the viscosity of water (0.001 Pa·s). The rest of the parameters were obtained through error minimization of a standard objective function. For a power law of the type η * = A *˙B, the fit to the experimental data resulted in A = 493 Pa·s. and B = −0.806. This power law fit resulted in identical viscosity predictions as those obtained  Table 1.

Parameter
Value by the Carreau-Yasuda model, and thus identical predictions of the shear thinning droplet deformation. Finally, we note here that there is an additional prediction denoted in the legend of the plot as power law with B = −0.73. This power law viscosity prediction curve will be revisited later.
For the Bautista Model in this part, we used the steady state viscosity (η * ss ) solution of equation 7 and 8 as reported in [27] in the following form: The steady state viscosity of the Bautista et al. model was fit by specifying the value of the combined parameter kτ B . The fitting parameter value that can give reasonable results for subsequent plots was obtained by fitting the last part of the viscosity-shear rate data (i.e., data at the higher end of the shear rate values). The value of the combined parameter kτ B giving the best fit was 1.51·10 −9 s·Pa −1 , and η * ∞ is assumed to be the viscosity of water. Compared to the fit of Carreau-Yasuda, the Bautista et al. steady state viscosity predictions do not fit the experimental results as well as the Carreau-Yasuda fit does.

Comparison of CV Model Implementations to Experimental Data
We will compare in this part theoretical predictions to experimental results of a shear thinning droplet suspended in a Newtonian matrix provided by Boufarguine et al. However, since the droplet size was not explicitly stated in that paper, we use the data reported in Figures 10 and 11 of the Boufarguine et al. paper, and we assume that data points reported in those two figures have a droplet size that is the same as that used in producing other plots in the same pa-per. We found from our comparison of Figures 10 and 11 that droplet size is roughly around 150 mm for data where the viscosity of the matrix phase fluid was 30 Pa·s, and hence this value is used in all of our subsequent comparisons to the experimental data from that paper. Figure 3 shows predictions of λ Lo (droplet extension immediately after the stop of step-shear), compared to λ La , which is the normalized droplet longest dimension after a given strain assuming affine flow as given by the following expression [22]:

Comparison with Boufarguine et al. Data Under Shear Deformation
where o is the applied strain in a step-shear experiment.  Figure 2. Note that the model predictions of the droplet behavior at applied shear rate of 16 s −1 agree with the experimental data only in the beginning of the step strain, then they diverge, where the model shows behavior akin to having an effective viscosity ratio value higher than 3.5, which is the condition necessary for a droplet to start showing tumbling behavior. In fact, if we calculate the apparent droplet viscosity according to the formula: we find that the model is predicting an apparent droplet viscosity value of 296 Pa·s, giving an effective viscosity ratio value of 9.85.   However, the CV-CY agrees with the experimental data very well, when compared to the CV-P, which highly underestimates the maximum that was observed experimentally. Lastly, the CV-B model does not show any deformation, because the viscosity of the droplet phase did not drop from its high zero-shear value during the step shear time.

Comparison with Delaby et al. Data Under Aflne Extensional Deformation
A question can be raised about the accuracy of the constrained volume model itself in predicting droplet deformation as was shown in Figure 3. For that purpose, we compare the predictions of the model to the experimental  Table 2.
whereε is the applied extensional rate andε * is the extensional rate inside the droplet. This equation is coupled with Eq. 10 and equations describing the shear thinning behavior of both components to estimate the extensional rate inside the droplet and λ Lo .

Comparison with Boufarguine et al. Relaxation Data
After Step-Shear

Overall Model Assessment in Describing Blends Composed of Shear Thinning Components
The

Improvement of the CV Model for the Case of Droplet Retraction
The decoupling approximation used in the CV model is the reason why the model does not display any retraction behavior when flow is stopped for a shear thinning droplet, although experimental data shows small relaxation followed by complete freezing of the flow. Thus, an estimation of the flow strength inside the droplet during relaxation is needed in order to calculate the viscosity of the retracting shear thinning droplet. Since the model is developed in terms of anisotropy rather than droplet dimensions, an expression was developed in terms of anisotropy invariants to describe droplet longest axis (L) based on the assumption of an axisymmetric prolate droplet. This expression is the following: where Iq and IIIq are the first and third invariants of anisotropy tensor, respectively. We assume that the speed of droplet retraction is equal to ∇v . r [44], where r is a vector that represents the location of a point in space at the interface of the droplet. If we take r to be the point at the tip of the prolate droplet and that relaxation induces flow that can be approximated as the reverse of uniaxial extensional flow, then it follows that the velocity gradient tensor can be approximated as: where s, on dimensional grounds, can be estimated as s ≈ −c1 L dL dt , and c 1 is introduced to the expression as a fitting constant. The value of the parameter c 1 was obtained by fitting the CV model to the experimental relaxation data displayed in Figure 6. The experimental droplet dimensions at the start of relaxation were used to calculate the anisotropy tensor components, which were used as the initial condition for the relaxation runs of the CV model. The experimental and theoretical droplet dimensions at the end of relaxation were used to calculate modeling error. Then, the error was minimized as a function of the fitting parameter c 1 and the best c 1 value was found. Figure 7 shows the results of this implementation, using c 1 = 20.6. At the end of the step shear and the start of retraction, the model shows minor change in droplet's dimensions taking place due to droplet retraction. As the driving force for droplet retrac-

Conclusions
The constrained volume model was adapted to the case where a shear thinning droplet is suspended in a Newtonian fluid. Three approaches were suggested to describe the case. Only two approaches showed reasonable results: the constrained volume model combined with the Carreau-Yasuda model and the constrained volume model combined with the power law model. The third approach, the constrained volume model combined with Bautista et al. model, showed a failure in predicting the shear thinning viscosity behavior of the polyoxyethylene water solution of Boufarguine et al. and failed in describing droplet deformation during startup of shearing experiment. The simplicity of the power-law model makes it suitable for describing droplet deformation if it is composed of purely shear thinning fluid. On the other hand, the Carreau-Yasuda model has the advantage over the power-law model that it can describe viscosity data that has, or lacks, a plateau viscosity. The decoupling approximation, which is built in the approach of original Doi-Ohta model, implies that droplet retraction does not have any effect on internal flow. This issue was addressed by developing an expression for flow strength inside the droplet in terms of anisotropy invariants in such a way that the larger the droplet deformation, the larger the flow strength inside the droplet. By using this expression, the relaxation data of a shear thinning droplet was excellently described.

A Calculation of Average Velocity Gradient Tensor Inside the Droplet
The procedure for obtaining the average velocity gradient inside the droplet from the anisotropy tensor, the imposed velocity gradient tensor and viscosity ratio is as follows [28]: 1. Obtain the area tensor from the anisotropy tensor, where the area tensor is defined as: 2. Obtain the eigenvalues of the A αβ tensor and arrange them in decreasing order. Save the order of arrangement. 3. It is necessary to obtain the fourth order tensor S.
The S tensor has five basic components that can be related to the arranged eigenvalues of the area tensor as follows: 7. Obtain the fourth order concentration tensor Btensor in contracted notation as follows: The other components of the fourth order B tensor in contracted notation have zero values.