A numerical model for the development of the morphology of disperse blends in complex flow

The aim of this study is to develop a constitutive model for disperse blends applicable in complex flows and to cast this model in a finite element framework. As the number of droplets in realistic conditions is extremely large, it is computationally intractable to model all droplets individually. Droplet populations are modeled that have macroscopically averaged morphological properties. These properties are the droplet stretch ratio, the unstretched droplet radius, the orientation vector, and the number of droplets per unit volume. The evolution equations of these properties vary based on the morphological state transitions. The current model describes the morphology evolution in complex geometries, assuming Newtonian mixture constituents and monodisperse droplet populations. The numerical model has been validated for simple shear flow. Results are discussed for Poiseuille flow and the eccentric cylinder flow.


Introduction
An important class of materials is polymer blends. The majority of polymer blends are immiscible, because mixing long polymer chains is thermodynamically unfavorable. Consequently, many polymer blends have a multiphase structure. There is a huge variety of structures that can be produced, such as dispersed droplets (dispersive blends), cocontinuous blends, and fibril-like structures. In this work, we focus on modeling dispersive blends. Dispersive mixing is an often used method in industry, not only in polymer blending, but also in, for example, food processing and the pharmaceutical industry. To obtain the desired material properties, it is crucial to be able to control the droplet size of the disperse phase. Depending on the emulsion properties and the processing history, the droplets can undergo changes as a result of deformation, orientation, breakup, and coalescence.
Many authors have contributed to the modeling and experimental work of the deformation and breakup of single Patrick D. Anderson p.d.anderson@tue.nl 1 Polymer Technology, Department of Mechanical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands droplets (see for example Minale (2010)). The case of small deformation and Newtonian liquids has been investigated by Taylor (1932) and Taylor (1934). Tomotika (1935) analytically investigated the breakup of an infinitely long viscous filament embedded within a matrix of another viscous liquid. This work was expanded upon by Mikami et al. (1975) and Khakhar and Ottino (1987). Cox (1969) described the deformation of an ellipsoidal droplet in a general time-dependent flow field. Doi and Ohta (1991) follow the approach of the influential paper of Batchelor (1970) by describing the morphology in terms of an interface tensor, quantifying the expansion, and orientation of complex interfaces under the influence the macroscopic velocity field and the interfacial tension. Another model describing the morphology with a shape tensor has been developed by Maffettone and Minale (1998). Ellipsoidal drop shape predictions have been verified by threedimensional visualization by Guido and Greco (2001). Both the Doi-Ohta model and the Maffetone and Minale model have been generalized using the GENERIC formalism from Grmela andÖttinger (1997), applied on immiscible blends by Grmela et al. (2001). Using the GENERIC formalism, Edwards et al. (2003) developed a range of rheological models that obey certain physical constraints, such as volume conservation of the disperse phase in dispersive polymer blends. Influential experimental work on the breakup conditions has been carried out by Grace (1982). Further experimental investigations were carried out by Bentley and Leal (1986a, b), Stone et al. (1986), Grizzuti and Bifulco (1997), and Vinckier et al. (1997). Excellent reviews on these topics have been given by Stone (1994) and Tucker and Moldenaers (2002). Development of a model to describe the morphology of immiscible Newtonian blends undergoing morphological changes from the initial structure to the final structure after deformation, breakup and coalescence has been done by Peters et al. (2001) and Almusallam et al. (2004).
The objective of this study is to develop a phenomenological model for predicting the morphology evolution of disperse blends in complex flows. Our starting point is the constitutive dispersive blend model of Peters et al. (2001). This model has been shown to be very successful for describing dispersive blend morphology in shear flow. However, this model has discontinuous solutions in time and space, hence is not attractive for implementation in a finite element formulation. The aim of this work is to extend this model, introduce smoothed field variables and apply our methodology for complex flow problems.

Droplet morphology
During processing, the droplet structure can change. This process is determined by the interplay between the viscous stress exerted by the surrounding matrix fluid which deforms the interface of the droplet, and the interfacial tension which tends to restore the droplet to a spherical shape. This process is governed by the dimensionless capillary number Ca, which is the ratio between the viscous forces and the interfacial tension forces. For shear flow, this is given by: where μ m is the dynamic viscosity of the matrix phase,γ is the shear rate, R 0 is the radius of the unstretched droplet and σ is the interfacial tension between the two fluid phases. Another important dimensionless quantity for describing this process is the viscosity ratio λ, as given by: where μ d is the dynamic viscosity of the droplet phase. Hydrodynamic forces exerted by the matrix will tend to cause deformation and orientation of the droplet. Depending on the capillary number and viscosity ratio, the droplet can become unstable as a result of the deformation. There is a critical capillary number Ca crit above which the droplet becomes unstable. This critical capillary number depends on the viscosity ratio and the flow type. Grace (1982) has measured this relation for shear flow and planar extensional flow. For simple shear flow, this is given by De Bruijn (1989): and for planar extensional flow, the critical capillary number is given by: Jansen et al. (2001) did measurements on the critical capillary number in concentrated emulsions for a wide range of viscosity ratios and volume fractions and observed that Ca crit decreases with increasing volume fractions. However, it was also observed that the Grace curve, which was determined for breakup of single droplets, can be used if the matrix viscosity μ m is replaced by an effective matrix viscosity. This makes sense, because in a (concentrated) emulsion, a droplet is not surrounded by pure matrix fluid, but by matrix fluid containing a large number of other droplets. The effective matrix viscosity is given by Choi and Schowalter (1975): where φ is the volume fraction of the droplet phase in the total fluid mixture. Droplets break up via mainly four mechanisms (Jansen et al. 2001;Van Puyvelde and Moldenaers 2005). The first one is necking (also called waist-thinning or binary breakup). In this case, a slightly deformed droplet simply splits up into two daughter droplets. The second breakup mechanism is filament breakup (also called capillary breakup), where a droplet has been stretched into a long thin filament and breaks up into a line of multiple daughter droplets under the influence of the Rayleigh-Plateau instability. The third breakup mechanism is end-pinching, where the ends of a filament break off as new daughter droplets, which can happen due to collision with neighboring droplets or when the flow is stopped at intermediate aspect ratios of the deformed droplet. The fourth breakup mechanism is tipstreaming, where the ends of a deformed droplet become cusp-like and small droplets are ejected from the tips. This occurs due to Marangoni stresses that arise as a result of a non-uniform surfactant distribution. In this work, only the first two mechanisms are considered (necking and filament breakup). Marangoni stresses do not occur, because we study only single populations of clean drops, which do not exhibit end-pinching and tip-streaming (Bazhlekov et al. 2006;Janssen and Anderson 2008).
Aside from deformation, orientation, and breakup, the last important phenomenon for describing droplet structures is coalescence. In coalescence, droplets collide with each other and merge into a single larger droplet. This is important because the aim of mixing is often to obtain droplet sizes that are as small as possible, whereas coalescence coarsens the morphology. The success of a coalescence event depends on two factors, namely, the frequency with which droplets collide and the probability that the film of matrix fluid between the two droplets is completely drained within the available process time. The drainage can be modeled assuming either fully immobile interfaces, partially mobile interfaces and fully mobile interfaces. Coalescence will be successful when the film thickness is reduced to a critical value h crit , as given by Chesters (1991): where A h is the Hamaker constant and R eq is the equivalent radius of the two differently sized droplets. In this work, it is assumed that R eq = R 0 . The analysis on the coalescence probability by Janssen (1993) shows that coalescence occurs mostly in regions with low Ca, because that is needed to allow enough interaction time for the film to drain completely, though at too low Ca, the probability of collisions occurring is very low, then nothing happens with the droplets.

Morphological quantities and state changes
Several approaches can be followed to describe the droplet structure in blends, such as sharp interface models, diffuse interface models and macroscopic models. In the first two types of models, the droplet interfaces are tracked directly, but this quickly becomes computationally intractable as the number of droplets becomes as large as what is seen in industrial blends. In the case of macroscopic models, a microstructural model of the droplet morphology development is used and translated to a number of macroscopically averaged quantities. We describe the blend structure using a macroscopic model and use the model of Peters et al. (2001) as a starting point. This model describes the blend structure phenomenologically and couples this to a stress contribution. We focus on the part describing the blend structure. It has been demonstrated to be very successful for describing the droplet structure in shear flow. However, the solutions are discontinuous in time and space, so it cannot easily be implemented in a finite element formulation in this form. Therefore, we introduce smoothed field variables that can be solved in a complex flow geometry. In every point in space x, there is a population of droplets with macroscopically averaged quantities. These quantities are the stretch ratio β(x), unstretched droplet radius R 0 (x), orientation vector m(x) and number of droplets per unit volume N d (x). The stretch ratio is defined as the length of the major axis of an extended ellipsoid droplet divided by the diameter of the unstretched droplet. The orientation vector is described using the b tensor, which is defined as the contravariant decomposition of the c tensor, with c = b · b T (Hütter et al. 2018). The tensor c is governed by a Giesekus model. The orientation vector m is defined as the eigenvector corresponding to the largest eigenvalue of the c tensor. Peters et al. (2001) use the Finger tensor B instead of c and the deformation gradient tensor F instead of b. Our definition for the orientation vector m is similar to how Peters et al. (2001) define it as the eigenvector corresponding to the largest eigenvalue of the Finger tensor B, which gives the direction in which the background flow field has the largest deformation.
The blend structure evolves under the influence of the flow field. The interaction is assumed to be one-way only, which means the flow field influences the morphology but not vice versa. Stretched droplets are modeled as cylinders, with a radius R that is determined by conservation of volume: The morphology is assumed to change only via four mechanisms, namely stretching into filaments, filament breakup (capillary breakup), necking (binary breakup) and coalescence (see Fig. 1). Which mechanism occurs depends on Ca, using the effective matrix viscosity μ e , the macroscopically averaged unstretched droplet radius R 0 and the effective shear rate, defined by: where D = 1 2 (∇u + (∇u) T ) is the symmetric part of the velocity gradient tensor L = (∇u) T , with u the velocity field. When Ca < Ca crit , droplet breakup cannot occur. In this regime, the morphology is described by coalescence. The rate of coalescence depends on the volume fraction of the droplet phase and the shear rate. When Ca is very large (Ca ≥ κ Ca crit ), the droplet first deforms into a long thin filament. The parameter κ depends on the flow type and has the value 2 for shear flow and 5 for planar elongational flow. For complex flows, the value for κ will be between 2 and 5. For simplicity, we assume κ = 2 for all cases in this work. This is valid for the cases of simple shear and Poiseuille flow. However, this assumption is not valid for complex flow problems with a significant extensional contribution, such as an eccentric cylinder flow problem. This topic will be investigated in a subsequent paper. When the filament radius R reaches a critical filament radius (R ≤ R crit ), Fig. 1 Illustration of the four considered morphological changes. From top to bottom: coalescence, necking, filament stretching, and filament breakup the filament breaks up via capillary breakup. The critical filament radius is given by Tjahjadi and Ottino (1991): where α 0 is the initial disturbance amplitude (which eventually grows sufficiently large to make the filament unstable), λ e = μ d /μ e is the effective viscosity ratio and m is the orientation vector, which is calculated as the eigenvector corresponding to the largest eigenvalue of c. α 0 is of the order O(10 −9 ) (Janssen 1993). Lastly, the intermediate Ca regime, where Ca crit ≤Ca < κ Ca crit , is described by necking.

Morphology evolution equations
The independent macroscopic variables describing the morphology are the stretch ratio β, the unstretched droplet radius R 0 , and the b tensor (from which the orientation vector is derived). The number of droplets N d can be derived from the droplet radius via conservation of volume. The evolution equations are given by: where the right-hand side of the b tensor is given by Hütter et al. (2018): This is the Giesekus model written in terms of b instead of the conformation tensor. We define the b tensor in this way to prevent it from growing indefinitely. With this method, b will eventually relax towards a certain steady value, depending on the parameters τ G and α G . For τ G →∞ in Eq. 13, b goes to the deformation gradient tensor F and c goes to the Finger tensor B. The functions f 1 and f 2 vary based on the mechanism of morphological change.
In the case of coalescence, the right-hand side functions are given by: with φ the volume fraction of the disperse phase in the blend, τ c a numerical smoothing parameter and .
The expression for f 2 is derived by Peters et al. (2001), using the calculation of the droplet radius during coalescence from Janssen (1993). The value of β + is the value that the stretch ratio β would discontinuously jump to according to the Peters model. In that model, the stretch ratio is not directly used in the stress calculation, but an intermediate step of calculating an effective stretch ratio β eff via an ODE of the form of Eq. 14 is performed. This step smoothes the discontinuous jump in time. We follow the same approach but apply the smoothing step directly to the stretch ratio field variable β itself. This smoothing in time effectively also smoothes discontinuous jumps along the streamline direction. However, there is no smoothing of discontinuous jumps in crosswind direction, as is seen in the results obtained for Poiseuille flow in "Poiseuille flow". In the case of necking, the right-hand side functions are given by: with τ n a numerical smoothing parameter and The expression for f 2 is taken from the Peters model and the smoothing step for the stretch ratio with f 1 as given by the right-hand side term of Eq. 10 is similar to the case of coalescence as given by Eq. 14.
In the case of filament stretching, the right-hand side functions are given by: with The quantity Ca ext is the extensional capillary number, as it is calculated based on D: mm (which is a measure of extension) instead of the equivalent shear rateγ . Since a stretching droplet does not change its undeformed radius R 0 , we have set f 2 = 0. The expression for f 1 is derived by Stegeman et al. (1999). In the case of shear flow, the droplet will initially be oriented in the direction of highest shear (D: mm is maximal). As the shear is maintained, the droplet rotates away from this direction and D: mm decreases to a very small value, depending on α G . In the case of filament breakup, the right-hand side functions are given by: with τ b a numerical smoothing parameter and In this case, both β and R 0 are smoothed as explained for Eq. 14. Equation 26 is important for the selection of the smoothing parameter. When smoothing is performed too slowly (i.e., with a too large value for the smoothing parameter τ b ), the capillary number evolves significantly differently, which means the droplet population will experience the morphological changes of filament stretching, filament breakup, necking, and coalescence very differently. Therefore, the smoothing parameter is chosen not too large.
In this work, we set it to be at least one order of magnitude larger than the numerical time step but not much higher. Earlier, it was mentioned that D: mm can go toward a very small value. If this happens, then R crit and R + 0 will go to infinity (see Eqs. 9 and 28). However, even if R crit is large but finite, there is still the issue that, because the droplet size after filament breakup is calculated according to Eq. 9, this can lead to the situation that the daughter droplets after the filament breakup are actually larger than the parent droplet before the breakup event. An example is demonstrated in "Influence of shear rateγ ". The reason is probably that Eq. 9 was derived for a filament in a stretching flow (with nearly constant extensional rate), which means it is not valid for too small extensional rates D: mm, where the flow can be considered to be practically stagnant. This requires separate modeling of filament breakup under quiescent conditions, which we have not done in this work. This topic will be investigated in a subsequent paper.

Finite element formulation
First, we introduce the logarithmic quantities s = log(β) and v = log(R 0 ). The reason is that s and v have the useful property that they are allowed to become negative, whereas β and R 0 could become negative due to numerical errors when simulated directly, which would lead to meaningless results. For example, Eq. 7 would take the square root of a negative number, which stops the numerical simulations. Eqs. 10 and 11 are replaced by: respectively. In other words, the stretch ratio and unstretched droplet radius are described according to: where β = exp(s) and R 0 = exp(v). Before these can be implemented in a finite element framework, these equations are rewritten in a weak formulation, by multiplying the equations with test functions w 1 and w 2 and integrating over the whole domain , to obtain: Here (.,.) denotes the inner product on . Similarly, the weak formulation of the tensor b becomes: The tensor b is solved component-wise. In this work, w 1 , w 2 and the components of W are chosen from the same function space. Flow fields are generated by solving the velocity u and the pressure p using the incompressible Stokes equations without body forces on different geometries. To elaborate, the flow field is determined from the continuity equation and momentum balance. Currently, the morphology is not coupled back to the velocity. The velocity field is calculated once, after which the morphology evolution is determined using this known velocity field. This backcoupling may be an interesting topic for future work.

Numerical implementation
The time-integration method is first-order semi-implicit, which means the material derivative is solved implicitly and the right-hand side function is evaluated explicitly. This leads to the following temporal discretization: where t is the time step. Putting the unknowns on the left hand side and the known quantities on the right-hand side results in: The convection terms are stabilized by using SUPG (streamline upwind Petrov-Galerkin) test functions, which means an extra stabilization term is added in the streamline direction. It amounts to replacing the test function w by w + τ s u ·∇w, where τ s is the SUPG parameter. This parameter is given by: where h is a characteristic element size (Hughes et al. 1986) and U is a characteristic velocity. There are several choices possible for U. In this work, U is calculated using the magnitude of the velocity u in every integration point. We found this to be the most stable choice. First the flow field is generated by solving the momentum and continuity equations, then the morphology equations are solved according to the following scheme. At each time step, the problem of b is solved first, followed by the problem of s (or β) and finally the problem of v (or R 0 ). Afterwards, during the next time step, the code loops back to b (see Fig. 2).
The velocity u is discretized using second-order polynomial basis functions, while p,s,v, and b are discretized using first-order polynomial basis functions.

Validation with simple shear flow
This section shows results for simple shear flow simulations with varying values of the smoothing parameter τ and the shear rateγ . We set the numerical smoothing parameters for filament breakup (τ b ), necking (τ n ) and coalescence (τ c ) equal, and define this value as τ . The main benefit of simple shear flow is that the shear rate is identical across the whole domain, which means the convection term in the morphology evolution equations drops out. This makes it possible to study the morphology evolution purely under the influence of deformation. Furthermore, the morphology can be simulated without using numerical smoothing, so results can be compared with the unsmoothed results by Peters et al. (2001).

Influence of numerical smoothing parameter τ
The first result (see Fig. 3) is a reproduction of the result already obtained by Peters et al. (2001). In this simulation, the time step is taken as t = 10 − 4 s and no numerical smoothing is applied. The material properties are as given in Table 1. The material parameters are taken as those of a 10% PIB-90% PDMS blend (same as Peters et al. (2001)). These values are used for all the simulation results in this work. First, filament stretching occurs, as can be seen from Fig. 3 The stretch ratio β over time t, without numerical smoothing. The shear rate isγ = 3s −1 and the time step is t = 10 − 4 s the gradual increase of the stretch ratio. Then suddenly filament breakup occurs, where the extremely extended filament breaks up into many smaller droplets. After some further breakup due to necking, a dynamic equilibrium between necking and coalescence is reached, which can be seen by the occasional upward spike indicating that the capillary number is momentarily larger than the critical capillary number. The sharp transitions of the stretch ratio are instantaneous (the up and down parts each span only one time step), as no numerical smoothing has yet been applied.
However, these instantaneous jumps are not physical, so we decide to numerically smooth these discontinuous processes. With a smoothing time of τ = 0.1 s, the stretch ratio evolution has also been simulated (see Fig. 4). It is seen that the extreme jumps have indeed become smoother. The upward spikes in the dynamic equilibrium regime are no longer visible after the smoothing, as they only encompassed a relatively small part of one "period". Note that the number of spikes in this final regime is not fixed (in the case without numerical smoothing), but dependent on the time step, because the time step determines how far the single necking step puts the capillary number Ca below the critical capillary number Ca crit . With a smaller time step, Ca  Fig. 4 The stretch ratio β over time t, with numerical smoothing τ = 0.1 s. The shear rate isγ = 3s −1 decreases less below Ca crit after the instantaneous jump, so coalescence will make Ca reach Ca crit sooner, so upward spikes due to necking will occur more frequently in that case.

Influence of shear rateγ
A simulation has been performed for a shear rate ofγ = 9s −1 (see Fig. 5 for the droplet size as function of time). There was no numerical smoothing. It is seen that the droplet size evolves not as expected. Instead of the expected sequence of filament stretching, then filament breakup followed by necking until the dynamic equilibrium state, the system switches between filament stretching and filament breakup while making the droplet larger with every filament breakup step. This is because the droplet size after filament breakup is always slightly larger than the critical filament radius R crit , as given by Eq. 28, so R 0 keeps increasing because R crit keeps increasing due to orientation of the orientation vector m with the flow field. The orientation vector is directly determined from the c tensor, which is not representative of the newly formed droplet after the breakup event.
To prevent b from growing indefinitely, the Giesekus relaxation term is added to Eq. 13. However, applying this relaxation of b still does not eliminate the problem of the unexpected morphology evolution sequence, as the relaxation time τ G is usually chosen large compared to the timescale of the physical process. The issue seems to be that b is not linked to R crit for low stretching rates. Equation 9 is not valid for too small values of D : mm, which occurs in simple shear flow when the droplet orients away from the direction of the background flow with the largest stretch (at 45 • ).
A separate morphological transition for filament breakup under quiescent conditions should be taken into account to model this situation, which we have not done in this work. We instead introduce another method to handle this situation. Since b and R crit affect only the part of filament stretching and filament breakup, it is decided to put an upper limit R crit,max on R crit , in such a way that the capillary number after breakup obeys Ca < κCa crit , so that the next morphological state will always be necking (see Eq. 43). R crit,max is given by: Including this upper limit, the expected sequence is obtained (see Fig. 6). Using Eq. 43, R crit in Eq. 9 will not go to infinite in the instance that D: mm goes to zero. The effect of the shear rate on the morphology development is investigated (see Figs. 7 and 8). As expected, the higher shear rate ofγ = 9s −1 leads to faster filament stretching than the shear rate ofγ = 3s −1 . In the densely colored region, the (unsmoothed) β rapidly alternates between the value for necking and coalescence. Furthermore, filament breakup is delayed, as the filament requires more deformation to become unstable. However, no filament stretching occurs at the shear rate ofγ = 1s −1 . This simulation immediately goes into the necking regime and continues until the dynamic equilibrium between necking and coalescence is reached. All in all, it is observed that as expected, the final droplet size decreases with increasing shear rate. We have validated our phenomenological model for simple shear flow and in the subsequent sections, we show results of our blend model in Fig. 6 The unstretched droplet radius R 0 , the filament radius R and the critical filament radius R crit as a function of time for a shear rate ofγ = 9s −1 . The solid black line indicates Ca = κCa crit and the dashed black line indicates Ca = Ca crit . No numerical smoothing has been applied. An upper limit has been imposed on R crit so that necking occurs after filament breakup Poiseuille flow and an eccentric cylinder flow, both of which are relevant for practical applications.

Poiseuille flow
In this section, we show results for Poiseuille flow, which essentially consists of many layers of simple shear flow at different shear rates. Poiseuille flow is a classical example Fig. 7 The stretch ratio over time for three different shear rates, without numerical smoothing Fig. 8 The unstretched radius over time for three different shear rates, without numerical smoothing of a pressure-driven flow, where a layered droplet structure arises near the walls, which is relevant for many industrial applications. The domain is a rectangle, with a mesh consisting of two elements in the (horizontal) x-direction and N y elements in the (vertical) y-direction. The time step is t = 10 − 3 s and the smoothing parameter is τ = 10 − 1 s. The top and bottom boundaries are no-slip boundaries. The left and right boundaries are coupled through a periodic boundary condition and a flow rate of Q = 10 m 2 /s is imposed. The length in y-direction is L y = 4 m and the length in x-direction is set equal to the width of two elements. The elements are squares, with biquadratic interpolation for velocity, bilinear interpolation for pressure and all the morphology variables. The initial droplet size is R 0 = 9.5 · 10 − 6 m and the initial stretch ratio is β = 1.

Comparison between standard and logarithmic versions of the blend model
As was discussed before, our model does not solve for the stretch ratio β and unstretched radius R 0 directly, but their logarithmic counterparts s and v. This section demonstrates why the latter is the better choice.
In Poiseuille flow, the shear rate is highest near the noslip walls. Here, the droplets first undergo stretching into filaments, until the filament radius R decreases below the critical filament radius R crit , after which filament breakup occurs. When this occurs, the stretching filaments in the adjacent layer are not yet breaking up. During filament breakup, the blend model describes a nearly discontinuous jump in both β and R 0 . The numerical smoothing technique only smoothes discontinuities in the streamline direction, not perpendicular to it. This leads to the occurrence of the Gibbs phenomenon (see Fig. 9), because the continuous Fig. 9 The stretch ratio in Poiseuille flow simulated using the standard formulation (red line) and derived from the logarithmic formulation (blue line) at time = 11.1 s. There are 50 elements in y-direction finite element basis functions cannot capture such a discontinuity. In the region of 0.5 < y < 1 m, it can be seen that there is an almost discontinuous jump between a region where filament breakup has already occurred and a region where filaments are still stretching. Of particular interest is the downward peak of the red line. When the mesh is too coarse or the discontinuity is too large, this peak can reach a negative value. This is undesirable for two reasons.
First of all, a negative stretch ratio has no physical meaning and secondly, inserting a negative β into Eq. 7 leads to taking the square root of a negative number, which leads to a floating point exception error, making the numerical simulation stop. The logarithmic formulation with s and v has the benefit that β = exp(s) and R 0 = exp(v) cannot become negative.
Another issue with the standard version of the blend model is that it seems to display the Gibbs phenomenon more strongly than the logarithmic version (see Fig. 10). The evolution of the local morphology depends directly on the latest state, so these extreme spikes can cause the morphology to evolve differently than expected (see Fig. 11). It is observed that the Gibbs phenomenon is much weaker when the mesh is refined four times (see Fig. 12). In this case, the results for the standard and logarithmic versions are observed to be very similar. In the outer regions, there is dynamic equilibrium between necking and coalescence, and in the center there is a region where coalescence dominates. The logarithmic formulation was found to be particularly useful for improving the numerical stability of the morphology solution, for relatively coarse meshes.

Mesh convergence
In this section, it is investigated how the solution at a very long time scale (2 · 10 6 time steps, t = 2 · 10 3 s) relates to the reference solution that is obtained by combining simple shear flow solutions at the same shear rate values as encountered in the Poiseuille flow problem. This comparison is useful, because the Poiseuille flow is affected by the oscillations of the Gibbs phenomenon, whereas the reference solution shows what the solution in every point should be, based on the local shear rate. The shear rates range from Fig. 11 The stretch ratio in Poiseuille flow simulated using the standard formulation (red line) and derived from the logarithmic formulation (blue line) at time = 44.1 s. There are 50 elements in y-direction Fig. 12 The stretch ratio in Poiseuille flow simulated using the standard formulation (red circles) and derived from the logarithmic formulation (blue line) at time = 11.1 s. There are 200 elements in y-directioṅ γ = 10 −6 s −1 toγ = 3.75s −1 . Note that the minimum shear rate is not taken asγ = 0, because D would then be a zero tensor and R crit in Eq. 9 would go to infinity. Therefore,γ = 10 −6 s −1 has been taken as a value to investigate the morphology evolution at very low shear rates.
It is observed that the droplet radius corresponds very well to the reference solution (see Fig. 13) for N y = 1000. As expected, the high shear rates near the walls have resulted in the smallest droplet sizes due to capillary breakup. Closer to the center, there is a region that remains nearly the initial droplet size. In this region, the capillary number Ca is too small to cause breakup, yet too large to allow for enough film drainage, so coalescence is unsuccessful here. Even closer to the center, coalescence is observed to be successful, as the droplets become slightly larger than the initial droplet size. Lastly, in the center itself, the shear rate is so low that not many collisions take place, so coalescence is again not (very) successful there.
The solution of the droplet size is observed to converge linearly (see Fig. 14). However, the solution on a coarse mesh (N y = 50) shows erroneous behavior in the center of the domain (y = 2). A lot of coalescence appears to be occurring while none would be expected. This phenomenon occurs up to a certain element size (around the element size at N y = 200). Further mesh refinement at that point seems to give nearly linear convergence. Based on the linear interpolation used for the morphology variables, quadratic convergence would be expected, though this would only hold for smooth functions. The discontinuous solution is not smooth, for which only linear convergence can be achieved (Canuto et al. 1988). The solutions of the stretch There are 1000 elements in y-direction ratio at a long time scale were also compared for N y = 50 and N y = 1000 (see Fig. 15). The convergence plot for the stretch ratio (see Fig. 16) has been generated by comparing the simple shear flow reference solution to the Poiseuille flow solution averaged over the final 1000 time steps. The averaging is done, because there is no true steady state solution, but a dynamic equilibrium between necking and coalescence. The convergence appears to be nearly linear throughout the whole range of N y . To conclude this section on Poiseuille flow, we have motivated our description of the blend morphology with logarithmic quantities and we

Problem description and mesh convergence
In this section, we investigate the spatial and temporal evolution of the blend morphology in an eccentric cylinder flow, which can be considered as a simplified case of the complex flow in a dynamic mixer (Tjahjadi and Ottino 1991). The problem consists of one cylinder positioned off-center inside another larger cylinder. In principle, both cylinders can rotate independently, however this work focuses on the case where only the outer cylinder rotates with an angular velocity of . The outer cylinder has radius R outer and the inner cylinder has radius R inner (see Fig. 17).
The distance between the center of the inner and outer cylinder is defined as the eccentricity . The first example shows the case where = 1 s − 1 , R outer = 1 m, R inner = 0.3 m, and = 0.3 m. The outer cylinder rotates clockwise. The effect of the element size has been investigated with the same simulation parameters as described in the previous section. Meshes are generated using Gmsh (Geuzaine and Remacle 2009). The starting point is a mesh with 1082 triangular elements (see Fig. 18). Quadratic triangular elements are used, where the velocity is interpolated quadratically and the pressure and morphology variables     Table 2). The relative minimum element size h min is calculated as the square root of the smallest element area divided by R outer , where R outer is used as the global length scale of the eccentric cylinder flow problem. The same approach is followed for calculating the relative maximum element size h max . For this convergence study, we focus on the morphology at t = 80 s (the morphology no longer changes much over time) along the vertical line C1 between the top of the inner and outer cylinder (see Fig. 19 for the stretch ratio profile and Fig. 20 for the unstretched radius profile), as this line crosses both regions of interest of vortex flow and outer flow, and this line is relatively far away from the stagnation point which is not expected to yield useful information regarding mesh convergence. It is observed that M0 and M1 yield morphology results that are very different from the most refined mesh M4. M4 shows significantly different behavior near the boundary of the inner cylinder, which suggests some boundary layer effect. The remaining part of the solution of M4 seems to be approximated quite well by M2 and M3. Since the computational time of M2 is a factor of 16 times faster than the finest mesh and the region of the boundary effect seems very narrow, M2 is considered good enough for further simulations. For other geometries (different eccentricities ), meshes are generated with a relative minimum element size h min of O(10 −3 ).
To comment further on the boundary layer effect that is observed for M4, this raises an interesting point about how our blend model behaves when there are very small length scales, as are encountered in chaotic mixing flows. These length scales cannot be captured with a finite element mesh. As is shown in Figs. 19 and 20, M4 is observed to describe a thin boundary layer region that is smoothed by the coarser meshes. Although this example is not chaotic, we expect the numerical method to behave in a similar way for the small length scales in chaotic flow, meaning that small scale details are smoothed, but the numerical model will not fail. This would suggest that information is lost due to this smoothing. However, our macroscopic model does not contain information about the morphology at such small length scales to begin with, so the limitation is in the theoretical model itself and not in the numerical method. Basically, there is no point in simulating length scales smaller than the maximum droplet size, because smaller details are absorbed into the macroscopic description of the droplet populations. We do expect that, if the fine spatial layering would be captured, that breakup and coalescence would eliminate it. This is because two adjacent thin layers would experience approximately the same local shear rate and our model describes the final morphology according to the Grace curve, which breakup and coalescence would evolve towards. Droplet migration across streamlines has no effect in our model, because we assume all droplets to advect along the direction of the background flow field.
Other models by Muzzio et al. (1992) and Florek and Tucker (2005) are excellent for describing the stretching of the droplet morphology in globally chaotic flows. However, these models neglect interfacial tension and thereby omit breakup and coalescence events. The strength of our model lies in the ability to capture these topological changes.

Morphology evolution
The morphology evolution is shown for the example of = 1 s − 1 , R outer = 1 m, R inner = 0.3 m, and = 0.3 m. In the outer regions, the flow rotates clockwise, while in the inner region above the inner cylinder there is a clockwise rotating vortex (see Fig. 21 for the streamline pattern and Fig. 22 for the shear rate pattern). The shear rate is observed to be largest below the inner cylinder, while being small above the inner cylinder and at the lower region near the outer cylinder.
The morphology development over this flow field with initial stretch ratio β = 1 and unstretched radius R 0 = 10 − 5 m has been simulated (see Figs. 23 and 24). The time step is t = 0.01 s, the smoothing parameter is τ = 0.1 s, the Giesekus parameters are α G = 0.01 and τ G = 10 6 s. The simulation was run until time t end = 60 s. It is seen that after 0.16 rotations, the stretch ratio first grows near the bottom of the inner cylinder (where the shear rate is largest) and is advected clockwise around the vortex above the inner cylinder. During the stretching process, it makes sense that the unstretched radius of the droplets does not change. After 0.56 and 0.72 rotations, the distinctive flow fields within and outside of the vortex can clearly be observed. Droplets in both regions stretch and rotate within their own regions, where the droplets in the outer region stretch faster because of passing the region of large shear rates below the inner cylinder. Some of the droplets can already be observed to have broken up. There seems to be an outward moving breakup front. Eventually, after 1.43 rotations, the droplets in the inner vortex have also broken up. Afterwards, after around three rotations, the final morphology is reached.

Effect of
Simulations were done for three values of the offset of the inner cylinder, with the values = 0.2 m, = 0.3 m, and = 0.5 m. The other parameters are = 1 s − 1 , R outer = 1 m, and R inner = 0.3 m. The initial values are β = 1 and R 0 = 10 − 5 m. The mean droplet size and its standard deviation along curve C1 were plotted over time (see Fig. 25).
It is seen that the geometry with the inner cylinder closest to the outer cylinder produces the finest morphology, with the highest amount of uniformity and in the shortest amount of time. It makes sense that the finest morphology is obtained, because the narrow region between the two cylinders results in larger shear rates than in the other two flow problems. The high amount of uniformity is due to the large area of the vortex above the inner cylinder, which easily distributes droplets over a large portion of the domain. The short process time is also related to the large shear rates, which leads to droplets stretching faster and breaking up faster.

Effect of
Simulations were done for three values of the angular velocity of the outer cylinder, with the values = 0.5 s − 1 , = 1 s − 1 and = 2 s − 1 . The other parameters are = 0.3 m, R outer = 1 m, and R inner = 0.3 m. The initial values are β = 1 and R 0 = 10 − 5 m. The time evolution of the mean droplet size and its standard deviation along curve C1 is investigated (see Fig. 26).
Increasing appears to have a similar effect to increasing . The reason is the same as in the case of varying . For higher , larger shear rates are generated in the narrow region between the two cylinders. The higher also improves the uniformity of the morphology, as it makes the flow faster and therefore, the efficiency of mixing is improved. We have hereby shown our model to be useful for predicting the spatial and temporal development of a disperse blend morphology in a complex flow geometry.

Fig. 25
The mean unstretched droplet radius R 0 and its standard deviation over the line C1 for varying . The statistical data describe R 0 along C1 at a certain value of the time t

Fig. 26
The mean unstretched droplet radius R 0 and its standard deviation over the line C1 for varying . The statistical data describe R 0 along C1 at a certain value of the time t

Conclusions
We developed a numerical model to describe the evolution of the morphology of disperse blends in complex flows, using a Stokes flow and the model by Peters et al. (2001) as a starting point. Our phenomenological model was implemented in an inhouse FEM code. For validation, our model was first compared with the Peters model for simple shear flow, with good agreement. Filament breakup under nearly quiescent conditions was found to be an important morphology state, but has not yet been taken into account in our model, which is a topic for future work. The next step was to simulate the morphology development along the cross-section of a Poiseuille flow problem. In the Poiseuille flow example, the occurrence of the Gibbs phenomenon was observed. These Gibbs oscillations were found to possibly result in negative stretch ratios when β was simulated directly. This justifies the choice of calculating the stretch ratio with the variable s = log(β).
After simulating Poiseuille flow, the model was used to study the morphology development in the eccentric cylinder flow problem, where the outer cylinder rotates with angular velocity . The effect of the position of the inner cylinder and were investigated. Putting the two cylinders closer to each other and increasing were the best choices for obtaining the finest morphology, with the highest amount of uniformity in the shortest amount of processing time. In both cases, the cause was the increased shear rate in the narrow region between the inner and outer cylinder. While the original model by (Peters et al. 2001) has already been validated for simple shear flow with experiments done by Vinckier et al. (1997), for future work it would be interesting to experimentally validate the blend model for Poiseuille flow and eccentric cylinder flow.
We have now obtained a predictive tool for estimating how to design the complex geometry of a mixing device in order to obtain a desired disperse blend morphology, for different process conditions and material properties.