# Chaotic mixing

An example of chaotic mixing

In chaos theory and fluid dynamics, chaotic mixing is a process by which flow tracers develop into complex fractals under the action of a fluid flow. The flow is characterized by an exponential growth of fluid filaments.[1] Even very simple flows, such as the blinking vortex, or finitely resolved wind fields can generate exceptionally complex patterns from initially simple tracer fields.[2]

The phenomenon is still not well understood and is the subject of much current research.

### Fluid flows

Two basic mechanisms are responsible for fluid mixing: diffusion and advection. In liquids, molecular diffusion alone is hardly efficient for mixing. Advection, that is the transport of matter by a flow, is required for better mixing.

The fluid flow obeys fundamental equations of fluid dynamics (such as the conservation of mass and the conservation of momentum) called Navier–Stokes equations. These equations are written for the Eulerian velocity field rather than for the Lagrangian position of fluid particles. Lagrangian trajectories are then obtained by integrating the flow. Studying the effect of advection on ﬂuid mixing amounts to describing how different Lagrangian ﬂuid particles explore the ﬂuid domain and separate from each other.

A fluid flow can be considered as a dynamical system, that is a set of ordinary differential equations that determines the evolution of a Lagrangian trajectory. These equations are called advection equations:

$\frac{d \vec x}{dt} = \vec v(\vec x, t)$

where $\vec v=(u, v, w)$ are the components of the velocity field, which are assumed to be known from the solution of the equations governing fluid flow, such as the Navier-Stokes equations, and $\vec x=(x, y, z)$ is the physical position. If the dynamical system governing trajectories is chaotic, the integration of a trajectory is extremely sensitive to initial conditions, and neighboring points separate exponentially with time. This phenomenon is called chaotic advection.

2-D Stationary flow: The velocity v is tangent to streamlines (dotted lines). For a stationary flow (left), Lagrangian trajectories X coincide with streamlines. For an unstationary flow (right), the position of streamlines changes with time and trajectories do no coincide with streamlines, hence better mixing is possible.

Dynamical systems and chaos theory state that at least 3 degrees of freedom are necessary for a dynamic system to be chaotic. Three-dimensional flows have three degrees of freedom corresponding to the three coordinates, and usually result in chaotic advection, except when the flow has symmetries that reduce the number of degrees of freedom. In flows with less than 3 degrees of freedom, Lagrangian trajectories are confined to closed tubes, and shear-induced mixing can only proceed within these tubes.

This is the case for 2-D stationary flows in which there are only two degrees of freedom $x$ and $y$. For stationary (time-independent) flows, Lagrangian trajectories of fluid particles coincide with the streamlines of the flow, that are isolines of the stream function. In 2-D, streamlines are concentric closed curves that cross only at stagnation points. Thus, a spot of dyed fluid to be mixed can only explore the region bounded by the most external and internal streamline, on which it is lying at the initial time. Regarding practical applications, this configuration is not very satisfying.

Blinking vortex flow: (a) The flow is created by the periodic alternated rotation of rods (b) The trajectory consists of arcs of circles. The particle can wander through the fluid domain (c) A small black disk of fluid is stretched into a large number that occupy the whole width of the fluid domain.

For 2-D unstationary (time-dependent) flows, instantaneous closed streamlines and Lagrangian trajectories do not coincide any more. Hence, Lagrangian trajectories explore a larger volume of the volume, resulting in better mixing. Chaotic advection is observed for most 2-D unstationary flows. A famous example is the blinking vortex flow introduced by Aref,[3] where two fixed rod-like agitators are alternately rotated inside the fluid. Switching periodically the active (rotating) agitator introduces a time-dependency in the flow, which results in chaotic advection. Lagrangian trajectories can therefore escape from closed streamlines, and visit a large fraction of the fluid domain.

### Shear

A flow promotes mixing by separating neighboring fluid particles. This separation occurs because of velocity gradients, a phenomenon called shearing. Let $\vec{x}_1$ and $\vec{x}_2$ be two neighboring fluid particles, separated by $\delta\vec{x}=\vec{x}_2-\vec{x}_1$ at time t. When the particles are advected by a flow $\vec{v}$, at time $t+\delta t$ the approximate separation between the particles can be found through Taylor expansion :

$\frac {\mathrm d} {\mathrm dt} (\vec{x}+\delta \vec{x}) \approx \vec{v} + \nabla \vec{v} \cdot \delta \vec{x}$

hence

$\delta x(t+\delta t) \approx \delta x(t) + \delta t (\delta x \cdot \nabla) \vec v$

and

$\frac {\mathrm d} {\mathrm dt} \delta \vec x \approx \nabla \vec v \cdot \delta \vec x$

The rate of growth of the separation is therefore given by the gradient of the velocity field in the direction of the separation. The plane shear ﬂow is a simple example of large-scale stationary flow that deforms fluid elements because of a uniform shear.

### Lyapunov exponents

If the flow is chaotic, then small initial errors, $\delta \vec x$, in a trajectory will diverge exponentially. We are interested in calculating the stability—i.e., how fast do nearby trajectories diverge? The Jacobi matrix of the velocity field, $\nabla \vec{v}$, provides information about the local rate of divergence of nearby trajectories or the local rate of stretching of Lagrangian space.

We define the matrix H such that:

$\frac {\mathrm d} {\mathrm dt} \boldsymbol{H} \equiv \nabla \vec{v} \cdot \boldsymbol{H}, \qquad \boldsymbol{H} (t=0)=\boldsymbol{I}$

where I is the identity matrix. It follows that:

$\delta \vec{x} (t) \approx \boldsymbol{H} \cdot \delta \vec{x}_0$
(a) In chaotic advection, neighboring trajectories separate exponentially fast with time. (b) Numerical integration of two trajectories initially very close, for a 2-D Stokes flow where a rod stirs fluid by moving on a figure-eight path. The trajectories separate fast because of chaotic advection, until they are completely uncorrelated.

The finite-time Lyapunov exponents are defined as the time average of the logarithms of the lengths of the principal components of the vector H over a time t:

$\boldsymbol{H^T} \cdot \boldsymbol{H} \cdot \delta \vec{x}_{0i} = h_i \cdot \delta \vec{x}_{0i}$
$\lambda_i(\vec{x},t) \equiv \frac {1} {2 t} \ln {h_i(\vec{x},t)}$

where $\lambda_i(\vec{x},t) \geq \lambda_{i+1}(\vec{x},t)$ is the ith Lyapunov exponent of the system, while $\delta \vec {x}_{0i}$ is the ith principal component of the matrix H.

If we start with a set of orthonormal initial error vectors, $\{\delta \vec x_{0i}\}$ then the matrix H will map them to a set of final orthogonal error vectors of length $\{\sqrt{h_i(\vec{x},t)}\}$. The action of the system maps an infinitesimal sphere of inititial points to an ellipsoid whose major axis is given by the $\sqrt{h_1(\vec{x},t)}$ while the minor axis is given by $\sqrt{h_N(\vec{x},t)}$, where N is the number of dimensions.[4][5]

This definition of Lyapunov exponents is both more elegant and more appropriate to real-world, continuous-time dynamical systems than the more usual definition based on discrete function maps. Chaos is defined as the existence of at least one positive Lyapunov exponent.

In a chaotic system, we call the Lyapunov exponent the asymptotic value of the greatest eigenvalue of H:

$\lambda = \lim_{t \to \infty} \lambda_1(\vec{x},t)$

If there is any significant difference between the Lyapunov exponents then as an error vector evolves forward in time, any displacement in the direction of largest growth will tend to be magnified. Thus:

$|\delta \vec x| \approx |\delta \vec x_0| e^{\lambda_1 t}.$

The Lyapunov exponent of a flow is a unique quantity, that characterizes the asymptotic separation of fluid particles in a given flow. It is often used as a measure of the efficiency of mixing, since it measures how fast trajectories separate from each other because of chaotic advection. The Lyapunov exponent can be computed by different methods:

• by following one single trajectory for very long times and computing $\lambda = \lim_{t \to \infty} \lambda_1(\vec{x},t)$.
• or by following an ensemble of trajectories for a given period of time, and computing the ensemble average: $<\lambda>_{trajectories}$

The equivalence of the two methods is due to the ergodicity of the chaotic system.

### Filament growth versus evolution of the tracer gradient

The following, exact equation can be derived from an advection-diffusion equation (see below), with a diffusion term (D=0) of zero:

$\frac{\mathrm d \nabla q}{\mathrm d t} = -\nabla q \cdot \nabla \vec v$

In parallel with the definition of the Lyapunov exponent, we define the matrix $\boldsymbol{H}^\prime$, as follows:

$\frac{\mathrm d \boldsymbol{H^\prime}}{\mathrm d t} = -\nabla \boldsymbol{H^\prime} \cdot \nabla \vec v \qquad \boldsymbol{H^\prime} (t=0)=\boldsymbol{I}$

It is easy to show that:

$\boldsymbol{H^\prime}=\boldsymbol{H}^{-1}$

If we define $\lbrace h_i^\prime\rbrace$ as the squared lengths of the principal components of the tracer gradient matrix, $\boldsymbol{H}^\prime$, then:

$h_i^\prime=1/h_i$

where the $\lbrace h_i^\prime\rbrace$'s are arranged, as before, from largest to smallest. Therefore, growth in the error vector will cause a corresponding decrease in the tracer gradient and vice versa. This can be understood very simply and intuitively by considering two nearby points: since the difference in tracer concentration will be fixed, the only source of variation in the gradients between them will be their separation.[4][6]

Contour advection is another useful method for characterizing chaotic mixing. In chaotic flows, advected contours will grow exponentially over time. The figure above shows the frame-by-frame evolution of a contour advected over several days. The figure to the right shows the length of this contour as a function of time.

The link between exponential contour growth and positive Lyapunov exponents is easy to see. The rate of contour growth is given as:

$\frac{\mathrm d L}{\mathrm d t} = \int | \nabla \vec v \cdot \mathrm d \vec s |$

where $\mathrm d \vec s$ is the path and the integral is performed over the length of the contour. Contour growth rates will approximate the average of the large Lyapunov exponents:[4]

$L \approx L_0 \exp(\bar \lambda_1 t)$

### Poincaré sections

In chaotic advection, a fluid particle travels within a large region, and encounters other particles that were initially far from it. One can then consider that a particle is mixed with particles that travel within the same region. However, the region covered by a trajectory does not always span the whole fluid domain. Poincaré sections are used to distinguish regions of good and bad mixing.

The Poincaré map is defined as the transformation

\begin{align} \boldsymbol{M} \colon \vec{x_i}(t_i)&\to \vec{x}_{i+1}(t_{i+1}=t_i+T,\vec{x_i}). \end{align}

$\boldsymbol{M}$ transforms a point-like particle into the position of the particle after a time-interval T. Especially, for a time-periodic flow with period T, applying the map several times to a particle gives the successive positions of the particle period after period. A Poincaré section is built by starting from a few different initial conditions and plotting the corresponding iterates. This comes down to plotting the trajectories stroboscoped every T.

Poincaré sections for two versions of the figure-eight protocol. The chaotic spans the entire domain for one protocol (right), while elliptic islands are visible for the other protocol (left).

Let us consider an example. The figure presented here (left part) depicts the Poincaré section obtained when one applies periodically a figure-eight-like movement to a circular mixing rod. Some trajectories span a large region: this is the chaotic or mixing region, where good mixing occurs. However, there are also two "holes": in these regions, the trajectories are closed. These are called elliptic islands, as the trajectories inside are elliptic-like curves. These regions are not mixed with the remainder of the fluid. For mixing applications, elliptic islands have to be avoided for two reasons :

• Fluid particles are unable to cross the boundaries of the islands (except by slow diffusion), resulting in segregation.
• Mixing inside these regions is not efficient because trajectories are closed and therefore not chaotic.

Avoiding non-chaotic islands requires understanding the physical origin of these regions. Generally speaking, changing the geometry of the flow can modify the presence or absence of islands. In the figure-eight flow for instance, for a very thin rod, the influence of the rod is not felt far from its location, and almost circular trajectories exist within the loops of the figure-eight. With a larger rod (right part of the figure), particles can escape from these loops and islands do not exist any more, resulting in better mixing.

With a Poincaré section, the mixing quality of a flow can be analyzed by distinguishing between chaotic and elliptic regions. This is a crude measure of the mixing process, however, since the stretching properties cannot be inferred from this mapping method. Nevertheless, this technique is very useful for studying the mixing of periodic flows and can be extended to a 3-D domain.

### Fractal dimension

Through a continual process of stretching and folding, much like in a "baker's map," tracers advected in chaotic flows will develop into complex fractals. The fractal dimension of a single contour will be between 1 and 2. Exponential growth ensures that the contour, in the limit of very long time integration, becomes fractal. Fractals composed of a single curve are infinitely long and when formed iteratively, have an exponential growth rate, just like an advected contour. The Koch Snowflake, for instance, grows at a rate of 4/3 per iteration.

The figure below shows the fractal dimension of an advected contour as a function of time, measured in four different ways. A good method of measuring the fractal dimension of an advected contour is the uncertainty exponent.

### Evolution of tracer concentration fields in chaotic advection

In fluid mixing, one often wishes to homogenize a species, that can be characterized by its concentration field q. Often, the species can be considered as a passive tracer that does not modify the flow. The species can be for example a dye to be mixed. The evolution of a concentration field $c$ obeys the advection-diffusion equation, also called Convection–diffusion equation:

$\big. \frac{\partial q}{\partial t} = \nabla \cdot D \nabla q - \vec{v} \cdot \nabla q.$

Compared to the simple diffusion equation, the term proportional to the velocity field $\vec{v}$ represents the effect of advection.

When mixing a spot of tracer, the advection term dominates the evolution of the concentration field at the beginning of the mixing process. Chaotic advection transforms the spot into a bundle of thin filaments. The width of a dye filament decreases exponentially with time, until an equilibrium scale is reached, at which the effect of diffusion starts to be significant. This scale is called the Batchelor scale. It is defined as the square root of the ratio between the diffusion coefficient and the Lyapunov exponent

$w_B = \sqrt{\frac{D}{\lambda}}$

where $\lambda$ is the Lyapunov exponent and D is the diffusion coefficient. This scale measures the balance between stretching and diffusion on the evolution of the concentration field: stretching tends to decrease the width of a filament, while diffusion tends to increase it. The Batchelor scale is the smallest lengthscale that can be observed in the concentration field, since diffusion smears out quickly any finer detail.

When most dye filaments reach the Batchelor scale, diffusion begins to decrease significantly the contrast of concentration between the filament and the surrounding domain. The time at which a filament reaches the Batchelor scale is therefore called its mixing time. The resolution of the advection–diffusion equation shows that after the mixing time of a filament, the decrease of the concentration fluctuation due to diffusion is exponential, resulting in fast homogenization with the surrounding fluid.