<figure><img src="image_1.png" /></figure>

HAL
open science

# Microstructure evolution during liquid-liquid laminar mixing: a kinetic theory approach

Francisco Chinesta, Malcom Mackley

► To cite this version:

Francisco Chinesta, Malcom Mackley. Microstructure evolution during liquid-liquid laminar mixing: a kinetic theory approach. International Journal of Material Forming, 2008, 1 (1), pp.47 - 55. 10.1007/s12289-008-0007-2 . hal-01007049

HAL Id: hal-01007049

https://hal.science/hal-01007049

Submitted on 2 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Microstructure evolution during liquid-liquid laminar mixing: a kinetic theory approach

F. Chinesta • M. R. Mackley

**Abstract** Many material forming processes involve liquid-liquid mixing, and in general the induced properties are strongly dependent on the resulting microstructure. Microstructure during liquid-liquid mixing exhibits a morphology characterized by a length scale usually much smaller than the one associated with the macroscopic flow. In this context a microstructure description allowing to characterize the microscopic morphology, its characteristic size, shape and orientation as well as the time evolution of the specific interface area, seems to be necessary in order to qualify and even quantify, the ability of flows to perform mixing, leading to the definition of optimal flows to maximize some desired criteria. The approach based on the definition of area tensor is a promising description of such phenomena, being its main drawback the necessity of introducing a closure relation to derive the equation governing its time evolution, whose impact on the computed solution can be in some cases significant. In this paper we propose a new description which considers the area tensor description as starting point, but defines its evolution in a kinetic theory framework, avoiding the introduction of any closure relation.

**Keywords** Mixing · Microstructured fluids · Kinetic theory description · Area tensor

F. Chinesta (✉)
LMSP UMR 8106 CNRS-ENSAM-ESEM,
151 Boulevard de l'Hôpital,
F-75013 Paris, France
e-mail: francisco.chinesta@paris.ensam.fr

M. R. Mackley
Department of Chemical Engineering, University of Cambridge,
Pembroke Street,
Cambridge CB2 3RA, UK
e-mail: mrm5@cam.ac.uk

## Introduction

Many engineering processes benefit from good bulk mixing. In such processes, determination of mixing rate is of importance in terms of understanding and prediction of mixing time. That understanding allows defining new flows and the associated processes for improving the mixing rate, but one could also try to optimize other parameters related to the microstructure as for example the characteristic length, the shape and orientation of the discrete phase, ...

Mixing rates have often been quantified by a mixing time, i.e. the time required to reach a certain degree of homogeneity (see [13]). Most studies have based on the concentration evolution at a small number of sample points [6]. However the sensibility of the results to the number and locations of the probes is its main limitation. An alternative way of quantifying the mixing rate can be derived from the consideration of only one of the basic mechanisms of mixing processes: the increase of the material interface due to the fluid mechanics in absence of interfacial tension and molecular diffusion. The molecular diffusion leads to smooth concentration gradients across the interface. However this mechanism only becomes significant when the interface has increased significantly. It is therefore expected that the overall mixing rate will be closely linked to the rate of mechanical stretching of the interfacial area. This mechanism has been widely analyzed [11, 14]. This approach was revisited in Mackley and Neves Saraiva [9] using a Lagrangian and concentration based numerical approach.

The approaches quantifying the mixing from the increase of the material interface have two important drawbacks: (i) the first one is related to the difficulty of introducing other additional physics as the one related to the surface tension, and (ii) sometimes the microstructure description needs
other information (morphology, characteristic length, shape and orientation of the discrete phase —inclusions-, ...) that the interface evolution cannot provide. The evolution of the morphology can then be modeled by predicting the change in the local morphological measure due to this velocity or deformation gradient. We will call such an approach, which treats some local characteristic morphological measure as a field variable, a micromixing analysis.

The micromixing approach has been used extensively for modeling passive mixing, where interfacial energy is negligible and the two phases have identical viscosities. In passive mixing the global velocity field can be found independently of the microstructure and then used to evolve the mixture structure, described with some area tensor, very rich from the morphological and microstructural points of view. See Wetzel and Tucker III [15] for an excellent description of this approach, in which the genesis of this kind of approaches is retraced precisely. The main difficulties noticed in that work [15] were related to the necessity of introducing a closure relation in order to derive an evolution equation for the second order area tensor. Different closure relations have been proposed, being their impact, in general, unpredictable.

In this paper we propose a new modeling based on the use of the area tensors, but that uses a kinetic theory description avoiding the introduction of closure relations. This approach is very close to the original one proposed by Doi and Ohta [5] but differs in the way of incorporating additional physics, as for example the surface tension. In the original Doi and Ohta model the surface tension is directly introduced in the preaveraged model which makes use of the orientation tensors and which needs for the introduction of a closure relation (the quadratic closure was considered in Doi and Ohta [5]). In the technique here proposed, the physics introduced in the preaveraged model is transferred again to the kinetic theory scale that allows computing the solution of the model without the introduction of any closure relation.

In what follows, only the case of passive mixing (fluids with same densities and viscosities) will be addressed. The case of active mixing is a work in progress.

## Description of microstructure morphology from the use of area tensors

We summarize in this section the main concepts related to the morphology description of microstructured fluids from the use of the area tensors, which define the basis of the kinetic theory approach that will be elaborated in next section. The interested reader can refer to Wetzel and Tucker [15] and the references therein for more details. Let's Ω be the domain in which the flow problem is

defined. Points in Ω will be referred by x which consists in a vector in 2D or 3D. In the numerical examples here considered only 2D models will be considered because in such case exact solution can be sometimes easily obtained. In order to quantify the morphology at any point x∈Ω a microscopic representative volume V(x) is considered centered at that point. This volume results small with respect to the macroscopic scale (related to the variation of the velocity field in Ω) but large enough with respect to the characteristic size of the microstructure. Let's S(x) be the interface within V(x). The area tensor is then defined as:

$$A(\mathbf{x}) = \frac{1}{V(\mathbf{x})} \int_{S(\mathbf{x})} \mathbf{n} \otimes \mathbf{n} \, dS \quad (1)$$

where n represents the unit vector defined on the interface S(x) and that we assume is pointing towards the continuous phase (the discrete one being the one having lower volume fraction). The symbol ⊗ denotes the tensor product, i.e. the components of the tensor product are given by: $(n \otimes n)_{i,j} = n_i n_j$. This area tensor is symmetric and has different appealing properties that we are summarizing in what follows:

1. The first property concerns its trace (sum of the diagonal components of A(x)) that we symbolize by $\mathrm{Tr}(A(x))$ that taking into account the normality of n results:

$$ \begin{aligned} \mathrm{Tr}(A(\mathbf{x})) &= \frac{1}{V(\mathbf{x})} \int_{S(\mathbf{x})} (\mathbf{n}_1^2 + \mathbf{n}_2^2 + \mathbf{n}_3^2) dS \\ &= \frac{1}{V(\mathbf{x})} \int_{S(\mathbf{x})} dS = \frac{S(\mathbf{x})}{V(\mathbf{x})} = S_v(\mathbf{x}) \end{aligned} \quad (2) $$

where $S_v(x)$ represents the specific surface related to point x, and whose maximization is usually searched in mixing processes.

2. If we define the volume fraction of the disperse phase as φ then we can define a characteristic length of the microstructure at point x from the volume of the discrete phase $V_d(x) = \phi V(x)$:

$$L(\mathbf{x}) = \frac{V_d(\mathbf{x})}{S(\mathbf{x})} = \frac{\phi V(\mathbf{x})}{S(\mathbf{x})} = \frac{\phi}{S_v(\mathbf{x})} = \frac{\phi}{\mathrm{Tr}(A(\mathbf{x}))} \quad (3)$$

3. In passive mixing the strain rate varies linearly at the scale of V(x), so that an initially spherical discrete domain will deform into an ellipsoid for any state of strain. This characteristic ellipsoid provides a convenient way to interpret the area tensor. The microstructure
shape and orientation can be easily deduced from the
normalized area tensor $\tilde{A}(x)$ defined as

$$
\tilde{A}(\mathbf{x}) = \frac{\mathbf{A}(\mathbf{x})}{\mathrm{Tr}(\mathbf{A}(\mathbf{x}))} \qquad (4)
$$

and whose eigenvalues allow computing the length of
the ellipsoid axes, being their orientation given by the
associated eigenvectors. See Wetzel and Tucker [15] for
more details concerning the relation between that
ellipsoid and the area tensor.

Area tensor evolution equation

The interest of computing the area tensor, from a given
initial condition, at each point of the flow domain and at
each time has been justified in the previous section.
However, its evaluation needs the establishment of the
partial differential equation governing its evolution. This
section focuses in the establishment of this evolution
equation and the closure problematic that its solution
involves.

To derive the evolution equation of $\mathbf{A}(\mathbf{x})$ we consider the
time derivative of Eq. 1 taking into account that both $\mathbf{n}$ and
$d\mathbf{S}$ evolve in time according to:

$$
\frac{dn}{dt} = \dot{n} = -(\mathrm{Grad}v)^{\mathrm{T}}n + (\mathrm{Grad}v:(n \otimes n))n \quad (5)
$$

expression initially proposed by Chella and Ottino [3]; and

$$
\frac{dS}{dt} = \dot{dS} = -(\mathrm{Grad}v:(\mathbf{n} \otimes \mathbf{n}))dS \quad (6)
$$

where “:” denotes the tensor product twice contracted, i.e.
**a**: **b** = Σ<sub>i</sub> Σ<sub>j</sub> **a**<sub>i,j</sub> ⋅ **b**<sub>i,j</sub>.

Introducing Eqs. 5 and 6 into the time derivative of Eq. 1
it results:

$$
\frac{d\mathbf{A}}{dt} = -(\mathbf{Grad}v)^{\mathrm{T}}\mathbf{A} - \mathbf{A}\mathbf{Grad}v + \mathbf{Grad}v : \mathbf{A} \quad (7)
$$

where the dependence of $\mathbf{A}(\mathbf{x})$ on $\mathbf{x}$ has been omitted for the
sake of clarity, and where $^4\mathbf{A}$ denotes the fourth order area
tensor defined by:

$$
{}^4\mathbf{A}(\mathbf{x}) = \frac{1}{V(\mathbf{x})} \int_{S(\mathbf{x})} \mathbf{n} \otimes \mathbf{n} \otimes \mathbf{n} \otimes \mathbf{n} \, dS \quad (8)
$$

The main difficulty in solving Eq. 7 is precisely related
to the existence of that fourth order area tensor. Thus, a
closure relation expressing that tensor as a function of the
second order area tensor is needed. As proved in Wetzel
and Tucker [15] considering a finite strain formulation, no
general exact closure relation exists, and then the validity of
any proposal must be checked carefully because its impact
on the computed solution is almost unpredictable.

Illustrative example

In this section we are considering a simple shear flow
where both kinematics and microstructure are assumed
defined in a two-dimensional space. Thus, the flow
kinematics is defined from:

$$
\mathbf{v} = \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} \dot{\gamma} y \\ 0 \end{pmatrix} \qquad (9)
$$

The surface orientation vector **n** results in the polar
coordinates system:

$$
\mathbf{n} = \begin{pmatrix} \cos \varphi \\ \sin \varphi \end{pmatrix} \qquad (10)
$$

Thus, Eq. 5 reduces to:

$$
\dot{\varphi} = -\dot{\gamma} \cos^2(\varphi) \qquad (11)
$$

proving the ability of this flow to align the interface on the
shear direction, which results in $\varphi = \frac{\pi}{2}$ and consequently
$\dot{\varphi}(\varphi = \frac{\pi}{2}) = 0$. Now, if the quadratic closure relation
(introduced by Doi and Ohta [5] and analyzed in Wetzel
and Tucker [15]) is assumed

$$
{}^4\mathbf{A} = \frac{1}{S_v} \mathbf{A} \otimes \mathbf{A} \Rightarrow {}^4\mathbf{A}_{i,j,k,l} = \frac{1}{S_v} \mathbf{A}_{i,j} \times \mathbf{A}_{k,l} \quad (12)
$$

then the evolution equation of the area tensor, Eq. 7, taking
into account that Gradv : $^4\mathbf{A} = \frac{1}{\mathrm{Tr}(\mathbf{A})}(\mathrm{Grad}v : \mathbf{A})\mathbf{A} =
(\mathrm{Grad}v : \mathbf{A})\tilde{\mathbf{A}}$, results:

$$
\begin{pmatrix} \bullet \\ \mathbf{A}_{11} \\ \mathbf{A}_{12} \\ \mathbf{A}_{22} \end{pmatrix} = \begin{pmatrix} \bullet & \gamma \frac{A_{12}A_{11}}{A_{11}+A_{22}} \\ \bullet & \gamma \left( \mathbf{A}_{11} - \frac{A_{12}^2}{A_{11}+A_{22}} \right) \\ \bullet & \gamma \left( 2\mathbf{A}_{12} - \frac{A_{12}A_{22}}{A_{11}+A_{22}} \right) \end{pmatrix} \quad (13)
$$

It is easy to prove that the quadratic closure relation
becomes exact when the interface within the representative
volume is fully aligned in a certain direction. Let's $\hat{\mathbf{n}}$ be
such direction. In this case, the area tensor reduces to

$$
\mathbf{A} = \frac{1}{V} \int_S \hat{\mathbf{n}} \otimes \hat{\mathbf{n}} \, dS = \frac{S}{V} (\hat{\mathbf{n}} \otimes \hat{\mathbf{n}}) = S_v (\hat{\mathbf{n}} \otimes \hat{\mathbf{n}}) \quad (14)
$$

whereas the fourth order one results

$$
\begin{align}
{}^4\mathbf{A} &= \frac{1}{V} \int_S \hat{\mathbf{n}} \otimes \hat{\mathbf{n}} \otimes \hat{\mathbf{n}} \otimes \hat{\mathbf{n}} \, dS = \frac{S}{V} (\hat{\mathbf{n}} \otimes \hat{\mathbf{n}} \otimes \hat{\mathbf{n}} \otimes \hat{\mathbf{n}}) = \nonumber \\
&= S_v (\hat{\mathbf{n}} \otimes \hat{\mathbf{n}}) \otimes (\hat{\mathbf{n}} \otimes \hat{\mathbf{n}}) = \frac{1}{S_v} \mathbf{A} \otimes \mathbf{A} \tag{15}
\end{align}
$$

which proves that in this case the quadratic closure relation
becomes exact. The numerical integration of Eq. 13 leads to
the specific surface growing in time, except when the
surface is initially aligned in the shear direction because in that case the time derivative of all the components of the area tensor vanishes, and then the initial solution is the solution at any time. From the physical point of view this fact is obvious because an interface located on the shear direction is not stretched. For any other initial orientation the specific surface is increasing the more and more and no steady state is reached. This result is also obvious because one could expect that an interface element initially oriented out the shear direction will be infinitely stretched because it takes an infinite time to orient in the shear direction. We are proving this reasoning. Consider an interface element whose normal vector is defined at the initial time $t=0$ by an angle $\varphi > \frac{\pi}{2}$. The shear flow that we are considering tend to orient this element in the flow direction (local alignment) defined by $\varphi^{l.a.} = \frac{\pi}{2}$, at the same time as it is stretched. The time $t$ required to reach a full alignment in the flow direction results from:

$$ \frac{d\varphi}{dt} = -\dot{\gamma} \cos^2(\varphi) \Rightarrow \frac{d\varphi}{\cos^2(\varphi)} = -\dot{\gamma} dt \\ \Rightarrow \tan(\varphi^{l.a.}) - \tan(\varphi) = -\dot{\gamma} t \quad (16) $$

which implies:

$$ t = \frac{-\tan(\varphi^{l.a.}) + \tan(\varphi)}{\dot{\gamma}} = \infty \quad (17) $$

This result proves that an infinite time is needed to reach a fully local alignment of the interface elements on the flow direction. Now, we are proving that during this time the interface is infinitely stretched. For this purpose we start from the expression of the surface stretching in the finite strain framework:

$$ \mathbf{n} \, dS = \det(\mathbf{F}) \mathbf{F}^{-T} \mathbf{n}^{(0)} dS_0 \quad (18) $$

where $\mathbf{F}$ represents the gradient of deformation tensor (describing the mapping between the reference and the present coordinates), $\mathbf{n}^{(0)}$ the unit vector defining the initial area $dS_0$. We assume the initial elementary area normal to the flow direction $(\mathbf{n}^{(0)})^T = (1 \ 0)$. For computing the surface stretching only the expression of $\mathbf{F}$ is needed, that can be computed by solving the equation governing its evolution:

$$ \begin{aligned} \dot{\mathbf{F}} &= \mathrm{Grad}\mathbf{v}\mathbf{F} = \begin{pmatrix} 0 & \dot{\gamma} \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \mathbf{F}_{11} & \mathbf{F}_{12} \\ \mathbf{F}_{21} & \mathbf{F}_{22} \end{pmatrix} \\ &= \dot{\gamma} \begin{pmatrix} \mathbf{F}_{21} & \mathbf{F}_{22} \\ 0 & 0 \end{pmatrix} \end{aligned} \quad (19) $$

As at the initial time the present coordinates coincides with the initial ones, the initial condition for $\mathbf{F}$ results:

$$ \mathbf{F}(t=0) = \mathbf{F}^{(0)} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \quad (20) $$

Thus, the solution of Eqs. 19 and 20 writes:

$$ \left. \begin{aligned} \mathbf{F}_{21} &= \mathbf{F}_{21}^{(0)} = 0 \\ \mathbf{F}_{22} &= \mathbf{F}_{22}^{(0)} = 1 \\ \mathbf{F}_{11} &= \mathbf{F}_{11}^{(0)} + \int_{0}^{t} \dot{\gamma} \mathbf{F}_{21}(t) dt = 1 \\ \mathbf{F}_{12} &= \mathbf{F}_{12}^{(0)} + \int_{0}^{t} \dot{\gamma} \mathbf{F}_{22}(t) dt = \dot{\gamma} t \end{aligned} \right\} \Rightarrow \mathbf{F}(t) = \begin{pmatrix} 1 & \dot{\gamma} t \\ 0 & 1 \end{pmatrix} \quad (21) $$

and Eq. 18 results:

$$ \left( \begin{array}{c} \mathbf{n}_1 \\ \mathbf{n}_2 \end{array} \right) dS(t) = \left( \begin{array}{c} 1 \\ -\dot{\gamma} t \end{array} \right) dS_0 \\ \Rightarrow \left\{ \begin{array}{l} dS(t) = dS_0 \sqrt{1 + \left(\dot{\gamma} t\right)^2} \\ \mathbf{n}^T = \frac{1}{\sqrt{1+(\dot{\gamma} t)^2}} \frac{-\dot{\gamma} t}{\sqrt{1+(\dot{\gamma} t)^2}} \end{array} \right. \quad (22) $$

which proves that the surface has not a long time limit, i.e.

$$ \left\{ \begin{array}{l} \lim_{t \to \infty} \left( \frac{dS(t)}{dS_0} \right) = \infty \\ \lim_{t \to \infty} (\mathbf{n}) = \begin{pmatrix} 0 \\ -1 \end{pmatrix} \end{array} \right. \quad (23) $$

proving that the elementary surface rotates in the clockwise until aligning in the flow direction, at the same time as it is subjected to an infinite stretching. The effect of the closure relation will be analyzed later.

## Introducing other physical effects

The formalism associated with the evolution of the area tensor has other advantages. One of them is the simplicity for incorporating other physical effects as for example the surface tension. One could expect that in absence of strain the disperse microstructure evolves towards an isotropic state consisting of coalescent microspheres, minimizing the specific area. Thus, Doi and Ohta [5] proposed introducing in Eq. 7 a source terms to take into account such effects:

$$ \frac{d\mathbf{A}}{dt} = -(\mathrm{Grad}\mathbf{v})^{\mathrm{T}}\mathbf{A} - A\mathrm{Grad}\mathbf{v} + \mathrm{Grad}\mathbf{v} : {}^{4}\mathbf{A} \\ -a \frac{\sigma}{\eta} S_{v}^{2} \left( (\tilde{\mathbf{A}} - \frac{I}{3}) + b \frac{I}{3} \right) \quad (24) $$

where $a$ and $b$ are two material parameters, $\sigma$ is the surface tension, $\eta$ the fluid viscosity (we are assuming that both fluids have the same viscosity) and $I$ the unit tensor. Thus, in absence of strain, Eq. 24 reduces to:

$$ \frac{d\mathbf{A}}{dt} = -a \frac{\sigma}{\eta} S_v^2 \left( \left( \tilde{\mathbf{A}} - \frac{I}{3} \right) + b \frac{I}{3} \right) \quad (25) $$

proving that the microstructure evolves towards an isotropic state induced by the presence of the term $\tilde{\mathbf{A}} - \frac{I}{3}$ in which the interface area reduces because the presence of
$-\frac{ab\sigma}{\eta} S_v^2 \frac{I}{3}$. Even if this model predicts a null specific area at infinite time, fact that motivated different corrections [8], we will consider the Doi-Ohta model in the kinetic theory approach described in the next section.

## A kinetic theory model of micromixing

To define a kinetic theory model we start introducing the area distribution function $\Psi(x, t, n)$, given at each point in the physical domain $x \in \Omega$ and for any time $t$, the specific surface oriented in the direction $n$. The area tensor can be obtained from the area distribution function according to:

$$A(x, t) = \int_{S(O,1)} n \otimes n \Psi(x, t, n) dn \quad (26)$$

where $S(O, t)$ represents the surface of the unit sphere centered at the origin. Now, the trace of the previous equation reads:

$$ \begin{aligned} \mathrm{Tr}(A(x,t)) &= S_v(x,t) = \int_{S(O,1)} \mathrm{Tr}(n \otimes n) \Psi(x,t,n) dn \\ &= \int_{S(O,1)} \Psi(x,t,n) dn \end{aligned} \quad (27) $$

The equation governing the evolution of the area tensor defined in Eq. 26 involves the time derivative of the area orientation vector $n$ given by Eq. 5 as well as the time derivative of the distribution function $\Psi$. The expression of the time derivative of that distribution function should allow deducing Eq. 7. It is easy to prove that such time derivative must have the following expression:

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial n}(\dot{n}\Psi) - (\mathrm{Grad}v:(n\otimes n))\Psi \quad (28) $$

that is close to the Fokker-Planck equation usually employed for describing the microstructure in particle suspensions. However in the present case because this function is not subjected to a normality condition ($\Psi$ is not a probability density function) there is the additional last term in the right member of Eq. 28 accounting for the interface area growing induced by its stretching.

If Eqs. 28 and 26 are solved instead the ones related to Eq. 7 and its associated closure relation, the area tensor can be computed without the necessity of introducing any closure relation, whose impact could be in some cases significant.

## Numerical issues

Equation 28 involves a material derivative with the associated advection effects in the physical coordinates:

$$ \frac{d\Psi}{dt} = \frac{\partial\Psi}{\partial t} + v \cdot \nabla_x \Psi \quad (29) $$

where $\nabla_x$ denotes the gradient differential operator in the physical coordinates, i.e.

$$ \nabla_x^T = \left( \frac{\partial}{\partial x_1} \quad \frac{\partial}{\partial x_2} \quad \frac{\partial}{\partial x_3} \right) \quad (30) $$

In flows involving a homogeneous microstructure this material derivative reduces to the partial one. In what follows only the case of flows leading to homogeneous solutions will be considered.

Moreover, in the conformation space related to the orientation coordinates, Eq. 28 also involves advection effects. Thus, for solving this equation an appropriate stabilization technique must be employed, e.g. discontinuous Galerkin or SUPG finite elements among many other possibilities [12]. However, all these stabilization techniques introduce some amount of numerical diffusion altering the computed solution. For this reason, in this work we are considering a particles technique that leads to very accurate solutions when the number of considered particles is large enough.

The particles technique starts considering that the initial area distribution function $\Psi(t=0, n) \equiv \Psi^{(0)}(n)$ can be discretized using N Dirac's masses initially located at $^{(i)}n^{(0)}$ allowing describing the initial microstructure morphology from:

$$ \Psi^{(0)}(n) = \sum_{i=1}^{i=N} \psi_i^{(0)} \delta(n - {}^{(i)}n^{(0)}) \quad (31) $$

with

$$ A^{(0)} = \int_{S(O,1)} (n \otimes n) \Psi^{(0)}(n) dn = \sum_{i=1}^{i=N} \psi_i^{(0)} ({}^{(i)}n^{(0)} \otimes {}^{(i)}n^{(0)}) \quad (32) $$

Now, the location of each particle $^{(i)}n$, and their weights $\psi_i$, can be updated at any time $t_{p+1} = (p+1) \cdot \Delta t$ according to Eq. (28). First of all, we consider the following splitting of Eq. 28:

$$ \left\{ \begin{aligned} \frac{\tilde{\Psi} - \Psi^{(p)}}{\Delta t} &= -\frac{\partial}{\partial n}(\dot{n}\Psi) \\ \frac{\Psi^{(p+1)} - \tilde{\Psi}}{\Delta t} &= -(\mathrm{Grad}v:(n\otimes n))\Psi \end{aligned} \right. \quad (33) $$

where the rotation velocity $\dot{n}$ is given by Eq. 5. The first equation in 33 allows to update the particle orientations:

$$ ^{(i)}n^{(p+1)} = ^{(i)}n^{(p)} + \dot{n}(^{(i)}n^{(p)}) \cdot \Delta t \quad (34) $$

and the second one the particles weights:

$$ \psi_i^{(p+1)} = \psi_i^{(p)} - (\mathrm{Grad}v:(^{(i)}n^{(p)} \otimes ^{(i)}n^{(p)})) \cdot \Delta t \quad (35) $$

that results identical to the Lagrangian strategy considered in Wetzel and Tucker [15] for generating reference data or
the one considered in Mackley and Neves Saraiva [9] for computing the evolution of the interface area.

As we prove later in the “Introducing other physical effects and related numerical issues” subsection, this kinetic theory framework allows an easy introduction of other physical effects avoiding the introduction of closure relations.

## Numerical example

In this section we are considering again the 2D numerical example involving the simple shear flow defined in the “Illustrative example” subsection. We are considering that the initial morphology could be accurately described by:

$$ \Psi^{(0)}(\varphi) = \sum_{i=1}^{i=N} \psi_i^{(0)} \delta(\varphi - \varphi_i^{(0)}) \quad (36) $$

where $\varphi_i$ represents the polar description of the init vector $^{(i)}n$ and

$$ \left\{ \begin{array}{l} \varphi_i^{(0)} = (i - 1) \frac{2\pi}{N} \\ \psi_i^{(0)} = \frac{1}{N} \end{array} \right. \qquad (37) $$

which results in the following second order area tensor

$$ \begin{aligned} A^{(0)} &= \int_{0}^{2\pi} \begin{pmatrix} \cos^2(\varphi) & \cos(\varphi) \cdot \sin(\varphi) \\ \cos(\varphi) \cdot \sin(\varphi) & \sin^2(\varphi) \end{pmatrix} \Psi^{(0)}(\varphi) d\varphi = \\ &= \frac{1}{N} \sum_{i=1}^{i=N} \begin{pmatrix} \cos^2(\varphi_i^{(0)}) & \cos(\varphi_i^{(0)}) \cdot \sin(\varphi_i^{(0)}) \\ \cos(\varphi_i^{(0)}) \cdot \sin(\varphi_i^{(0)}) & \sin^2(\varphi_i^{(0)}) \end{pmatrix} \\ &= \begin{pmatrix} 0.5 & 0 \\ 0 & 0.5 \end{pmatrix} \end{aligned} \quad (38) $$

describing a homogeneous isotropic microstructure with a unit specific surface.

Now, the orientation of particles and their corresponding weights are updated according to Eqs. 34 and 35 allowing to compute the area tensor at any time $t_p$:

$$ A^{(p)} = \frac{1}{N} \sum_{i=1}^{i=N} \begin{pmatrix} \cos^2(\varphi_i^{(p)}) & \cos(\varphi_i^{(p)}) \cdot \sin(\varphi_i^{(p)}) \\ \cos(\varphi_i^{(p)}) \cdot \sin(\varphi_i^{(p)}) & \sin^2(\varphi_i^{(p)}) \end{pmatrix} \quad (39) $$

When this technique is applied to a shear flow characterized by $\dot{\gamma} = 0.1 \text{ s}^{-1}$ both the evolution of the area tensor and the evolution of the specific surface can be computed. Figure 1 compares the solution obtained integrating the evolution equation of the area tensor (Eq. 7) with a quadratic closure relation and the one obtained using the just described particles algorithm. Figure 2 compares

<figure><img src="image_2.png" /><figcaption>Fig. 1 Evolution of the components of the area tensor in a simple shear flow by using both the area tensor and the kinetic theory descriptions</figcaption></figure>

the specific surface evolution obtained by using both techniques. Significant deviations are noticed which as expected goes up as the time increases.

Now, we are applying the just described strategy for solving the problem considered in the “Introducing other physical effects” subsection, using the particles technique, whose solution could be compared with the exact one. If instead of considering the isotropic initial distribution given by Eq. 37 we consider the alternative one in which all the particles are initially oriented in the flow direction:

$$ \left\{ \begin{array}{l} \varphi_i^{(0)} = 0 \\ \psi_i^{(0)} = \frac{1}{N} \end{array} \right. \qquad (40) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2 Evolution of the specific surface in a simple shear flow by using both the area tensor and the kinetic theory descriptions</figcaption></figure>
the solution can be computed with the same accuracy independently of the size of the particles population N. Now, the exact solution is calculated from Eq. 22:

$$ \left\{ \begin{array}{l} dS(t) = dS_0 \sqrt{1 + (\dot{\gamma}t)^2} \\ \mathbf{n}^T = \frac{1}{\sqrt{1+(\dot{\gamma}t)^2}} \quad \frac{-\dot{\gamma}t}{\sqrt{1+(\dot{\gamma}t)^2}} \end{array} \right. \qquad (41) $$

The first conclusion is that when the time is large enough the previous expression can be approximated by;

$$ \left\{ \begin{array}{l} \mathbf{n}^T \approx \begin{pmatrix} 0 & -1 \\ -1 & 0 \end{pmatrix} \\ dS(t \to \infty) \approx \dot{\gamma} \, t dS_0 \end{array} \right. \qquad (42) $$

leading to a linear evolution of the specific surface:

$$ A(t \to \infty) \approx \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix} \dot{\gamma} \, t S_0 \Rightarrow S_v(t \to \infty) \approx \dot{\gamma} \, t S_0 \qquad (43) $$

Now, the exact solution can be computed for any time according to Eq. 41:

$$ S_v(t) = S_0 \sqrt{1 + (\dot{\gamma}t)^2} \qquad (44) $$

Figure 3 compares the evolution of the specific surface computed using the particles strategy and exact solution given by Eq. 44. This comparison confirms the expected excellent accuracy of the particles strategy.

The steady state solution of Eq. 28 can be obtained analytically for the simple shear flow considered in this

<figure><img src="image_4.png" /><figcaption>Fig. 3 Evolution of the specific surface: computed solution using the particles strategy versus the exact one</figcaption></figure>

section. Thus, the evolution equation for the distribution function writes:

$$ \begin{align} \frac{\partial}{\partial \varphi} (\dot{\gamma} \cos^2(\varphi) \Psi) - \dot{\gamma} \sin(\varphi) \cos(\varphi) \Psi &= 0 \nonumber \\ \Rightarrow \cos^2(\varphi) \frac{\partial \Psi}{\partial \varphi} - 3 \sin(\varphi) \cos(\varphi) \Psi &= 0 \Rightarrow \nonumber \\ \Rightarrow \ln(\Psi) &= C - 3 \ln(\cos(\varphi)) \Rightarrow \Psi = \frac{\tilde{C}}{\cos^3(\varphi)} \tag{45} \end{align} $$

The integral of that function diverges because in the neighborhood of $\varphi = \frac{\pi}{2}$ the cosine can be approximated by the angle $\varphi - \frac{\pi}{2} \equiv \tilde{\varphi}$, i.e. $\cos(\varphi) \approx \varphi - \frac{\pi}{2} \equiv \tilde{\varphi}$ for $\varphi \approx \frac{\pi}{2}$, and the integral of $\frac{1}{\tilde{\varphi}^2}$ in any interval containing $\tilde{\varphi} = 0$ diverges if $q \ge 1$, as it is the case here. This result proves that the area does not reach a finite steady value.

Introducing other physical effects and related numerical issues

After proving the ability of the kinetic theory approach for solving passive mixing models, and the accuracy of the particle strategy for solving those models, we are proving that other physical effects, and in particular the one related to the surface tension, can be easily introduced.

As previously, we are considering the expression of the time derivative of the area tensor (that in the kinetic theory framework is given by Eq. 26), looking for the time derivative of the distribution function $\Psi$ leading to the expression (24).

It is easy to prove that the associated kinetic model is given by:

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial \mathbf{n}} (\dot{\mathbf{n}}\Psi) - (\mathrm{Grad}\mathbf{v}:(\mathbf{n}\otimes\mathbf{n}))\Psi + \frac{\partial}{\partial \mathbf{n}}\left(D\frac{\partial\Psi}{\partial \mathbf{n}}\right) + \mathrm{F} \qquad (46) $$

where the diffusion coefficient and the source terms result:

$$ D = \frac{a}{4} \frac{\sigma}{\eta} S_v = \frac{a}{4} \frac{\sigma}{\eta} \left( \int_{S(O,1)} \Psi(\mathbf{n}) d\mathbf{n} \right) \qquad (47) $$

and

$$ F = -\frac{a}{4\pi} \frac{\sigma}{\eta} (S_v)^2 b = -\frac{a}{4\pi} \frac{\sigma}{\eta} \left( \int_{S(O,1)} \Psi(\mathbf{n}) d\mathbf{n} \right)^2 b \qquad (48) $$

Equation 46 is an advection-diffusion equation which can be accurately solved using numerous stabilized discretization strategies. In the case of complex flows this equation is defined in a high-dimensional space involving physical and conformation coordinates (x and n respectively) and time. Despite this apparent difficulty several strategies have been
recently proposed to take into account efficiently the high dimensionality and the extremely large discrete systems resulting from the discretization [1, 2, 4, 7, 10].

When surface tension exists the impact of the closure relation could be more significant for some classes of rheometric flows. Thus, we are considering again a 2D model (without any practical interest but that allows illustrating the results with more clarity) in a small oscillation shear flow, whose kinematics is given by:

$$ \dot{\gamma}(t) = \dot{\gamma}_0 \cos(\omega t) \quad (49) $$

with $\dot{\gamma}_0 = 0.1s^{-1}$ and $\omega=0.1s^{-1}$. The material parameters are chosen (arbitrarily) as $\frac{a\sigma}{\eta} = 0.1$ and $b=0.1$. The initial condition is given by an isotropic orientation distribution that reads:

$$ A^{(0)} = \begin{pmatrix} 0.5 & 0 \\ 0 & 0.5 \end{pmatrix} \quad (50) $$

The solution of the transient Fokker-Planck equation (46) leads to the time evolution of the area tensor $A(t)$:

$$ A(t) = \begin{pmatrix} A_{11}(t) & A_{12}(t) \\ A_{12}(t) & A_{22}(t) \end{pmatrix} \quad (51) $$

The first numerical experiment is carried out in absence of flow, because in this case as the closure relation does not work the results computed using the area tensor evolution equation and the one computed within the kinetic theory framework should be the same. Figure 4 depicts the evolution of the different components of the area tensor: solid lines represents the solution obtained solving the evolution equation of the area tensor and broken lines the solution obtained from the solution of the Fokker-Planck equation. As expected, both solutions are superposed. The

green curve represents the $A_{11}(t)$ and $A_{22}(t)$ components of the area tensor computed using the both just referred strategies. The blue curve represents both $A_{12}(t)$ solutions. As expected, the isotropy is maintained whereas the specific surface $S_v(t) = A_{11}(t) + A_{22}(t)$ is decreasing because the bubbles coalescence.

Now, starting from the same initial conditions, isotropic orientation distribution and unit specific surface, the small amplitude shear oscillation flow is applied. Figure 5 illustrates the evolution of the different area tensor components computed using both strategies: the one based on the evolution equation of the area tensor (which involves the quadratic closure relation) and the one computed within the kinetic theory framework which does not require any closure relation. As expected the gap between both solutions exist and goes up as the time increases.

## Conclusions

In this paper we proposed a new framework for analyzing mixing processes. In this approach the morphology at the microscopic scale is described from the area tensor which contains a lot of information: specific interface, characteristic length, shape and orientation, ... However classical descriptions based on the use of this area tensor fail because the necessity of introducing a closure relation whose impact could be significant and in general unpredictable. In the approach here introduced the area tensor has been reformulated from a distribution function, whose evolution is governed by an advection equation, which allows solving the problem without the necessity of introducing any closure relation. Moreover, this approach allows the incorporation of other physical effects as illustrated in the last section with the

<figure><img src="image_5.png" /><figcaption>Fig. 4 Evolution of the area tensor in absence of flow</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5 Evolution of the area tensor in a small amplitude shear oscillation flow</figcaption></figure>
introduction of surface tension effects. Efficient solvers for
the high-dimensional resulting advection–diffusion equation
based on model reduction or separated representations are
under progress.

References

1. Ammar A, Ryckelynck D, Chinesta F, Keunings R (2006) On the reduction of kinetic theory models related to finitely extensible dumbbells. J Non-Newton Fluid Mech 134:136–147

2. Ammar A, Mokdad B, Chinesta F, Keunings R (2006) A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. J Non-Newton Fluid Mech 139:153–176

3. Chella R, Ottino JM (1985) Stretching in some classes of fluid motions and asymptotic mixing efficiencies as a measure of flow classification. Arch Ration Mech Anal 90:15–42

4. Chinesta F, Chaidron G, Poitou A (2003) On the solution of Fokker–Planck equations in steady recirculating flows involving short fiber suspensions. J Non-Newton Fluid Mech 113:97–125

5. Doi M, Ohta T (1991) Dynamics and rheology of complex interfaces. J Chem Phys 95:1242–1248

6. Hall KR, Godfrey JC (1968) The mixing rates of highly viscous Newtonian and non-Newtonian fluids in a laboratory sigma-blade mixer. Trans Inst Chem Eng 46:205–212

7. Keunings R (2004) Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory. In: Binding DM, Walters K (eds) Rheology reviews. British Society of Rheology, UK, pp 67–98

8. Lee HM, Park OO (1994) Rheology and dynamics of immiscible polymer blends. J Rheol 38:1405–1425

9. Mackley MR, Neves Saraiva RMC (1999) The quantitative description of fluid mixing using Lagrangian- and concentration-based numerical approaches. Chem Eng Sci 54:159–170

10. Öttinger HC, Laso M (1992) Smart polymers in finite element calculation. Int Congress on Rheology, Brussels, Belgium

11. Ottino JM (1989) The kinematics of mixing: stretching, chaos and transport. Cambridge University Press, UK

12. Pironneau O (1989) Finite methods for fluids Wiley, New York

13. Rielly CD, Smith DLO, Lindley JA, Niranjan K, Philips VR (1994) Mixing processes for agricultural and food materials: Part 4, assessment and monitoring of mixing systems. J Agric Eng Res 29:1–18

14. Roberts EPR, Mackley MR (1995) The simulation of stretch rates for the quantitative prediction and mapping of mixing within a channel flow. Chem Eng Sci 50:3727–3746

15. Wetzel ED, Tucker CL III (1999) Area tensors for modeling microstructure during laminar liquid–liquid mixing. Int J Multiph Flow 25:35–61