Numerical simulations of the polydisperse droplet size distribution of disperse blends in complex flow

The blend morphology model developed by Wong et al. (Rheologica Acta, 2019), based on Peters et al. (J Rheol 45(3):659–689, 2001), is used to investigate the development of the polydispersity of the disperse polymer blend morphology in complex flow. First, the model is extended with additional morphological states. The extended model is tested for simple shear flow, where it is found that the droplet size distribution does not simply scale with the shear rate, because this scaling does not hold for coalescing droplets. Subsequently, the model is applied to Poiseuille flow, showing formation of distinct layers, which occurs in realistic pressure-driven flows. Finally, the model is applied on an eccentric cylinder flow, where histograms are made of the average droplet size throughout the domain. It is observed that outer cylinder rotation results in narrow distributions where the small droplets are relatively large, whereas inner cylinder rotation results in broad distributions where the small droplets are significantly smaller than in the case of outer cylinder rotation. Eccentricity seems to only have a minor effect if the maximum shear rate is held constant. The flow profile and history in combination with the maximum shear rate strongly determine how the polydisperse droplet size distribution develops.


Introduction
A common method for creating polymer materials with targeted properties is to blend multiple homopolymers. The majority of polymer blends are immiscible, because mixing long polymer chains is thermodynamically unfavorable, leading to a multiphase structure (Lipatov 2002;Tucker and Moldenaers 2002). Depending on the volume fractions of the blend constituents, there can be a disperse or co-continuous morphology. In this work, we focus on modelling disperse blends. The morphology undergoes changes due to deformation, breakup and coalescence. generalized in the well-known Doi-Ohta model (Doi and Ohta 1991). Maffettone and Minale (1998) developed another model that describes the morphology with a shape tensor. Ellipsoidal drop shape predictions have been verified by three-dimensional visualization by Guido and Greco (2001). Droplet deformation and breakup have been investigated experimentally by Grace (1982), Bentley and Leal (1986b), Bentley and Leal (1986a), Stone et al. (1986), Grizzuti and Bifulco (1997), Vinckier et al. (1997), and Iza and Bousmina (2000). Coalescence has been modelled by Chesters (1991). 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). Based on Peters et al. (2001), Wong et al. (2019) have developed a numerical model to simulate the monodisperse blend morphology development in complex flow geometries.
For the material properties of the solidified polymer blend, it is not only important to know the droplet size, but also the droplet size distribution (polydispersity) (Wu 1988;Premphet and Paecharoenchai 2002;Caserta et al. 2004;Diop and Torkelson 2015). The objective of this study is to describe the development of polydisperse droplet size distributions in complex flow geometries using the work from Wong et al. (2019). We first introduce a number of extensions to this model and then use the extended model to investigate the shear rate dependence of the polydispersity in simple shear flow and the formation of distinct layers in Poiseuille flow. Finally, we use eccentric cylinder flow as an example of complex flow and show how polydispersity can be influenced through alterations to the geometry and flow protocol.

Morphology model
In this section, we begin by summarizing the features of our previous blend morphology model (Wong et al. 2019). This is followed by three extensions to our original model. Subsequently, we show the equations of this extended model. Finally, we show the set of equations as they are solved in our numerical framework.

Recap of the previous model
If we consider a blend consisting of a single droplet and assume the droplet to be ellipsoidal, we can quantify the model using the radius of the unstretched droplet R 0 , the stretch ratio β = L/B = L/2R 0 , where L is the major axis and B is the minor axis of the ellipsoidal droplet respectively, and the orientation vector m. In practical blends, the number of droplets is very large, which makes it computationally unattractive to track every droplet individually. We follow the approach of Peters et al. (2001) and model the morphology macroscopically. Droplets are locally grouped together in populations, that are described by an average R 0 , β and m, along with the number of droplets per unit volume N d . There is one droplet population for every spatial point x. The morphology over a whole geometry is described by macroscopic fields of the average population variables R 0 (x), β(x), m(x) and N d (x). These variables evolve under the influence of the background velocity field u(x). We calculate this flow field using the Stokes flow equations and assume that the morphology does not couple back to the flow field calculation. The morphology contributes an extra stress term to the equation for conservation of momentum. This extra stress can be split into an elastic term and a viscous term. The viscous term is taken into account by using an effective matrix viscosity instead of the pure matrix viscosity. We neglect the elastic term, and assume that this contribution is negligible in the cases that we are interested in, i.e., closed domains with no free surfaces. We aim to model the morphology development in practical mixers and in these situations, it is usually preferable to obtain as small as possible droplet sizes through rapid filament stretching. This implies that we are mostly dealing with very high capillary numbers and relatively small time scales. In the case of very high capillary numbers, where the stresses from the background flow field dominate the interfacial stresses of the microstructure, the interfaces tend to passively follow the flow. Interfacial stresses tend to become active at larger time scales, when droplets are very small. In this situation, the microstructure evolves relatively slowly. In our study of practical mixers, we are mostly not interested in these time scales. In short, in case of very high capillary numbers and relatively small time scales, we believe it is a valid assumption to neglect the elastic stress contribution of the microstructure on the momentum conservation equation.
In our original model, we first described the stretch ratio β and undeformed droplet radius R 0 with the following partial differential equations: where f 1 and f 2 are semi-empirical functions based on others' modelling and experimental work on droplet deformation, breakup and coalescence. Our model exhibits discontinuous jumps in directions perpendicular to streamlines. Numerically, this leads to the occurrence of the Gibbs phenomenon, where the solution shows oscillatory behavior around the discontinuity. Consequently, this may lead to negative values for R 0 and β, which is unphysical and makes the numerical simulations stop. We found that using the logarithmic variables v = log(R 0 ) and s = log(β) in the model and find R 0 and β by exponentiation, improved that stability of the model substantially (Wong et al. 2019). With the logarithmic variables, the partial differential equations become: with the extra division in the right-hand side coming from the chain rule of differentiating a natural logarithm. We describe the orientation vector m as the eigenvector corresponding to the largest eigenvalue of the conformation tensor c = b·b T , where b is the contravariant decomposition of c and is described by a Giesekus model, following the approach by Hütter et al. (2018). It can be thought of as something similar to the deformation gradient tensor from continuum mechanics with an additional relaxation term. The tensor c has a meaning similar to the Finger tensor from continuum mechanics. We calculate b using: which is precisely the well-known Giesekus-model, but written in terms of the b-tensor instead of the more conventional conformation tensor. The fourth morphology variable, the number of droplets per unit volume N d , is not described as an independent field variable, but is determined using conservation of volume based on R 0 : The right-hand side functions f 1 and f 2 depend on the morphological state, which depends on the capillary number. 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. In our case, the matrix surrounding a droplet contains other droplets. To take this into account, we replace μ m with an effective matrix viscosity μ e , as defined by Choi and Schowalter (1975): where φ is the volume fraction of the droplet phase in the total fluid mixture and λ is the viscosity ratio of the pure components. We define an effective viscosity ratio according to: where μ d is the viscosity of the droplet phase. In Wong et al. (2019), we took four morphological states into account, namely: coalescence, necking, filament stretching and dynamic filament breakup.

Modifications to the morphology model
In this work, we make three additions to the blend morphology model. The first addition is making a distinction between droplets and filaments. In the original model, the morphological state was only determined by the most recent local capillary number, so in theory a long filament could be advected into a lower shear rate region and instantaneously exhibit binary breakup (necking), which is not physical. We classify any droplet population with a stretch ratio β larger than 5 as filaments. This choice is made because, following Peters et al. (2001), we assigned β = 4 for the case of necking (Wong et al. 2019), so β > 5 can only be reached if there was Ca ≥ κCa crit (filament stretching) in the recent history. The second addition consists of two additional morphological states, namely static filament and static filament breakup. This had been missing from our original model. It is important to have a state for static filaments, in case the flow would be stopped. Suppose long filaments have been created through the state of filament stretching, which occurs at large capillary numbers. Stopping the flow would result in a capillary number of 0, which would make the droplet population undergo coalescence. Coalescence immediately following filament stretching would be unphysical.
The third addition to our original model is a more detailed classification of the shear and elongational components of a complex flow. In our previous work, we assumed κ = 2 and Ca crit to be only the shear flow part of the Grace curve. This was obviously done for simplicity, but this needs to be improved upon, because many practical flow problems involve significant elongational flow. Taking only the upper limit of the Grace curve (the shear flow part) will likely mispredict whether or not a droplet (population) will experience breakup.

Extended morphology model
A new step in the simulation procedure is to test whether a filament should be classified as stretching or quiescent. We quantify this using the stretching efficiency e f , which is defined according to Ottino (1989): where D = 1 2 (∇u + (∇u) T ). The stretching efficiency can be viewed as a ratio between the realized elongational rate and the maximum obtainable elongational rate when the filament would be perfectly aligned with the flow direction with the maximum deformation. This can be seen as follows. Firstly, the numerator D : mm corresponds to the elongational rate that a stretching droplet experiences in its axial direction. This dyadic product can be interpreted as follows: the tensor D is a three-dimensional mathematical object that contains all the spatial components of the strain rate. The dyadic product selects the strain rate component pointing in the direction of m on a surface of which the normal vector points in the direction of m, in other words, we have obtained the stretching rate in the axial direction of the elongated droplet. This direction does not necessarily have to be the direction of principle strain rate, which is the direction in which the flow field exerts the largest strain rate. The magnitude of this principle strain rate is equal to the denominator in the definition of the stretching efficiency, namely √ D : D. Therefore, the stretching efficiency is the ratio between the realized elongational rate and the maximum obtainable elongational rate.
In order to describe a quiescent filament (static filament) and filament breakup under quiescent conditions (static filament breakup), we introduce two additional field variables. These are the filament radius before breakup R f and the Rayleigh disturbance amplitude α f (the subscript "f" stands for "filament"), both of which determine when a quiescent filament becomes unstable and breaks up. As demonstrated in our previous work, discontinuities in the field variables in the direction perpendicular to the streamlines leads to the Gibbs phenomenon, where the solution shows oscillations around the discontinuity. These oscillations may lead to negative values for the field variables. Therefore, we describe R f and α f with logarithmic variables v f = log(R f ) and a f = log(α f ). They are described by the following partial differential equations: where f 3 and f 4 are right-hand side functions that depend on the morphological state. Note that there is no multiplication by 1 α f in Eq. (13). This is because the Rayleigh disturbance amplitude α f is described with an exponential function (Tomotika 1935), so the differential equation for the logarithmic a f is written in such a way that there is no division by itself in the right-hand side, in contrast with the other quantities β, R 0 and R f .
In our previous work (Wong et al. 2019), we assumed for simplicity that Ca crit and κ were always equal to the values for simple shear flow, as given by Grace (1982).
However, it is known that droplets cannot break up in shear flow for λ e > 4, whereas the same droplets could be broken up by extensional flow. This example demonstrates the importance of taking the Grace curve for (planar) extensional flow into account as well. The Grace curve gives the critical value Ca crit as a function of the viscosity ratio above which the deforming droplet becomes unstable and will eventually break up. The curves are given for the two extreme cases of simple shear flow and planar extensional flow. To reiterate, for simple shear flow, κ = 2 and the critical capillary number is given by De Bruijn (1989): and for planar extensional flow, κ = 5 and the critical capillary number is given by: In complex flow, we calculate the critical capillary number by interpolating between these two expressions. To do this, we need a quantity to describe the amount of shear and extension that is present in the velocity field. We use the ratio of the scalar deformation rate D s and scalar vorticity s , as defined by Hulsen (1988): where L = (∇u) T is the velocity gradient tensor. If this number is ζ = 1, shear flow is dominant. If ζ < 1, then extension is dominant, with ζ = 0 representing pure extension. Lastly, for ζ > 1, rotation is dominant, but we assume the Grace curve for simple shear to hold for this regime. In reality, when rotation is dominant, no Grace curve exists, because a rotating droplet does not deform and therefore does not break up. This also means that our assumption is not important, because when rotation is dominant, the local flow conditions will not require the Grace curve to be called upon. Using the quantity ζ , we interpolate complex flow types according to: As was mentioned before, the right-hand side functions f 1 , f 2 , f 3 and f 4 depend on the morphological state. The criteria for the six considered morphological states are shown in Fig. 1, which shows the solution procedure. The filament breakup criterion uses the filament radius R. For this, we assume that a stretched filament can be approximated as a cylinder, with equal volume as an initially spherical droplet with radius R 0 . It is defined as: The blue boxes indicate the original model from Wong et al. (2019) and the yellow boxes indicate the extensions from this work. In the following subsections, the right-hand side functions are shown for every morphological state individually.

Coalescence
In the case of coalescence, the equation for f 2 remains the same as in Wong et al. (2019). In this state, f 3 = 0 and f 4 = 0. We made a change to the equation for the logarithm of the stretch ratio s. The equations for coalescence are: with β + given by: where D = (L − B)/(L + B) is the dimensionless drop deformation parameter. 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 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 our previous work, we used the expression from Cox (1969): This expression is not valid for small values of the viscosity ratio λ e . It can be derived that for small λ e , D ≈ Ca. Coalescence occurs at Ca < Ca crit . According to Grace (1982), Ca crit 1 for small values of λ e . This means D ≥ 1 would be a common occurrence, which leads to a negative base in the exponent of Eq. (26). This is unphysical and we update our model with the more general expression given by Choi and Schowalter (1975): For small values of λ e , this results in: so: The result is smaller than 1, so there will not be a problem with Eq. (26) for small values of λ e .

Necking
For the case of necking, we use the same equations as in Wong et al. (2019), following Peters et al. (2001), with the addition of f 3 = 0 and f 4 = 0: with τ n a numerical smoothing parameter and

Filament stretching
For filament stretching, f 1 and f 2 are again the same as before, with the addition of f 4 = 0 and an additional expression for f 3 : with The expression for f 1 has been derived by Stegeman et al. (1999).

Dynamic filament breakup
Dynamic filament breakup is given by the same equations as in Wong et al. (2019) for the state of filament breakup, based on Tjahjadi and Ottino (1991): with τ b a numerical smoothing parameter and

Static filament
In the newly introduced state of a static filament, the droplet size and stretch ratio do not change. The only change occurs in the Rayleigh disturbance amplitude α f , so only f 4 is not zero: where m is the dimensionless disturbance growth rate as given by Tomotika (1935).

Static filament breakup
The second newly introduced morphological state of static filament breakup is described by the following equations: with where x m is the dominant dimensionless disturbance wavelength as given by Tomotika (1935) and τ s a numerical smoothing parameter set equal to the estimated time required for breakup of a filament under quiescent conditions (Janssen 1993):

Polydispersity
The goal of this work is to describe the development of polydisperse droplet size distributions in complex flows.
We model the distribution of droplet size with a probability density function p, as defined by: where V tot is the total volume of the disperse phase and V a→b is the volumetric fraction of droplets with a radius in the interval between R a 0 and R b 0 . We discretize this distribution into a number of bins, where the value of p is assumed to be constant within a bin width.

Probability distribution prediction using discrete bins
We start from a distribution of droplet sizes and assume there is a finite maximum R max 0 and minimum R min 0 that are not exceeded during flow. The total probability is initially: i.e., the area under the curve p init (R 0 ) is 1. As a result of the flow field, the distribution will be different, but we still assume all droplets to be within the finite interval where R 0 0 = R min 0 and R N 0 = R max 0 . Within a single interval we approximate the probability function with a constant value per interval p i , with i = 1, 2, ..., N. We aim for conservation of total probability, i.e., we aim for conservation of: Therefore, we initialize with: or: The algorithm for the transfer from p i init , i = 1, 2, ..., N to p i , i = 1, 2, ..., N must preserve total probability, which means:

Numerical procedure
We first select a range of droplet sizes R 0 that contains all the initial droplet sizes: [R min 0,init , R max 0,init ]. Then we divide this interval into n poly bin values. Note that N ≥ n poly , because all droplets that become smaller through breakup and larger through coalescence should be captured by the entire interval [R min 0 , R max 0 ] that we take into consideration. Every subpopulation is simulated separately according to the solution procedure as is shown in Fig. 1. They are assumed to be independent and not influencing each other. This means that different subpopulations do not interact, which is a strong assumption. We expect this to matter relatively little in the early stages, when droplets are mainly stretching into long filaments. These interactions are expected to be most relevant in the case of coalescence. Coalescence between unequally sized droplets can be described by using an average droplet size (Janssen 1993). We would then have to calculate interactions between every subpopulation with every other subpopulation, which we chose not to implement due to the large computational cost and because we are mostly interested in shorter time scales at which coalescence is relatively insignificant. After simulating every individual subpopulation, a certain initial distribution p init is chosen, and the available monodisperse simulations are then used to postprocess how the selected shape of p develops over time and space. The procedure for calculating the distribution p is as follows: -Start from a certain initial p init , for example, a constant value over a certain range of R 0 values. -Calculate the initial values for V i V tot in every bin. -For every time step, determine the new bins where every initial subpopulation has moved to and add the original value of the relative volume of the subpopulation to the new bin. -After all the relative volumes V i V tot have been added up, the new distribution p is calculated by dividing all the bin values of the relative volume by the bin interval lengths (R i+1 0 − R i 0 ). In this work, we always choose p init to be a constant value over the range specified by R min 0,init and R max 0,init and zero outside of this range. As should hold for any probability density function, the integral of p over the entire domain is always equal to 1, because the sum of all relative volumes is equal to 1.

Numerical implementation
The blend model has been implemented in a finite element framework, which solves the blend morphology in space and time according to the flow chart in Fig. 1. The velocity is solved with quadratic polynomial basis functions and the pressure is solved with linear polynomial basis functions according to Stokes flow equations. This is solved once in the beginning, since we assume steady state flow and no feedback from the blend morphology on the flow field. The deformation rate tensor is extracted from the velocity solution and plugged into the blend morphology model. The morphology variables are solved with linear polynomial basis functions. The time integration is first-order semiimplicit, which means that the right hand side functions f 1 to f 4 are evaluated explicitly using the known morphology quantities from the previous time step. The advection term is solved implicitly and stabilized using Streamline-Upwind Petrov-Galerkin (SUPG) test functions. The details are are described by Wong et al. (2019). This procedure is followed for every subpopulation independently at every time step, since we assume no interaction between them.

Influence of shear rate
We use the same simulation parameters (see Table 1) as in our previous work (Wong et al. 2019), as was taken from Peters et al. (2001).
To investigate the effect of the shear rateγ on the polydisperse droplet size distribution, we study an example of simple shear flow with n poly = 100 subpopulations in the range between R min 0,init = 10 −7 m and R max 0,init = 10 −4 m and three shear rate values ofγ = 5 s −1 ,γ = 25 s −1 anḋ γ = 125 s −1 . The 100 subpopulations are initialized with values of R 0 that are spaced logarithmic equidistant within the specified range.
In this study, we investigate whether increasing the shear rate merely speeds up the same polydispersity development or actually exhibits different kinetics. A higher shear rate is expected to speed up the droplet breakup processes, so to be able to compare the kinetics of the three shear rates, the time variable must be non-dimensionalized. By multiplying the time with the shear rate we obtain the dimensionless straiṅ γ t. The three shear rates are simulated for 160 strain units and histograms are made forγ t = 20,γ t = 40,γ t = 80 andγ t = 160 (see Fig. 2).
Comparing the three shear rates, there are a few similarities and differences. Forγ t = 20, it is observed that for all three cases, droplets that are slightly smaller or slightly larger than the Grace curve value quickly gather in a single bin that contains the Grace curve value.
Forγ t = 40, differences in the coalescence behavior are observed. In the case ofγ = 125 s −1 , no coalescence is observed, because the Grace curve value is smaller than the smallest droplet size of R 0 = 10 −7 m in the initial distribution, therefore all droplets exhibit breakup for this shear rate. In the case ofγ = 25 s −1 coalescence seems to occur relatively strongly, as all the droplets smaller than the Grace curve value rapidly enlarge and join the single bin containing the Grace curve value. For the case oḟ γ = 5 s −1 , coalescence is observed to occur much slower, which is expected because a lower shear rate results in a lower collision frequency of the small droplets.
Scaling of dR 0 dt withγ is not observed when coalescence occurs. By "scaling", we mean that increasing the shear rateγ by a factor of five does not speed up coalescence by a factor of five. The observation that the coalescence kinetics for a higher shear rate are not simply the scaled up version of those for a lower shear rate can be explained by the reciprocal shear rate 1/γ not being the correct timescale for non-dimensionalization in the case of coalescence. The shear rate is only the correct timescale when deformation is dominant over the interfacial tension (Ca ≥ Ca crit ), which is the opposite of the conditions for coalescence.
For bothγ t = 40 andγ t = 80, it is observed that the droplets that are slightly larger than the Grace curve value behave similarly. These droplets are in the necking regime, which scales with the shear rateγ , so this part of the distribution does indeed scale with the dimensionless strain unitγ t. Lastly, forγ t = 160, it is observed that the very large droplets still either have only decreased to a specific value or have remained the same size. This is a consequence of the description of stretching or static filaments in our model. The large droplets are all in the regime of Ca ≥ κCa crit , which means they initially undergo filament stretching. The relatively small droplets in this regime exhibit the conventional sequence of filament stretching, followed by dynamic filament breakup and necking towards the Grace curve value (see Fig. 3). The other droplets follow a pattern such as the example shown in Fig. 4. Filament stretching is cut short because the shear flow has rotated the droplet so much away from the direction of greatest deformation by the flow field that the stretching efficiency e f is considered low enough to enter the part of the flow chart of static filament and static filament breakup. This large rotation is caused by describing the deformation affinely (τ G is taken very large). We speculate that this too A limitation of our model is that the final droplet size that arises from static filament breakup can be larger than κR 0,crit . In reality, the droplets should again deform and break up, but our model does not capture this effect, because the b-tensor is never reset, so the amount of deformation and orientation vector m of these daughter droplets is not described accurately. Resetting the b-tensor and the effect this has on the morphology is beyond the scope of this paper. Because the orientation vector remains unchanged, the stretching efficiency also remains unchanged and our model can therefore not describe a second round of filament stretching. We do not propose a solution to this limitation at this point, because this issue does not seem to arise in our later simulations in complex flow (eccentric cylinder flow), so we assume that this is not an issue for flow problems that are not purely shear flow.

Poiseuille flow
Poiseuille flow is a model problem for practically relevant flow problems that are pressure-driven and exhibit a layered droplet morphology. We simulate a rectangular domain with N x = 2 and N y = 99 rectangular biquadratic elements, where −1 ≤ y ≤ 1 m. The flow rate isV = 10 3 m 2 /s, which results in shear rates in the range of 0 ≤γ ≤ 5 [s −1 ]. To demonstrate the layeredness of the polydisperse morphology, we visualized the droplet size distribution with n poly = 100 along the cross-section (the y-direction of the domain (see Fig. 5a) at time t = 12 s, with R min 0,init = 10 −7 m and R max 0,init = 10 −5 m. Several distinctive regimes can be recognized. In the core, there is a region where Ca > Ca crit , where coalescence occurs, but this is a slow process, so not much change is observed at such a small time scale. Slightly beyond the dashed curve, which represents Ca = Ca crit according to the local shear rate, some droplets have begun to breakup (decrease in size) due to necking. When looking at even larger values of y, a flat region is observed where droplets are stretching into filaments. At the largest values of y (near the walls of the domain), a sharp decrease of droplet size is observed due to dynamic filament stretching.
Histograms have been made of the droplet size distributions at three positions, namely y = 0.1 m, y = 0.5 m and y = 0.9 m (see Fig. 5). At y = 0.1 m, the majority of the population has not changed significantly at this time scale. At y = 0.5 m, it is seen that the largest droplets have not yet changed in size, so are still undergoing filament stretching, while the smaller droplets that are closer to the Grace curve have started to accumulate around the predicted value from both sides via coalescence and necking. At y = 0.9 m, it is seen that the largest droplets have all broken up via filament stretching, which is a fast process and it is also seen that the smallest droplets at this position are larger than those in the distribution at y = 0.5 m, because here smaller droplets are experiencing coalescence under the influence of a larger shear rate. This Poiseuille example has demonstrated how small variations in the shear rate profile can lead to a drastically different layered morphology in a pressure-driven flow situation.

Problem description and model validation
The eccentric cylinder flow problem consists of a smaller cylinder with radius R inner positioned off-center within a larger cylinder with radius R outer , with the distance between their center points being given by the eccentricity (see Fig. 6). The two cylinders can rotate independently, with angular velocities of inner and outer , respectively. We first use this geometry to do a qualitative validation of our blend morphology model for a complex flow problem that combines shear and extensional flow. To this end, we make use of the experiments performed by Tjahjadi and Ottino (1991) on droplets stretching and breaking up in an eccentric flow setup. It is technically an example of distributive mixing, whereas we study dispersive mixing, but we chose this example of a flow that is more complex than simple shear flow to qualitatively validate where the large and small droplets appear and the effect of varying the viscosity ratio. The matrix phase is corn syrup 1632 with a dynamic viscosity of μ m = 32.9 Pa s. The droplet phase is a mixture of No. 40 oxidized castor oil and 1-Bromonaphtalene, mixed in varying ratios to obtain a dynamic viscosity μ d for a specified viscosity ratio, λ = 0.010 for case 1 and λ = 0.40 for case 2. We use for our two cases the same droplet and pure matrix phase viscosities as Tjahjadi and Ottino (1991), but it has to be noted that our calculations actually do not use the pure matrix viscosity, but the effective matrix viscosity μ e . The geometry is described by the outer cylinder radius R outer = 7.62 cm, inner cylinder radius R inner = 2.54 cm and eccentricity = 2.29 cm. Tjahjadi and Ottino (1991) use one initial drop with a radius of approximately 0.5 cm. In our model we cannot look at a situation with only one drop, and we initialize across the entire domain monodisperse droplet populations with an average droplet size of R 0,init = 0.5 cm. The applied flow protocol consists of a specified number of periods, where a period consists of: first 1 4 period clockwise rotation of the inner cylinder with inner = − 1 2 s −1 , followed by 1 2 period counterclockwise rotation of the outer cylinder with outer = 1 6 s −1 , continuing with another 1 4 period clockwise rotation of the inner cylinder with inner = − 1 2 s −1 . We simulate two periods, with one period corresponding to 150 s. The results for the case of λ = 0.010 are shown in Fig. 7 and for the case of λ = 0.40 are shown in Fig. 8.
Qualitative similarities are nicely observed between the experiment and simulation. The dark blue regions in the simulations with the small droplets correspond almost exactly with the regions in the experiments that contain the small droplets generated from filament breakup. These regions contain particle trajectories along which relatively large shear rates are experienced throughout the flow history. On the other hand, the red regions in the simulations with the large droplets correspond almost exactly with the regions in the experiments that are completely devoid of droplets. These regions experience mostly small shear rates and therefore, breakup does not occur there. In the experiments, these regions are empty, because the timescale is too short for any droplets to migrate there. It is also interesting to note that our simulations successfully capture Fig. 7 Comparison between the experiment by Tjahjadi and Ottino (1991) (Figure 13, with permission from J. Fluid Mech.) and the simulation of the unstretched droplet size R 0 in (m) with our model after two periods for viscosity ratio λ = 0.010. (a) Experiment. (b) Simulation the quality that λ = 0.40 results in much smaller droplets than λ = 0.010, which is to be expected based on the Grace curve that predicts smaller critical capillary numbers Ca crit , which is the value separating coalescence from breakup, around λ = 1. With this, we have qualitatively validated our blend model for an example of a complex flow.

Polydispersity and element-area weighted histograms
In the subsequent sections, we investigate the distribution of the mean value of R 0 throughout the whole domain, with R inner = 0.03 m, R outer = 0.1 m and = 0.03 m, except in "Influence of ". We do that by generating the probability density p in every node in the mesh, then calculating the mean R 0 in every node, then finally generating histogram data of these nodal values to find the global mean value of R 0 for every time step. However, this is not a fair comparison, because the computational mesh in such problems is usually not uniform, but denser near the inner cylinder and coarser near the outer cylinder.
Therefore, many more nodes from the inner part of the domain are sampled for the histogram data. We try to mitigate this issue by weighing the nodal values using the area of the surrounding mesh elements. The resulting histogram data is more mesh-independent, because the outer region has fewer nodes than the inner region. The weighted histograms represent a more uniform sampling of the morphology variables throughout the entire domain.

Influence of n poly
In this section, we investigate the influence of the value of n poly on the polydisperse droplet populations. We simulate the case of rotating the outer cylinder with outer = 1 s −1 , with R 0 initially distributed logarithmically between R min 0,init = 10 −7 m and R max 0,init = 10 −4 m. For the evolving polydispersity, we focus on the point midway between the bottom of the outer and inner cylinder. This region has the highest shear rates, so it is expected that most dynamical processes occur here. The distribution of R 0 after 10 rotations of the outer cylinder is shown in Fig. 9. We cannot compare the distributions with an exact solution and it is computationally extremely expensive to generate a sufficiently fine distribution to serve as pseudoexact solution, so we compare the volume-averaged R 0 within this sample point (see Fig. 10).
It is observed that taking n poly = 100 subpopulations appears to be good enough to describe the polydisperse

Influence of outer
In this section, we study the influence of the rotational speed of the outer cylinder outer on the distribution of the mean R 0 of the polydisperse droplet populations throughout the domain. The mesh consists of 35,008 nodes and 17,312 triangular elements, with quadratic interpolation for the velocity and linear interpolation for the pressure and blend morphology variables. The case of outer = 1 s −1 is compared to the case of outer = 2 s −1 and outer = 0.5 s −1 . The temporal evolution of the global mean value of R 0 is shown in Fig. 11, with weighted histograms of four values of the time shown in Fig. 15.
As expected, it is seen that droplets break up faster for higher outer , due to the larger shear rates throughout the domain for the same flow pattern. Larger shear rate values lead to larger values for the capillary number Ca, which leads to faster stretching and thinner filaments before breakup, resulting in smaller droplet sizes after breakup. Another consequence of a larger outer is that advection is faster, so the influence of regions with large shear rates is propagated faster throughout the domain. From the four histograms, it is observed that the last three frames (t ≥ 90 s) of outer = 2 s −1 are almost identical, so as expected, for Fig. 11 Global mean of R 0 as function of time for varying outer large outer , the droplet size distribution quickly reaches its final state, which is determined by the values of Ca crit . It is important to note that the entire distribution consists of relatively small droplets. This is in contrast with the reference case of outer = 1 s −1 , where it is clearly seen that at first (t = 90 s), a relatively wide distribution is created. The bigger droplets still need more time to stretch and break up, but eventually they do break up and the distribution narrower. Interestingly, in the case of outer = 0.5 s −1 , the width of the distribution seems to be nearly constant in the later stages of the process (t ≥ 90 s). This is because the shear rates are so low that it takes a very long time to break up the larger droplets, while the same low shear rates also cannot break up droplets to sizes as small as those of outer = 1 s −1 and outer = 2 s −1 . Therefore, the range of the distribution stays relatively the same, but the mean value is constantly decreasing over time.

Influence of inner
In this section, we study the influence of the rotational speed of the inner cylinder inner on the distribution of the mean R 0 of the polydisperse droplet populations throughout the domain. The case of outer = 1 s −1 is compared to the case of inner = 1 s −1 and inner = 2.815 s −1 . The value of inner = 2.815 s −1 was chosen such that the maximum shear rate corresponds exactly to the maximum shear rate in the reference case. The temporal evolution of the global Fig. 14 Global mean of R 0 as function of time for outer = 1 s −1 , inner = 3 s −1 and the time-periodic alternating case between outer = 1 s −1 and inner = −3 s −1 mean value of R 0 is shown in Fig. 12, with weighted histograms of four values of the time shown in Fig. 16.
It is seen that inner = 1 s −1 is very weak and takes a very long time for any significant change to take place. In the case of inner = 2.815 s −1 , it is seen that a very broad distribution is created, though a significant fraction of the distribution remains on the large side near the initial value. The large droplets that do not seem to break up are in the regions with low shear rates. The flow profile in the case of only rotating the inner cylinder is such that droplets are advected along trajectories with almost constant shear rates. Therefore, droplets that begin in low shear rate regions will never experience high shear rates anywhere along their particle trajectories, so they do not experience any change. The shear rates are too low to deform the droplets, but the droplets are too large for coalescence. The droplets that do break up appear to form a very broad and relatively flat distribution. This is because this flow profile results in an almost linear progression of the shear rate away from the inner cylinder, with the highest values near the inner cylinder and and the lowest values near the outer cylinder.
Other simulations that are not shown here have shown this same trend, where rotating the inner cylinder results in a relatively broad and flat distribution with a peak near the initial value, and where rotating the outer cylinder results in a relatively narrow distribution around a certain droplet size. The minimum droplet size resulting from rotating the inner cylinder is much smaller than from rotating the outer cylinder, because those droplets are formed in the region of highest shear rate and due to the flow profile, they never experience smaller shear rates, so will not experience any significant coalescence to become larger droplets.

Influence of
In this section, we study the influence of the eccentricity of the inner cylinder on the distribution of the mean R 0 of the polydisperse droplet populations throughout the domain. The mesh of the two cases with = 0.05 m consists of 32,260 nodes and 15,938 quadratic triangular elements. The case of outer = 1 s −1 and = 0.03 m is compared to the case of outer = 1 s −1 and = 0.05 m and outer = 0.5396 s −1 and = 0.05 m. The value of outer = 0.5396 s −1 was chosen such that the maximum shear rate corresponds exactly to the maximum shear rate in the reference case. The temporal evolution of the global mean value of R 0 is shown in Fig. 13, with weighted histograms of four values of the time shown in Fig. 17.
At first sight, it may seem that increasing the eccentricity leads to faster break up and smaller droplets, but this is in fact a result of the increased maximum shear rate that results from this change to the flow geometry. This is illustrated by the results from outer = 0.5396 s −1 and = 0.05 m, which is the case where the maximum shear rate is exactly matched to the reference case. In this case, the histograms at all times are very similar to the reference case in shape and magnitude. This suggests that the distribution of the mean R 0 of the polydisperse populations is mostly determined by the combination of the maximum shear rate and the flow profile.

Influence of a chaotic flow protocol
In this section, we study the influence of alternatingly rotating the outer and inner cylinder on the distribution of the mean R 0 of the polydisperse droplet populations throughout the domain. The case of outer = 1 s −1 is compared to the case of inner = 3 s −1 and a timeperiodic case that alternates between outer = 1 s −1 and inner = −3 s −1 . A positive rotational speed is defined as counterclockwise rotation and a negative rotational speed is defined as a clockwise rotation.
The flow protocol is the same as the one given by Tjahjadi and Ottino (1991), but with higher values for the rotational speeds. We define the flow field according to inner = −3 s −1 for 0 ≤ t ≤ 37.5 s, outer = 1 s −1 for 37.5 < t ≤ 112.5 s, inner = − 3 s −1 for 112.5 < t ≤ 187.5 s, inner = 1 s −1 for 187.5 < t ≤ 262.5 s and inner = − 3 s −1 for 262.5 < t ≤ 300 s. The temporal evolution of the global mean value of R 0 is shown in Fig. 14, with weighted histograms of four values of the time shown in Fig. 18.
The time-periodic protocol seems to take the best aspects of both individual flow profiles. From the outer rotating cylinder, it takes the property of breaking the larger droplets earlier than in the case of the inner rotating cylinder. From the inner rotating cylinder, it takes the property of creating a much smaller minimum droplet size compared to the case of the outer rotating cylinder. In the last frame, it is observed that the chaotic flow protocol has the very small droplet sizes without the stagnant large droplets (Figs. 15,16,17 and 18).

Conclusions
The extended blend model was first applied to simple shear flow, where it was found that the polydispersity distribution does not simply scale with the shear rate, because this scaling does not hold for the small droplets that exist in the coalescence regime. The model was then applied to Poiseuille flow, showing formation of a layered blend morphology. Subsequently, the model was applied on eccentric cylinder flow, where histograms were made of the average droplet size throughout the domain. It was observed that outer cylinder rotation results in narrow distributions where the small droplets are relatively large, whereas inner cylinder rotation results in broad distributions where the small droplets are significantly smaller than in the case of outer cylinder rotation. Outer cylinder rotation broke up droplets sooner than inner cylinder rotation if the maximum shear rate was held constant. Eccentricity did not seem to have any significant effect if the maximum shear rate was held constant. A timeperiodic chaotic flow protocol, based on Tjahjadi and Ottino (1991), turned out to break up the large droplets effectively like outer cylinder rotation, yet also obtain the very small droplets after breakup like inner cylinder rotation does, seemingly avoiding coalescence of the very small droplets. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.