# Mixing by stirring:
optimizing shapes and strategies

Maximilian F. Eggl<EMAIL_ADDRESS>Institute of Physiological
Chemistry, Johannes Gutenberg University of Mainz, Mainz, Germany  Peter J.
Schmid Department of Mechanical Engineering, KAUST, Thuwal, Makkah, Saudi
Arabia

###### Abstract

The mixing of binary fluids by stirrers is a commonplace procedure in many
industrial and natural settings, and mixing efficiency directly translates
into more homogeneous final products, more enriched compounds, and often
substantial economic savings in energy and input ingredients. Enhancements in
mixing efficiency can be accomplished by unorthodox stirring protocols as well
as modified stirrer shapes that utilize unsteady hydrodynamics and vortex-
shedding features to instigate the formation of fluid filaments which
ultimately succumb to diffusion and produce a homogeneous mixture. We propose
a PDE-constrained optimization approach to address the problem of mixing
enhancement for binary fluids. Within a gradient-based framework, we target
the stirring strategy as well as the cross-sectional shape of the stirrers to
achieve improved mixedness over a given time horizon and within a prescribed
energy budget. The optimization produces a significant enhancement in
homogeneity in the initially separated fluids, suggesting promising
modifications to traditional stirring protocols.

††preprint: APS/123-QED

##  I Introduction

The mixing of binary fluids is a fundamental process that underpins a wide
array of applications in fluid mechanics. From tracking ocean pollution [1]
where wind-driven mixing affects the distribution of buoyant plastic debris,
the study of the ventilation of buildings [2, 3] to improved mixing of fuel in
combustion engines [4], the importance of mixing – and the incentive to do it
well – can hardly be overstated. Particularly in industrial settings, mixing
has become increasingly important as most chemical engineering processes are
constantly optimized for energy efficiency, yield, output quality or time-to-
product. Recent examples of this effort can be seen in [5] using computational
tools to study industrial pre-mixers, in [6] investigating helical mixers and
in [7] modeling and optimizing gas stirred ladle systems. Far from exhaustive,
these examples illustrate the substantial interest and research directed
towards a better understanding of mixing and its optimization.

In particular the physical mechanisms responsible for good or acceptable
mixing have been the focus of a great many studies, in an effort to enhance
mixing efficiency and homogeneity of the final mixture. In an early
investigation [8] identified stretching, cutting and folding as the primary
building blocks of mixing, while [9] pioneered a theoretical framework for a
simplified analysis of viscous mixing by ignoring interfacial tension. In a
more recent effort, the fundamental process and the definition of mixing has
attracted the attention of mathematicians [10, 11, 12] who derived necessary
conditions for optimality in the viscous regime and proposed appropriate
measures of mixedness.

While earlier studies have concentrated on the high-viscosity regime, a
significant effort has also been directed towards understanding mixing in an
inertially dominated, chaotic and turbulent setting. Explorations into this
regime include work using direct numerical simulations to study passive scalar
mixing via shock-wave interactions over a wide array of governing parameters
such as Mach numbers, Reynolds number and Schmidt numbers [13]. Based on the
early landmark work of [14] who isolated the process of filament formation
through the stirrers’ shear action which subsequently fed into an effective
diffusion process, numerous studies suggested schemes and mechanistic
procedures to enhance mixing. A modern account of mixing by vortical motion
using a network formalism of interacting vortices has been given by [15].
Furthermore, laminar and turbulent fluid motion has been used and analyzed in
[16], where oscillatory wall motion has been used to enhance mixing, and, more
recently, [17] found that turbulent flow in the near-wall region accomplishes
a markedly higher degree of mixing in a passive scalar field than than either
chaotic or laminar flow.

A better understanding of the underlying mechanisms for effective mixing have
further encouraged investigations into optimizing, or at least enhancing, the
overall mixing process within prescribed bounds and constraints. Treatises on
this topic include the use of boundary feedback control based on open-loop
system characteristics [18, 19] to optimize the velocity profile imposed by a
double-periodic sine-mixer. The newly proposed mixing measures (for example,
the $H^{-1}$-norm, equivalent to the square-root of the variance of a low-pass
filtered passive tracer field) have been used in [20] to achieve improved
mixing. These studies corroborate earlier work by [21] who proved the
existence of an optimal flow and derived associated optimality conditions for
achieving it. In a further study of mixing enhancement, and building on
fundamental preparatory work [22, 23], [24] demonstrated the use of wall-
mounted blowing-and-suction strategies in a channel, coupled to an adjoint-
based optimization formalism, as an effective tool to enhance mixing – if an
appropriate mixing measure [see 10] is considered as an cost objective. This
study has been extended to three-dimensional and stratified flows by [25] and
[26] with similar findings. In a recent article [27], mixing by jets has been
demonstrated in an experimental setup, using machine learning to arrive at
enhanced mixing. Optimization of mixing via jet actuators has also been
pursued in [28], finding that a combined axial and helical actuation is far
more effective of reaching a homogeneous mixture than manipulating solely by
helical actuation. Finally, and similar in approach to this study, [29]
employed a gradient-based approach to enhance passive mixing by optimizing
static baffles.

In this article, we will pursue stirrer-based mixing enhancement by combining
a gradient-based and PDE-constrained optimization formulation with a modern
description of mixing measures. This approach continues earlier work [30, 31]
by extending the optimization space to (i) topology changes for the stirrers’
cross-sections and to (ii) the temporal velocity protocol along given paths.
In addition, we explore a parameter regime that has received less attention
from the mixing community: the inertial, but laminar regime, described by a
Reynolds number above the Stokes-regime, but below a critical value that would
induce turbulent fluid motion. In this chosen regime, we expect and observe a
rich flow behavior, characterized by strong shear and robust vortical motion,
ideal for seeking a sophisticated and perhaps non-intuitive mixing strategy,
while at the same time avoiding exceeding sensitivity due to turbulent
fluctuations (which would run counter to our gradient-based optimization
approach). We anticipate an optimized stirring protocol that, simultaneously
and/or progressively, takes advantage of advective, unsteady and diffusive
processes. Despite the restriction, this parameter regime covers many
applications of industrial and technological interest.

We will concentrate on an active stirrer system and utilize the computational
framework introduced in previous work to optimize the velocity of the stirrers
along a prescribed path as well as their cross-sectional area to obtain mixing
protocols that enhance mixing in a binary, incompressible and two-dimensional
fluid system. We enforce constraints driven by mechanical and material
limitations of the stirrer system, which translate into bounds on the
stirrers’ velocities and accelerations; moreover, we restrict the stirrers’
shapes to avoid multiply-connected cross-sections or delicate shapes. The
focus of this work is in the enhancement of mixing using a fixed time horizon
and energy budget and the description of the identified stirring protocol in
terms of a dynamic interplay of the vortical structures with themselves and
with the embedded stirrers.

The rest of this article is organized as follows; in §II we present the
governing equations, motivate our chosen measure of mixing, and then describe
mixing enhancement as an optimization problem. Following on from this, in §III
we briefly discuss the numerical implementation that arises from our governing
equations; this section is included for completeness, and special emphasis is
drawn on computational concepts arising in the present study. Subsequently in
§IV we introduce three distinct test cases, present the optimization results
and discuss the key features responsible for an enhancement in mixing. The
three optimization experiments have been designed as follows: case _(i)_ aims
to optimize the cross-sectional shape of the stirrers while keeping the
velocities of the stirrers on their respective paths fixed to a prescribed
protocol. The resulting enhanced shapes are then fixed for case _(ii)_ , which
then targets the velocities along the paths in an effort to further enhance
mixing. Lastly, case _(iii)_ presents an optimization where both the cross-
sectional shapes and path velocities are optimized simultaneously, rather than
successively. In §V we draw conclusions and provide an outlook for future
optimization studies.

##  II Theoretical background

###  II.1 Governing equations

The focus of this study is the mixing process of a binary, miscible and
Newtonian fluid by multiple stirrers on prescribed paths, and its optimization
by manipulating the cross-sectional shape of the stirrers and/or their
stirring strategy within specified constraints. A two-dimensional
configuration is considered. The process can be simulated by solving the fluid
equations of motion, augmented by a transport equation for a passive scalar
$\theta$ that identifies the binary fluid. We have

$\displaystyle\frac{\partial}{\partial
t}{\bf{u}}+{\bf{u}}\cdot\nabla{\bf{u}}+\frac{1}{C_{\eta}}\left(\chi{\bf{u}}-\chi_{k}{\bf{u}}_{s,k}\right)+\nabla
p-\frac{1}{Re}\nabla^{2}{\bf{u}}$ $\displaystyle=0,$ (1)
$\displaystyle\nabla\cdot{\bf{u}}$ $\displaystyle=0,$ (2)
$\displaystyle\frac{\partial}{\partial
t}\theta+\left((1-\chi){\bf{u}}+\chi_{k}{\bf{u}}_{s,k}\right)\cdot\nabla\theta-\nabla\cdot\left(\left[\frac{1}{Pe}(1-\chi)+\frac{\chi}{C_{\eta}}\right]\nabla\theta\right)$
$\displaystyle=0.$ (3)

with ${\bf{u}}$ as the velocity vector, $p$ as the pressure field and $\theta$
as a passive scalar (ranging from zero in one fluid to one in the other). The
governing equations have been cast in non-dimensional form using a
characteristic length $L_{0}$ and a velocity scale $u_{0}.$ In our case, we
choose the stirrer diameter and the initial control velocity of the stirrer as
our characteristic quantities. This choice introduces the Reynolds number $Re$
and the Péclet number $Pe,$ which express the kinematic viscosity and the
diffusion coefficient of the mixing fluid in non-dimensional form.
Furthermore, the system (1)-(3) contains terms that model the embedded
stirrers via a Brinkman penalization approach [see 32]. The multiple solid
stirrers, indexed by the subscript $k$, are described in shape by masks
$\chi_{k},$ equal to one for points occupied by the $k$-th stirrer and zero
outside of it. The mask $\chi$ (without subscripts) accounts for the overall
geometry, such as the domain boundaries. In our configuration, the stirrers
travel on circular paths throughout our domain. It thus seems prudent to
redefine the velocity of the $k$-th stirrer $\bm{u}_{s,k}$ using a polar
coordinate system. To this end, we introduce the vector-valued function
$\bm{l}$, which transforms between the polar and Cartesian coordinate systems
according to $\bm{l}(\phi)=(l_{1}(\phi),l_{2}(\phi))=(-\sin(\phi),\cos(\phi))$
with $\phi$ as the angle traversed along the path of the circle. The
parameterization of the velocity of the $k$-th stirrer thus becomes

$\displaystyle\bm{u}_{s,k}$ $\displaystyle=$
$\displaystyle\omega_{k}(t)r_{k}\bm{l}(\varphi_{k}(t)),$ (4)

where $r_{k}$ denotes the radius of the circular path for stirrer $k,$ and
$\varphi$ denotes the sum of the angles traveled along the path according to

$\displaystyle\varphi_{k}(t)=\int_{0}^{t}\omega_{k}({\tau})\ \text{d}{\tau}.$
(5)

The constant $C_{\eta}$ in the governing equations ensures the rapid
relaxation of the fluid variables towards the respective values imposed by the
stirrers or the geometry. The above formalism allows for an efficient
treatment of objects moving through a background grid on which the motion of
the surrounding fluid is described and solved. Details of the above approach,
together with its numerical implementation, can be found in [30]. The
particular set-up, shown in (1)-(3), imposes no-slip velocity boundary
conditions on the stirrers and no-flux Neumann conditions for the passive
scalar on the solid bodies.

The above governing equations provide a constraint on the optimization of
mixedness of the binary fluid via the stirrers. Besides this constraint, a
cost objective has to formulated that quantifies the degree of homogeneity of
the binary fluids. We follow earlier studies and use the mix-norm as a measure
of mixedness. Representing a Sobolev norm of negative fractional index, it is
defined as

$\displaystyle\|\theta\|_{\rm{mix}}$
$\displaystyle\equiv\frac{1}{|V|}\int_{V}\|\nabla^{-2/3}\theta({\bf{x}},t)\|\text{
d}V,$ (6)

with $V$ denoting our computational domain, and $|V|$ representing its size
(volume or area). This type of norm gives weight to larger structures and de-
emphasizes small-scale structures [10, 12, 24]. Used within an optimization
framework, it targets and promotes the breakdown of large scales into smaller
structures, as well as the generation of thin filaments. Smaller values of the
mix-norm tend to indicate a more homogeneous mixture.

To complete the optimization problem we have to decide on a time horizon $T$
over which the stirrers will be active. In addition, we have to specify the
total amount of energy $E_{0}$ exerted by the stirrers, as well as the limits
on the stirrers’ velocities, $\bm{u}_{s,k}$, and accelerations,
$\bm{a}_{s,k}$. For an argument for this set of constraints, the reader is
referred to [33]. The full optimization problem can then be stated as

$\displaystyle\min\left\\{\int_{0}^{T}\|\theta\|_{\rm{mix}}\ dt\right\\}$ (7)
$\displaystyle\hbox{subject to
equations}~{}(\ref{eq:GovEqu1})-(\ref{eq:GovEqu3})$ (8) $\displaystyle\hbox{
and}\quad\int_{0}^{T}\sum_{k}\|\bm{u}_{s,k}\|^{2}\ dt\leq E_{0}$ (9)
$\displaystyle\hbox{
and}\quad\bm{u}_{s,\rm{lower}}\leq\bm{u}_{s,k}\leq\bm{u}_{s,\rm{upper}}\qquad
k=1,2,$ (10) $\displaystyle\hbox{
and}\quad\bm{a}_{s,\rm{lower}}\leq\bm{a}_{s,k}\leq\bm{a}_{s,\rm{upper}}\qquad
k=1,2.$ (11)

This system has to be discretized, embedded within a gradient-based (direct-
adjoint) framework and coupled to a standard unconstrained optimization
routine.

##  III Computational framework

Since the computational framework, including the discretization of the
governing equations, has been presented in full in [30], we will only briefly
discuss the numerical details here for completeness. The starting point for
the computational implementation was the open-source software FLuSI [34, 35],
a parallelized Fourier-based pseudo-spectral code. The three-dimensional
incompressible Navier-Stokes equations are solved on a triple periodic domain
using an equidistant mesh. Time-stepping is performed adaptively by multi-
stage schemes and is based on numerical stability constraints. The interaction
between fluid and solid cells is modeled using a Brinkman-style penalization
method, which was discussed above in relation to the governing equations.

Based on the Fourier representation of the flow field variables, we can
formulate spatial derivatives as a multiplication by a Fourier discretization
matrix. We thus replace the derivative terms in the governing equation as
follows

$\displaystyle\frac{\partial}{\partial x_{i}}\rightarrow\bf{A}_{i}$ (12)

where the index $i$ denotes the respective coordinate direction or the
associated spatial wavenumber. Under this transformation we can state the
semi-discretized equations, where we have assumed the Einstein summation
convention over repeated indices, as

$\displaystyle\frac{d}{dt}{\bf{u}}_{i}+{\bf{u}}_{j}\circ\left[{\bf{A}}_{j}{\bf{u}}_{i}\right]+\frac{1}{C_{\eta}}\left(\chi\circ{\bf{u}}-\chi_{k}\circ{{\bf{u}}_{s,k}}\right)_{i}+{\bf{A}}_{i}p-\frac{1}{Re}{\bf{A}}_{j}{\bf{A}}_{j}{\bf{u}}_{i}$
$\displaystyle=0,$ (13) $\displaystyle{\bf{A}}_{j}{\bf{u}}_{j}$
$\displaystyle=0,$ (14)
$\displaystyle\frac{d}{dt}\theta+\left(({\bf{1}}-\chi)\circ{\bf{u}}+\chi_{k}\circ{{\bf{u}}_{s,k}}\right)_{j}\circ\left[{\bf{A}}_{j}\theta\right]-{\bf{A}}_{j}\left[\frac{1}{Pe}\left({\bf{1}}-\chi\right)+\kappa\chi\right]\circ{\bf{A}}_{j}\theta$
$\displaystyle=0.\phantom{123456}$ (15)

In the above expressions we denote by $\circ$ the Hadamard (elementwise)
product. The solution of this set of equations follows a typical operator-
splitting approach, where a pressure-delayed solution is computed, which
subsequently is corrected by a pressure Poisson equation of the form

${\bf{A}}_{j}{\bf{A}}_{j}p+{\bf{A}}_{i}\left(\bm{u}_{j}\circ\left[{\bf{A}}_{j}\bm{u}_{i}\right]\right)+{\bf{A}}_{j}\left[\frac{\chi}{C_{\eta}}\circ\bm{u}-\frac{\chi_{k}}{C_{\eta}}\circ\bm{u}_{s,k}\right]_{j}=0.$
(16)

The iterative progression of this system of equations forms the forward
problem for our mixing setup, resulting in temporal snapshots of the flow
field and the passive scalar.

###  III.1 Optimization framework using adjoint techniques

The above semi-discretized system has to be embedded in an overall
optimization scheme. This is accomplished by formulating an unconstrained
optimization problem consisting of our cost objective ${\cal{J}}$ (based on
the mix-norm, and defined by (7)) and Lagrange-multiplier terms enforcing all
our constraints. The so-called augmented Lagrangian is given as

$\displaystyle\mathcal{L}=\mathcal{J}-$
$\displaystyle\int_{0}^{T}\bm{u}_{i}^{{\dagger},H}M\left\\{\partial_{t}\bm{u}+\bm{u}_{j}\circ[\bm{A}_{j}\bm{u}]+\frac{\chi}{C_{\eta}}\circ\bm{u}-\frac{\chi_{k}}{C_{\eta}}\circ\bm{u}_{s,k}+\mathbf{A}p-Re^{-1}[\bm{A}_{j}\bm{A}_{j}\bm{u}]\right\\}_{i}$
$\displaystyle+p^{{\dagger},H}M\left\\{[\bm{A}_{j}\bm{A}_{j}]p+\bm{A}_{i}(\bm{u}_{j}\circ[\bm{A}_{j}\bm{u}_{i}])+\bm{A}_{j}\left[\frac{\chi}{C_{\eta}}\circ\bm{u}-\frac{\chi_{k}}{C_{\eta}}\circ\bm{u_{s,k}}\right]_{j}\right\\}$
$\displaystyle+\theta^{{\dagger},H}M\left\\{\partial_{t}\theta+(\bm{1}-\chi)\circ\bm{u}_{j}\circ[\bm{A}_{j}\theta]-\chi\circ(\bm{u}_{s,k})_{j}\circ[\bm{A}_{j}\theta]-\bm{A}_{j}([Pe^{-1}(\bm{1}-\chi)+\kappa\chi]\circ\bm{A}_{j}\theta)\right\\}$
$\displaystyle+\chi_{k}^{{\dagger},H}M[\chi_{k}-g_{k}(\bm{x},t)]+\omega_{k}^{{\dagger},H}M[\omega_{k}-z_{k}(t)]\
\text{d}t$ (17)

where we recognize the various terms of the governing equations, each
premultiplied by a Lagrange multiplier or an adjoint variable. We denote by
$z_{k},g_{k}$ the given previous values of $\omega_{k}$ and $\chi_{k}$,
respectively. All adjoint variables act as independent variables in the
augmented Lagrangian ${\cal{L}},$ and can be linked to sensitivities of the
cost functional ${\cal{J}}$ with respect to respective flow variables or
parameters. The scalar product chosen in the above expression is defined by
the positive-definite weight matrix $M$ which accounts for mesh geometry or
user-specified weights. The next step consists of taking the first variation
of ${\cal{L}}$ with respect to all independent variables and rendering these
variations zero (for an optimum). The first variation with respect to the
adjoint variables recovers the governing equations stated above. The first
variation with respect to the direct variables (i.e., flow velocities,
pressure and passive scalar) produces – after some cumbersome algebraic
manipulations – governing equations for the adjoint variables. They read

$\displaystyle\frac{d}{dt}{\bf{u}}^{{\dagger}}_{i}-\Pi^{{\dagger}}_{j}\circ[{\bf{A}}_{i}{\bf{u}}_{j}]-{\bf{A}}_{j}^{H}[{\bf{u}}_{j}\circ\Pi^{{\dagger}}_{i}]-\frac{\chi}{C_{\eta}}\circ\Pi^{{\dagger}}_{i}+\frac{1}{Re}{\bf{A}}_{j}^{H}{\bf{A}}_{j}^{H}{\bf{u}}_{i}^{{\dagger}}-({\bf{1}}-\chi)\circ\theta^{{\dagger}}\circ[{\bf{A}}_{i}\theta]$
$\displaystyle=0,$ (18) $\displaystyle{\bf{A}}_{j}^{H}\Pi_{j}^{\dagger}$
$\displaystyle=0$ (19)
$\displaystyle\frac{d}{dt}\theta^{{\dagger}}-{\bf{A}}_{j}^{H}[({\bf{1}}-\chi)\circ{\bf{u}}_{j}\circ\theta^{{\dagger}}]+{\bf{A}}_{j}^{H}([\frac{1}{Pe}({\bf{1}}-\chi)+\kappa\chi]\circ{\bf{A}}_{j}^{H}\theta^{{\dagger}})-{\bf{A}}_{j}^{H}[\chi_{k}\circ({\bf{u}}_{s,k})_{j}\
\circ\theta^{{\dagger}}]$ $\displaystyle=0$ (20)

where we have introduced the abbreviation
$\Pi^{{\dagger}}_{i}={\bf{u}}^{{\dagger}}_{i}+{\bf{A}}^{H}_{i}p^{{\dagger}}$.
We also obtain the terminal conditions (which result from integration by parts
in the time-domain) ,

$\displaystyle\bm{u}^{{\dagger}}(\bm{x},T)=0,\qquad\qquad\theta^{{\dagger}}(\bm{x},T)=\frac{2}{|V|}\left({\bf{A}}_{i}^{-2/3}\right)^{H}\left({\bf{A}}_{i}^{-2/3}\theta(\bm{x},T)\right).$
(21)

This set of equations has to be solved backward in time using the same
computational framework as introduced for the forward problem.

The final ingredient for the optimization problem is the first variation of
the augmented Lagrangian ${\cal{L}}$ with respect to the control variable. Two
cases shall be considered in this article: (i) the velocity of the stirrers
along their circular path, and (ii) the cross-sectional shape of the stirrers.
These two cases are treated by altering the variables $\omega_{k}$ and
$\chi_{k}$, respectively. Either case produces a gradient of the cost
functional that ultimately will be used to guide an unconstrained, gradient-
based optimization algorithm towards an enhanced solution with an improved
value of the cost functional (mixedness measure).

###  III.2 Gradient for velocity optimisation

The first case, optimizing the velocity of the stirrers along their circular
path, results from the first variation of the augmented Lagrangian with
respect to $\omega_{k}$ and yields an expression for the gradient that
involves adjoint and direct variables. After some algebra we obtain

$\displaystyle\omega_{k}^{\dagger}=r_{k}\left[l_{j}(\varphi(t))+\frac{\omega_{k}}{\dot{\omega}_{k}}\frac{\partial
l_{j}}{\partial\varphi}\right]\chi_{k}^{H}\left(\frac{\Pi^{{\dagger}}_{j}}{C_{\eta}}-(\theta^{{\dagger}}\circ[\bm{A}_{j}\theta])\right)$
(22)

for the adjoint rotational speed, which is directly linked to the sought-after
cost-functional gradient. Equation 22 then provides the gradient information
to our iterative optimization scheme, and, together with a line-search
routine, we march the current stirring protocol towards a more effective one.
The optimization terminates when no further progress can be made, or the
magnitude of the cost functional gradient falls below a user-prescribed
threshold.

###  III.3 Gradient for shape optimization

The optimization of the cross-sectional stirrer shape is somewhat more complex
and requires more computational steps. A detailed discussion of shape
parameterization, regularization techniques and their embedding into an
optimization scheme, was presented in [33], and we refer the reader to this
reference for further details. Here, we will briefly touch upon the main
concepts.

We pursue a shape parameterization that is computationally efficient, yet
sufficiently flexible to cover a rich design space. In view of the Fourier-
based computational framework for the governing equations, we also express the
cross-sectional shapes in terms of truncated Fourier series. We define each
shape by a two-dimensional parametric curve $f$ given by

$\displaystyle f_{x}(\alpha)$
$\displaystyle=\frac{a_{0}}{2}+\sum_{k=1}^{n}a_{k}\cos\left(k\alpha\right)-b_{k}\sin\left(k\alpha\right),$
(23) $\displaystyle f_{y}(\alpha)$
$\displaystyle=\frac{c_{0}}{2}+\sum_{k=1}^{n}c_{k}\cos\left(k\alpha\right)-d_{k}\sin\left(k\alpha\right),$
(24)

where $\alpha$ is the independent variable with $\alpha\in[0,\ 2\pi].$ The
cross-section is thus determined (and low-dimensionally parameterized) by
choosing sets of four coefficients according to, i.e.,
$f(\\{a_{i}\\},\\{b_{i}\\},\\{c_{i}\\},\\{d_{i}\\})=[f_{x},f_{y}].$ With the
curve thus defined, we need to determine the mask $\chi$, or, more
specifically, design an efficient way of determining whether a point on the
computational grid is contained within the closed parametric curve $f.$ To
this end, we introduce the winding number $w_{k}(\bm{x})$ of the $k$-th
stirrer which, for a discrete setting, is defined as

$\displaystyle w_{k}(\bm{x})$ $\displaystyle=\sum_{i=1}^{n}\phi_{i,i-1}.$ (25)

In the above expression $\bm{x}$ defines a point of interest, and $\phi_{i,j}$
denotes the angle of the arc that runs along the parametric curve of stirrer
$k,$ from $\bm{x}_{i}$ to $\bm{x}_{j}$ on the curve. As we sweep the full
parametric path, and obtain the full sum of angles, the total angle traversed
will be non-zero multiples of $2\pi$ if the point $\bm{x}$ is inside the
curve, or zero otherwise. When coupled to our stirrers’ mask $\chi_{k},$ the
grid projection can be stated as

$\chi_{k}(\bm{x},t)=\begin{cases}1,\hskip 12.0pt&w_{k}(\bm{x})=2\pi,\\\
0,\hskip 12.0pt&w_{k}(\bm{x})=0,\\\ \end{cases}$ (26)

where we have neglected the case of winding numbers higher than one.

As this formulation is discontinuous, simulation relying solely on this
definition of the mask $\chi$ would result in spurious oscillations and
numerical instabilities. Hence, one final modification consists of a smoothing
layer and, by extension, an applicable distance function. For simplified
geometries, such as in [35] or [30]), for example, it is sufficient to simply
measure the distance from the centroid of the cross-sectional shape. For more
complicated shapes, however, we must explicitly compute the minimum distance
from our chosen point to our curve, following

$\displaystyle d=\min_{\alpha}\sqrt{(x-f_{x})^{2}+(y-f_{y})^{2}}.$ (27)

The expression for the full mask $\chi_{k}$ is thus given as

$\chi_{k}(x,t)=\begin{cases}1,\hskip 12.0pt&w_{k}=2\pi\ \text{and}\ d>h,\\\
\frac{1}{2}\left(1+\cos\left(\frac{\pi(h-d)}{h}\right)\right),&w_{k}=2\pi\
\text{and}\ d<h,\\\ 0,\hskip 12.0pt&w_{k}=0.\\\ \end{cases}$ (28)

with $h$ as the grid size. With this geometric representation in place, we now
have all necessary ingredients to formulate the shape gradient for our
adjoint-based optimization framework.

For clarity, we present the formulation for only one of the four Fourier
coefficient $a_{i}$, but stress that the effectiveness of the shape
optimization rests in the ability to act on all control parameters
simultaneously. Since the augmented Lagrangian is not explicitly dependent on
the Fourier coefficients, we invoke a chained sequence of gradients to link
the mask $\chi_{k}$ to the Fourier coefficient $a_{i}.$ We have

$\displaystyle\frac{\partial\mathcal{L}}{\partial a_{i}}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\chi_{k}}\frac{\partial\chi_{k}}{\partial
d}\frac{\partial d}{\partial a_{i}}.$ (29)

The first term can be obtained from equation (17). The second term in this
equation derives from (28) and can be seen in analogous form in [30]. For the
last term, we note that equation (27) does include a minimisation, however the
derivative term in equation (29) is independent of this minimum, which
consequently does not affect the subsequent derivation. Given the form of the
parameterization, the derivative is in essence isolating the parameter
$\tilde{\alpha}$ within the minimum distance function, i.e.,

$\displaystyle\frac{\partial d}{\partial
a_{i}}=\frac{(f_{x}-x)\cos(i\alpha)}{d}\biggr{|}_{\tilde{\alpha}}.$ (30)

Combining these terms, we arrive at the full gradient. The composite nature of
the above formulation means that many Fourier coefficients (and thus complex
shapes) can be computed with relative ease and little computational effort:
the first two terms can be computed off-line for all Fourier coefficients,
while only the last term changes during the optimization and needs to be
continuously updated. Consequently, the optimization procedure is numerically
efficient and readily scalable, as the majority of the terms can be calculated
a-priori.

Following the computation of the shape gradient, we may have to impose further
constraints to avoid unrealistic or undesired solutions. These additional
constraints on the cross-sectional shapes will be dealt with via projections:
once an optimized shape has been determined, the resulting parametric curve is
projected onto the set of permissible solutions. This enforces feasibility of
the enhanced stirrers – of course, at the cost of full optimality. Perhaps the
most obvious constraint is that of area conservation. Both extremes, i.e.,
vanishingly thin and fragile stirrers or stirrers that occupy the bulk of the
computational domain, must be avoided. This is easily accomplished by dividing
the Fourier coefficients of the optimized shape, with an area of $A,$ by a
factor of $\sqrt{\frac{A}{A_{0}}}$, where $A_{0}$ refers to the area of the
original starting shape. In this manner, we only allow a re-arrangement of
material, rather than a reduction/augmentation.

A further imposed constraints strongly penalizes self-intersecting parametric
curves which would lead to topological pinch-off. In fact, a related
constraint avoids thin and frail features in the cross-sectional shapes that
may fracture in real-world settings. This latter constraint is enforced by an
untwisting and expansion routine based on intramolecular (e.g., Buckingham)
potentials distributed along the optimized parametric curve.

###  III.4 Summary of optimization procedure

Before turning to the presentation of the result we summarize the overall
optimization procedure. We start by initializing the mixing strategy with an
initial guess for the velocity (commonly a constant velocity along the
circular path) and shape of the embedded stirrers (commonly a circular cross-
section). Using this guess, we then solve the governing equations (1) - (3)
forward in time over the chosen time horizon $[0,T]$. Due to the nonlinear
nature of our governing equations, which results in a dependence of the
adjoint variables on the corresponding direct state, it is necessary to retain
or recover state variables at every time step (see comments below). Once we
have arrived at $t=T$, we initialize the adjoint variables with the
appropriate terminal condition and commence their backward propagation from
$t=T$ to $t=0$, injecting the direct state variables at the appropriate times.
The trajectory of the adjoint variables is then used to evaluate the
optimality condition and retrieve the gradient information of the cost
functional with respect to either the path velocity and cross-sectional shape
for each stirrer. These gradients can then be utilized by a standard
optimization routine, such as, e.g., steepest descent or conjugate gradient,
to arrive at updated, and improved, stirrer shapes and/or path velocity
signals. Finally, these updated variables are then projected onto the feasible
design space, such that the accepted variables do not violate any of our
imposed constraints. This forward-backward-update-projection loop is repeated
until either no further improvements can be achieved or further steps would
unavoidably violate constraints.

Finally, it is worth commenting on the computational cost associated with the
iterative direct-adjoint part of the optimization scheme. As alluded to above,
the nonlinear nature of the governing equations produces a link between the
adjoint governing equations and the direct state variables, which appear in
the form of varying coefficients in the otherwise linear adjoint equations.
Saving a large number of multi-dimensional arrays may exceed the available
memory resources. To mitigate this bottleneck we employ checkpointing, which
allows us to trade excessive memory requirements for slightly longer
simulation times. This strategy involves saving a set of state variables at
particular points, referred to as checkpoints, which are subsequently used as
initial conditions to recover the necessary flow fields in-between the
checkpoints. Only storage of flow fields for a typical interval between
checkpoints and the checkpoints themselves is required, which is typically far
smaller than storing the flow fields over the entire simulation horizon. An
efficient implementation of an optimized placement of checkpoints and the flow
field management during the backward sweep is given by the revolve library
[36].

##  IV Results

The iterative optimization algorithm is now applied to enhance mixing of a
binary fluid by stirrers on circular paths. We consider a two-dimensional
configuration where the binary fluid fully occupies a circular vessel, is
initially separated and confined to the lower and upper half, respectively.
Two stirrers on circular paths of different radii then move through the fluid
and cause the binary fluid to mix. For the simulations, we choose a Reynolds
number of $Re=100$ and a Péclet number of $Pe=10000.$ This choice of
parameters is a departure from earlier studies [30, 33] where advective
mechanisms dominated viscous diffusion in the optimal strategies. The current
choice requires a far greater effort to achieve mixing enhancement and thus
poses greater challenges to the optimization scheme. The optimization horizon
has been fixed to $T=12$ non-dimensional time units. Over this time span,
information on the evolving flow field and passive scalar is available for an
optimization of the mix-norm at the end of the time horizon. Beyond $t=12,$
the stirrers will cease to move and become passive elements within the mixing
vessel; rest momentum will continue to prevail, but increasingly be dissipated
by viscous diffusion. Despite our focus on this single parameter
configuration, we stress that the direct-adjoint optimization proves efficient
over a wide range of alternative $Re$-$Pe$-values. There are, however,
limitations on the choice of the optimization horizon $T,$ which we will
discuss later.

|
---|---
|

(a) unoptimized

(b) shape optimization

(c) velocity optimization

(d) combined optimization

Figure 1: Comparison of final snapshots at the end of the optimization horizon ($t=T=12$), visualized by the passive scalar $\theta:$ (a) unoptimized base case, (b) area-preserving shape optimization of stirrers’ cross-sections, (c) path-velocity optimization of the shape-modified stirrers along circular paths, (d) composite optimization of cross-sectional shapes and path-velocities. |  |  | iterations |  | $\|\theta\|_{{\rm{mix}},t=12}$
---|---|---|---|---|---
(i) | unoptimized base case |  | n/a |  | 0.168
(ii) | shape optimization |  | 11 |  | 0.161
(iii) | velocity optimization |  | 8 |  | 0.105
(iv) | combined optimization |  | 8 |  | 0.072
Table 1:  Summary of results for all numerical experiments. Number of
iterations taken by the direct-adjoint optimization algorithm are displayed,
where applicable, as well as the values of the final mix-norm of the passive
scalar $\theta$ at $t=12.$

To showcase the efficacy and applicability of the direct-adjoint framework, we
present three distinct optimization cases: _(i)_ optimization of the cross-
sectional shape of the moving stirrers, _(ii)_ optimization of the velocity
along a circular path of the shape-modified stirrers (from step _(i)_) and
_(iii)_ compound and simultaneous optimization of both the shape and path
velocity for both stirrers. The constraints include energy, velocity and
energy bounds for any path velocity optimization [see, e.g., 31], area-
conservation and untwisting for the shape optimization [see, e.g. 33] and both
sets of constraints for the third case. In all cases presented here, we
achieve enhanced mixing and a significant decrease in the mix-norm at the end
of the optimization horizon, attesting to the effectiveness of an objective
optimization approach to reach a more homogeneous mixture through unorthodox
stirrer designs and improved stirring protocols. The final snapshots (at
$t=T=12$) for each of the four scenarios – unoptimized base case, shape
optimization, path-velocity optimization of the optimized shape, and compound
optimization – are presented in figure 1. These snapshots show an increasing
amount of thin filaments appearing as a result of stirring action. This
qualitative analysis is further corroborated by the quantitative values listed
in table 1, where we observe a monotonic decrease of the mix-norm of the
passive scalar $\theta$ from the initial (unoptimized) base case to the
combined shape/velocity optimization. In what follows, we will present the
results of each optimization case in more detail and analyze the stirrer-
induced vortex motion that yields enhanced mixing strategies.

###  IV.1 The base case: cylindrical stirrers with uniform path velocities

Our reference case, against which we measure any mixing improvement, consists
of two stirrers of circular and identical cross-section that move along
concentric circular paths (of given outer and inner radius) with a constant
speed. Starting in the twelve-o-clock (for the outer radius) and six-o-clock
(for the inner radius) position, the two stirrers complete a full revolution
over the chosen time horizon. As the two stirrers plunge through the
interface, separating the two components of the binary fluid, they deform the
interface. Moreover, wake vortices, shedding continuously along their paths,
further cause and promote the breakdown into small-scale features which
ultimately dissipate (see Control.mp4 from the supplemental material).
Snapshots of the base-case scenario can be seen in figure 3 as the reduced-
size subplots attached to the corresponding optimized configurations. They
confirm a notable mixing of the binary fluid – but simultaneously leave room
for improvement.

###  IV.2 Optimization of the stirrers’ shape

(a) | (b)
---|---
Figure 2: Optimization of the cross-sectional shape of both stirrers, leaving
the velocity along the circular paths unchanged. The red and green curves
signify the initial and final configuration, respectively, with gray shapes
denoting the intermediate steps of the optimization procedure. Blue arrows
have been added to indicate the travel direction of the two stirrers. Black
arrows pinpoint features of interest, as they develop during the optimization
process. (a) The outer stirrer shows edges on the leading edge, up to the
location perpendicular to the travel direction, and a pronounced splitter-
plate-like appendage at the trailing egde. (b) The inner stirrer shows a more
complex shape, with a marked appendage perpendicular to its travel direction,
causing a significant increase in the cross-sectionally exposed area. (a) $\
t=2$(b) $\ t=4$(c) $\ t=6$(d) $\ t=8$(e) $\ t=10$(f) $\ t=12$ Figure 3:
Selected snapshots of the passive scalar $\theta$ after $11$ iterations of the
direct-adjoint looping, optimizing the cross-sectional stirrer shapes to
enhance the mix-norm at $T=12.$ The corresponding snapshots for the base case
with cylindrical stirrers are juxtaposed to the lower right of each snapshot
(in reduced size) for comparison.

We begin by presenting the first optimization case: optimizing the cross-
sectional shape of the two moving stirrers. Given that the velocity will
remain unchanged in this case, no considerations about energy constraints need
to be made. On the other hand, we remark that fixing both path and velocity
renders the rise of any collaborative strategies nigh impossible. As we will
confirm later, these collaborative strategies are highly effective in mixing
enhancement, and the limitation to only shape optimization will be reflected
in relatively modest improvements in the mix-norm over the optimization
horizon (see table 1). When optimizing the stirrers’ shape, the focus is
rather on mixing enhancement in the neighborhood of each stirrer. This
enhancement takes the form of sharp-edged appendages that (i) transport more
material between unmixed sections of the fluid, (ii) decrease the aerodynamic
cross-sectional profile of the stirrers, (iii) promote the generation of
vortices shedding from sharp corners and edges, and (iv) modify the trailing
edge of the stirrers to increase the complexity of the trailing wake. The
changes in shape during the optimization process are displayed in figure 2,
where we observe the realization of key features mentioned above. While the
outer stirrer concentrates on the generation of modified wake vortices via an
extended trailing edge and sharp edges perpendicular to the travel direction,
the inner (slower) stirrer adopts a bluff-body strategy by markedly extending
in the direction perpendicular to its travel direction. The resulting shape
will foster an increased shedding activity and an injection of substantial
vorticity into the flow. The effects induced by these modified stirrers play
an important role during various periods of the stirring process, and can be
most clearly observed in the animations (see the supplemental material); we
urge the reader to peruse the videos. These videos will also contain the
dynamics beyond the optimization horizon ($T=12$) and show the subsequent
passive mixing by rest inertia and diffusion. While we cannot fully do justice
to this intrinsically dynamic process, we nonetheless will illustrate the
stirring of a binary fluid using optimized cross-sectional stirrer shapes with
the help of the static images in figures 3 and 4.

$(a)$ | $(b)$ | $(c)$
---|---|---
Figure 4: Examples of increased passive scalar transportation with shape-
optimized stirrer. (a) Modification of the inner (low-velocity) stirrer,
showing an elongation perpendicular to its direction of motion. (b)
Interaction of shed vortices generated by the perpendicular edge of the outer
stirrer and the trailing edge protrusion, respectively, with arrows indicating
the path of the vortices. (c) Generation of thin filaments and vortex
interaction in the wake of the enhanced shape of the outer stirrer. See
ShapeOpt.mp4 from the supplemental material for an animation of the mixing
process.

Vortex generation and vortex interaction are the key fluid-dynamical tools for
optimal mixing. The rather low velocity of the inner stirrer (closer to the
center of the vessel) precludes the formation of significant wake vortices,
let alone instigate their interaction. In light of this impediment, the inner
stirrer enhances mixing by morphing into a scoop-like shape that aids in
transporting homogeneous material into high-shear zones simply by a plunging
action. While suboptimal to vortex-interactions (as will be shown later), this
strategy results in modest, but noticeable improvements of the mix-norm over
our specified optimization window.

As observed above, over the course of $11$ direct-adjoint iterations, the
shape of the inner stirrer develops an extended appendage that introduces a
paddle-like feature as well as an associated cavity (see figure 4(a)). The
extended reach of this appendage also yields, via the locked rotation of the
stirrer about its centroid along the circular path, a marginally larger
velocity for shedding vortices, when compared to the initial circular cross-
section. The effect of this topology change is marked and can be observed when
comparing the snapshots in figures 3(a) and (b) to their base-case
equivalents. Already at this early stage in the optimization window, the
interface between the binary fluids has deformed more significantly, and long,
thin fluid structures have formed that are absent in the corresponding base-
case solution. The formation of long and thin filaments is crucial in
preparation of a mostly diffusive process, once the stirring stops. At $t=8$
(see figure 3(d), or figure 4(a) for a close-up), the advantage of an enlarged
normal cross-section becomes evident: the wake structure of the inner stirrer
contains significantly finer structures when compared to the wake of the
initial shape.

The outer stirrer (closer to the vessel wall) does possess the necessary
velocity to take advantage of a wake-vortex-interaction strategy in order to
enhance mixing. As a consequence, and in contrast to the bluff inner stirrer,
it develops a more streamlined cross-sectional shape over the course of the
optimization process. Nevertheless, two edges are visible at about $\pm
90^{\circ}$ to the traveling direction which are responsible for a continuous
and vigorous vortex shedding (see figures 3(c-f)) and the injection of small-
scale structures into the stirred fluid.

We furthermore observe that the outer stirrer favors a slender protuberance
near its trailing edge, roughly aligned with the direction of motion. This
solid feature acts as a splitter plate for the wake, and promotes vortices of
increased complexity in the flow past the stirrer. The slightly asymmetric
positioning of this elongated appendage introduces oblique vortex pairs of
differing strengths (see. e.g., figures 3(c,d) and figure 4(b)) which
subsequently interact to form thin vortex filaments in the wake of the outer
stirrer (see figures 3(e,f), and figures 4(b,c) for a close-up). The
combination of multi-edged shedding of vortices and their interaction via the
split wake has been identified by the direct-adjoint scheme as an effective
means to decrease the mix-norm over the optimization window, to introduce
small-scale structures and to promote homogeneity of the passive scalar. In
this effort, the wake profile induced by the shape-modified stirrer is
considerably more complex than the wake profile for the cylindrical stirrers
(as can be seen from figures 1 and 4c).

In summary, the two stirrers have been optimized in their cross-sectional
shape according to the velocity regime in which they operate. The inner
stirrer, in a lower-speed environment, has adopted a drag-based strategy
whereby an appendage perpendicular to its direction of travel reaches into a
region of larger velocity (increased distance from the center) and induces
vortex shedding from its edge. An associated scooping motion, coupled to a
secondary motion brought about by a cavity on the stirrer, affects the mixing
towards lower mix-norm values. In contrast, the outer stirrer – with its
position in a higher-velocity range – selects a strategy of injecting vortices
into the binary fluid via an appropriately manipulated wake profile. Edges on
either side of the stirrer at strategic locations generate vortices which then
interact with an elongated trailing edge (and each other) to modify the wake
into a more convoluted vortex street with intensified smaller scales. Despite
these encouraging efforts to promote mixing, we see little evidence of
collaborative strategies between structures induced by the two stirrers. For
this characteristic to arise, we will need to allow for individual velocity
protocols for each stirrer; this will be the case for the next optimization.

We conclude this section with the remark that optimizing the shape for a fixed
velocity along the circular paths does not alter the energy budget (even
though an additional effort may be required due to an increased drag on the
stirrers). As such, optimization of the stirrers’ cross-sectional form may be
an attractive alternative to augment mixing efficiency when an increase in
energy is not feasible or advisable.

###  IV.3 Velocity optimization for enhanced stirrer shapes

After the cross-sectional shapes of the stirrers have been modified to enhance
mixing, the natural next step in our quest for an improved mixing strategy is
to optimize the path velocities of the modified stirrers along their circular
trajectories. As the optimization results of this section are even more
dynamic in nature than in the previous section, we once more prompt the reader
to consult the relevant animations in the supplemental material for a full
appreciation of the enhanced mixing process.

(a)$t\in[0,\ 3]$ |  | (b)$t\in[3,\ 6]$
---|---|---
(c)$t\in[6,\ 9]$ |  | (d)$t\in[9,\ 12]$
Figure 5:  Velocity trace of the two stirrers, visualized by position samples
at equispaced intervals in time. Each time the stirrers change direction, the
radius is reduced by a small percentage to avoid visual overlay. The true
stirrers, of course, remain on their circular paths. The entire optimization
horizon is broken down into to equal temporal segments.

Optimizing the velocity along the circular path introduces the possibility of
creating numerous start and stop vortices. While in the previous case, only
one start vortex (at $t=0$) and one stop vortex (at $t=T$) resulted, we are
now in a position to inject multiple vortical structures into the interior by
manipulating the acceleration of the stirrers within the optimization window.
In this manner, we produce ’hydrodynamic’ stirring elements, beyond the two
obvious solid stirrers. As a consequence of this degree of freedom, we have to
impose additional bounds on the energy, velocity and acceleration of the
stirrers to avoid divergences in the optimization scheme and, instead, arrive
at a converged and realizable mixing strategy (see [33] for details). But even
with these bounds in place, we achieve a dramatic decrease in the mix-norm,
indicating a superior mixing result.

Once stirrers are allowed to alter their position (and velocity) in time along
the path, cooperative strategies form where stirrers generate start-stop
vortices which get either sheared by the wake or split by the appendages of
the modified stirrers. Even after the stirrers cease to move, they often
position themselves in advantageous locations and act as obstacles for the
moving vortices. The consequence of either process is a swift generation of
small-scale fluid structures, combined with a rapid breakdown of the
heterogeneous mixture into a homogeneous blend.

Figure 5 displays the velocities along the circular paths for both stirrers.
The visualization is based on a equispaced sampling in time, such that widely
spaced symbols indicate fast motion, while densely spaced symbols point to
slower motion. To conceptualize reverse motion, we artificially decrease the
radius of the path by a small amount to avoid overlays. Figure 5 breaks down
the entire optimization horizon into four segments. The first segment
$t\in[0,\ 3]$ shows a co-rotating outer and inner cylinder, with a strong
initial acceleration of the outer cylinder followed by a noticeable slow down
and an initial shaking motion of the inner cylinder. This segment is followed
by an oscillatory motion (indicated by the varying density of the temporally
equispaced points) of the outer cylinder, and a similar motion of the inner
cylinder. The two later segments (see figure 5(c,d)) are characterized by a
convulsive and reverse motion which (as we will see below) introduces start-
stop vortices and interacts with the already generated vortices in the vessel.
This motion is maintained until the end of the optimization window, when the
two stirrers finally come to rest.

|  |
---|---|---
|  |

(a) $t=2$

(b) $t=4$

(c) $t=6$

(d) $t=8$

(e) $t=10$

(f) $t=12$

Figure 6: Selected snapshots of the final velocity optimization for the
modified stirrer shapes. (a) Initial plunging through the interface by the
outer stirrer, while the inner stirrer generates start-stop vortices. (b)
Oscillatory motion of the outer stirrer and acceleration of the inner stirrer.
(c,d) Reversal of both stirrers, allowing interaction of the stirrers with the
generated vortices. (e,f) Small-amplitude rapid oscillation of the outer
stirrer, while the inner stirrer gradually comes to rest. See VelOpt.mp4 from
the supplemental material for an animation of the mixing process.

With the modified shapes, the velocity optimization of the two stirrers yields
a substantial decrease in the mix-norm and thus a better mixing result. Figure
6 depicts snapshots of the mixing process. The outer stirrer accelerates
towards the interface and plunges through it. In the process, it creates wake
vortices which further distort the interface. Towards the first quarter of the
optimization window, the outer stirrer decelerates and engages in an
oscillatory forward motion near the bottom of the mixing vessel. This unsteady
motion induces enhanced vortex shedding which will ultimately produce thin
filaments in the domain. During the same time window, the inner stirrer has
generated start-stop vortices by initiating its trajectory with a rapid
oscillatory and low-amplitude motion. The thin cross-stream appendage is used
to full advantage in injecting vortical elements which will subsequently be
further distorted by the outer stirrer (see figure 6(b)).

$(a)$ | $(b)$ | $(c)$
---|---|---
Figure 7: Closeups of key features of the velocity-optimized stirring
protocal. (a) Shedding of stop-vortex by the outer stirrer, aimed for impact
on the inner stirrer. (b) Breakdown of large-scale vortical elements via the
appendage of the outer stirrer. (c) Generation of shed vortices by oscillatory
motion of the outer cylinder, near the end of the optimization horizon.

The goal of the optimization is the generation of a maximum of vortical
element – and their breakdown into smaller vortices – within the optimization
window. Beyond this time window, the introduced vortices will continue to
interact with each other, the vessel wall and the stationary stirrers and
cause further breakdown in scale due to rest inertia and ultimately diffusion.
In this effort, the direct-adjoint optimization scheme takes advantage of the
stirrers’ edges, appendages and cavities to devise a velocity profile along
the circular paths that maximizes vortex-vortex and vortex-body interactions
to promote the breakdown of scales and the enhancement of mixedness. The
optimal strategy is highly collaborative and, as a result, achieves
significantly better mixing results than a simple shape optimization (without
a concurrent velocity optimization). This tendency is the consequence of a far
larger design space that is exploited by the direct-adjoint framework to
produce more effective and efficient mixing protocols.

###  IV.4 Combined shape-velocity optimization protocol

The previous configuration consisted of a sequential application of two
optimizations, for the cross-sectional shape of the stirrers and the velocity
profile along circular paths for the modified shapes. In this effort, the
velocity profile has been tailored to the modified shape, but no feedback from
current results to the stirrers’ shapes, via the shape gradient, has been
allowed. In this final section, we ’unfreeze’ the stirrer shapes and hence
pursue a simultaneous optimization of the shapes and path velocities of the
two stirrers. In this sense, we seek a multi-component optimum from our
direct-adjoint framework.

Both gradients – the first based on the path velocity for each stirrer, and
the second based on the cross-sectional shape of each stirrer – aim at
improving the mix-norm over the chosen optimization horizon. As these
gradients enter the optimality condition, a choice has to be made as to how
much of the improvement will come from a modified shape and how much will stem
from an altered velocity protocol. Previous observations have confirmed a
noticeably stronger overall response to velocity changes; topology changes, on
the other hand, imposed a more subtle, but still perceptible effect on the
mixedness of the binary fluid. A straightforward application of the direct-
adjoint scheme then would favor an optimization that is chiefly focused on
velocity changes, while discounting modifications of the stirrers’ shapes. To
counterbalance this intrinsic tendency and investigate the complex interaction
of shape and velocity variations, we choose to underrelax the velocity-based
gradient and give equal prominence to the shape-based gradient.

(a) | (b)
---|---
Figure 8: Optimal cross-sectional shapes for both stirrers, using a combined
shape-velocity optimization. The red and green curves signify the initial and
final configuration, respectively, with gray shapes denoting the intermediate
steps of the optimization procedure. Black arrows pinpoint features of
interest, as they develop during the optimization process. (a) The outer
stirrer develops an irregular pentagon shape, with one edge more pronounced.
(b) The inner stirrer shows similar shape with one conspicuous extended blade.

Figure 8 shows the evolution of the cross-sectional shapes of both stirrers
over the course of the optimization. Starting from a circular cross-section,
the optimization of the outer stirrer (figure 8(a)) quickly develops a more
complex shape with five corners, one of them more pronounced than the
remaining four. The gray lines in the figure indicate a robust convergence
towards a final shape after eight iterations. The inner stirrer (figure 8(b))
undergoes a similar transformation from an initially circular cross-section to
a more complicated shape, with more edges and a blade-like and curved
protrusion. The gray lines, visualizing the intermediate optimization steps,
show that in this case a regularization step had to be invoked to maintain a
minimal thickness of the cross-sectional shape. Nonetheless, the convergence
appears to be robust in this case as well.

(a)$t\in[0,\ 3]$ |  | (b)$t\in[3,\ 6]$
---|---|---
(c)$t\in[6,\ 9]$ |  | (d)$t\in[9,\ 12]$
Figure 9:  Velocity trace of the two stirrers for the combined shape-velocity
optimization, visualized by position samples at equispaced intervals in time.
Each time the stirrers change direction, the radius is reduced by a small
percentage to avoid visual overlay. The true stirrers, of course, remain on
their circular paths. The entire optimization horizon is broken down into to
equal temporal segments.

The corresponding velocities for the two stirrers – which have been found
simultaneously with the evolving shapes – are displayed in figure 9 for the
final (eighth) iteration. Again, the optimization horizon has been broken down
into four temporal segments, for easier visualization. Figure 9(a) shows that,
in contrast to the previous case, the two stirrers run counter to each other,
with the outer stirrer moving counter-clockwise and the inner cylinder
clockwise. This choice strongly suggests the pursuit of a collaborative mixing
strategies where both stirrers cooperate to introduce dynamic structures into
the vessel during the early stages of the optimization horizon. Over the first
quarter of the full window, the outer stirrer progresses towards the
interface, but significantly slows down as it plunges through the interface
(as evidenced by the higher density of symbols on the left edge of figure
9(a)). This strategy causes shed vortices to overtake the stirrer and produces
an enhanced distortion of the interface as a consequence. The inner stirrer
follows a similar protocol, stopping and starting on its way towards the
interface. The counter-moving travel of the stirrers, together with the
extended shapes, produces a close encounter of the stirrers and a strong and
localized shear layer.

The remaining temporal segments provide evidence of further vortex shedding,
arising from the start-and-stop movement of the outer as well as inner
stirrer. The liberal use of multiple and rapid changes in direction attest to
the optimality of this strategy. The direct-adjoint scheme simply exploits
this strategy and launches a multitude of shed vortices that accomplish the
breakdown of heterogeneities during and – more importantly – beyond the
optimization horizon (when the physical stirrers no longer move).

|  |
---|---|---
|  |

(a) $t=2$

(b) $t=4$

(c) $t=6$

(d) $t=8$

(e) $t=10$

(f) $t=12$

Figure 10: Selected snapshots of the final combined shape-velocity
optimization. (a) Initial approach of the interface by both stirrers, with
deceleration-acceleration near the interface. (b) Acceleration of outer
stirrer, deceleration of inner cylinder. (c) Oscillatory motion of the outer
stirrer and near stagnation of inner stirrer. (d-f) Multiple reversals of
outer stirrer and slow motion of inner stirrer, as the stirrers come to their
respective resting position. See ComboOpt.mp4 from the supplemental material
for an animation of the mixing process.

The success of this strategy can be assessed by inspecting the dynamic process
over the course of the optimization window and beyond. Again, the animated
version of this process (in the supplemental material) furnishes a far better
perspective on this enhanced mixing. The still pictures in figure 10 show that
better mix-norm results have been obtained; even a qualitative assessment of
the generated structures suggests a superior final result, when compared to
previous cases. During the early movements of the stirrers, the close
encounter of the two stirrers near $t=2$ is noteworthy. It creates a strong
shear layer that is responsible for much of the subsequent interface
distortion. The intensity of this shear layer is further aggravated by the
elongated features of the stirrers which ensures an even closer brush of the
two elements. Already at $t=4$ (see figure 10(b)) the initially flat interface
is markedly contorted. The mixing process is subsequently supported by a
complex stop-start sequence of both stirrers which further brings about the
injection of spinning vortices into the stirred-up mixture. The interplay of
all vortical components, large and small, and the interaction with the solid
stirrers and the vessel wall causes a prevalence of thin filaments (see figure
10(e,f)) which ultimately succumb to diffusion.

$(a)$ | $(b)$ | $(c)$
---|---|---
Figure 11:  Closeups of key features of the shape-velocity-optimized stirring
protocol. (a) Narrow passage of two stirrers, creating a strong shear layer
and a subsequent roll-up of the interface. (b) Guiding and deflection of
vortices by the inner stirrer. (c) Breakdown of vortices as they interact with
the sharp protrusions of the outer stirrer.

Figure 11 highlights three key features from the optimal dynamic mixing
process. The first detail (figure 11(a)) showcases the narrow brush of the two
stirrers early during the optimization horizon, where the two shapes nearly
touch each other and pinch the interface between them. The resulting shear
layer continues the mixing process in this section of the vessel long after
the stirrers have moved towards other parts of the domain. The second feature
(figure 11(b)) demonstrates the impingement of vortical structures on the
bespoke shapes of the stirrers, in this case the inner stirrer. The blade-like
protrusion breaks down the scale of the incoming vortex, deflects its travel
direction and induces in its cavity a counter-rotating smaller vortex which
eventually detached from the cavity. A similar scenario can be observed in
figure 11(c) where a pair of sizeable stop-start vortices directly impinge on
the bluff edge of the lower stirrer. Near the time of impact, the stirrer
shifts abruptly to the left, intensifying the collision and reinforcing the
subsequent breakup into smaller scales. This general scenario of vortex-
generation and vortex-splitting can be observed throughout the latter part of
the optimization window, and constitutes a particularly effective mechanism
chosen by the optimization scheme to enhance mixing.

In summary, the composite and simultaneous optimization of the velocity
protocols and stirrer shapes yielded a promising strategy for mixing
enhancement. It is based on a cooperation of the two stirrers and an intricate
interaction of multiple vortices at the local as well as global level. It
should be kept in mind, however, that a concentration on the velocity-
component of the optimization may produce the bulk of the gains in mixedness.
Nonetheless, the transformed shapes of the stirrers can contribute additional,
albeit more modest benefits.

##  V Summary and conclusions

A direct-adjoint optimization framework has been applied to enhance the mixing
of a binary, incompressible, Newtonian fluid. Bound in a circular vessel, the
fluid is mixed by two stirrers moving on two concentric circles. The cross-
sectional shape of the stirrers, parameterized via Fourier series, and the
velocities along the circular paths are optimized – in series or combined – to
enhance the homogeneity of the passive scalar field, while observing a finite
time horizon, a finite energy budget, and bounds on the velocities and
accelerations [31]. Despite an exceedingly large design space, the gradient-
based procedure arrived at effective stirrer shapes and stirring protocols
that yielded substantially improved mixing measures, expressed in the mix-norm
of the passive scalar.

Four cases have been considered: (i) a base case with cylindrical stirrers on
prescribed circular paths, against which all subsequent optimization have been
compared, (ii) an optimization of the cross-sectional shape, while carrying
forward the previous path velocities, (iii) a pure velocity-optimization for
the modified stirrer shapes, and (iv) a combined and weighted shape-and-
velocity optimization of the mixing process. Starting from the reference case,
each subsequent optimization results in a lower mix-norm over the optimization
horizon of $T=12,$ and thus a more homogeneous mixture. The accomplished
increase in mixing efficiency is summarized in table 1 which lists the mix-
norm values for all simulations, together with the number of direct-adjoint
iterations.

We observe that a pure shape optimization, disregarding the velocities of the
stirrers, will result in a noticeable but rather moderate reduction in mix-
norm. Only a decrease of $4\%$ could be achieved by manipulating the cross-
sectional shapes of the stirrers. While this improvement may translate into a
more sizeable saving of energy and efforts within an industrial setting, a far
larger gain can be accomplished by targeting the velocities of the two
stirrers. The targeting of this degree of freedom introduces the option of a
collaborative strategy between the two stirrers – an option that was largely
unavailable during the pure shape-optimization. Vortex generation and
interaction, plunging, unsteady vortex shedding, vortex collisions with
themselves, the wall and the stirrers are at the optimizer’s disposal to cast
into a comprehensive mixing strategy. This optimization may lead (and has led)
to counter-intuitive and complex stirring protocols that aid in the breakdown
of the initial structures into small-scale elements that eventually dissipate.
The efficacy of this more dynamic stirring strategy is reflected in the
substantial drop in mix-norm (see table 1): by optimizing the velocities for
the modified shapes (case (iii)), we exhibit a further $34\%$ drop in mix-norm
from the previous case. By pursuing a combined optimization of the shapes and
velocities of the stirrers, we can extract an additional increase in mixing
efficiency, compared to the sequential optimization approach (case (iii)
following case (ii)). We experience another $31\%$ decrease in mix-norm from
case (iii) by using a simultaneous shape-velocity optimization. This results
in the lowest value (and thus the best mixing result) for all cases considered
in this study; it constitutes a $57\%$ improvement over the unoptimized base
case (case (i)).

The lowest mix-norm strategy (case (iv)) contains a variety of processes that
complement each other, interact and cooperate to maximally increase the
presence of thin filaments (in the passive scalar) and the dominance of
shearing structures. The desired state of the fluid system at the end of the
optimization window is characterized by a flow configuration that is maximally
prone to rest inertia and species diffusion. The evolution of the seeded
structures beyond the optimization horizon (see the supplemental material)
shows a continued breakdown of scales and the convergence towards a maximally
homogeneous mixture. The direct-adjoint optimization scheme proved capable and
robust in finding effective and efficient stirring strategies within a small
number of iterations.

Despite this success, several challenges and opportunities remain. First, and
foremost, is the inability of any gradient-based optimization approach applied
to a nonlinear partial differential equation to produce global optimality.
This shortcoming is due to the fact that we are dealing with a non-convex
optimization problem, for which any derivative-based scheme can only guarantee
a local optimum. While this issue may limit our results, is it important to
realize that we achieved a substantial improvement (a $57\%$ reduction in mix-
norm) which, under any industrial setting, would translate into vast savings
in energy and resources.

A further challenge lies in the extension of the optimization horizon. We have
chosen the parameters of our study (notably the Reynolds number, Péclet number
and time horizon) to allow a rich dynamic range of motion, but, at the same
time, avoid the influence of turbulence. As we rely on access to gradient
information, the presence of turbulent fluid motion will obscure the
sensitivity information needed for the shape and/or velocity optimization.
While techniques to overcome this predicament (such as least-squares
shadowing, see [37, 38]) are currently investigated, their cost-effectiveness
applied to systems of the present complexity is still an open problem.

Besides these challenges, the current study demonstrated a great amount of
flexibility in its approach, which prompts a rich set of possibilities for
further optimization studies.

These studies would include, for example, optimizing the wall shape and motion
– relinquishing the circular shape and its static nature. A moving/oscillatory
and corrugated wall could join the stirrers in the generation and manipulation
of vortical structures. A similar optimization could target the path of the
stirrers, from concentric circles to more complicated closed curves.
Introducing static baffles which break up the vortical motion may constitute
yet another interesting research direction.

An extension of the fluid model to non-Newtonian constitutive behavior is an
opportunity to describe a wide range of agents used in industrial mixing
[see,e.g., 39, 40], as is the treatment of three-dimensional geometries,
stirrers and fluid motion.

The presented direct-adjoint optimization approach has shown its mettle on a
variety of inertia-dominated mixing problems for a binary fluid. It is capable
of addressing many of the above-mentioned generalizations and extensions, and
future efforts in this direction are expected to further contribute to
enhanced mixing strategies via optimized stirring protocols.

##  VI Backmatter

#### Acknowledgements

We gratefully acknowledge discussions and exchanges with Profs. Colm
Caulfield, Jean-Luc Thiffeault, Kai Schneider and Dr. Florence Marcotte.

#### Funding Statement

M.E. gratefully acknowledges funding through the Joachim Herz Stiftung.

#### Declaration of Interests

The authors declare no conflict of interest.

#### Author Contributions

Both authors created the research plan, designed the numerical experiments and
wrote the manuscript. M.E. performed the code modifications and performed the
numerical simulations.

#### Data Availability Statement

N.A.

#### Ethical Standards

The research meets all ethical guidelines, including adherence to the legal
requirements of the study country.

#### Supplementary Material

Supplementary material (animations) are available at … .

## References

  * Kukulka  _et al._ [2012] T. Kukulka, G. Proskurowski, S. Morét-Ferguson, D. Meyer, and K. Law, The effect of wind mixing on the vertical distribution of buoyant plastic debris, Geophys. Res. Letts. 39 (2012)
  * Hunt and Linden [1999] G. Hunt and P. Linden, The fluid mechanics of natural ventilation – displacement ventilation by buoyancy-driven flows assisted by wind, Buildg. and Environm. 34, 707 (1999)
  * Linden [1999] P. Linden, The fluid mechanics of natural ventilation, Annu. Rev. Fluid Mech. 31, 201 (1999)
  * Peng  _et al._ [2020] Y. Peng, M. Gerdroodbary, M. Sheikholeslami, A. Shafee, H. Babazadeh, and R. Moradi, Mixing enhancement of the multi hydrogen fuel jets by the backward step, Energy 203, 117859 (2020)
  * Thompson  _et al._ [2017] J. Thompson, S. Walton, O. Hassan, S. Rolland, and J. Sienz, The use of cfd and multi-objective optimisation techniques to customise an industrial pre-mixer, Struct. and Multidisc. Optim. 55, 2339 (2017)
  * Mihailova  _et al._ [2018] O. Mihailova, T. Mothersdale, T. Rodgers, Z. Ren, S. Watson, V. Lister, and A. Kowalski, Optimisation of mixing performance of helical ribbon mixers for high throughput applications using computational fluid dynamics, Chem. Eng. Res. and Design 132, 942 (2018)
  * Mazumdar  _et al._ [2017] D. Mazumdar, P. Dhandapani, and R. Sarvanakumar, Modeling and optimisation of gas stirred ladle systems, ISIJ Int. 57, 286 (2017)
  * Spencer and Wiley [1951] R. Spencer and R. Wiley, The mixing of very viscous liquid, J. Colloid Sci. 6, 133 (1951)
  * Ottino [1989] J. Ottino,  _The kinematics of mixing: stretching, chaos, and transport_ (Cambridge University Press, 1989)
  * Mathew  _et al._ [2005] G. Mathew, I. Mezic, and L. Petzold, A multiscale measure for mixing, Physica D 211, 23 (2005)
  * Mathew  _et al._ [2007] G. Mathew, I. Mezic, S. Grivopoulos, U. Vaidya, and L. Petzold, Optimal control of mixing in stokes fluid flows, J. Fluid Mech. 580, 261 (2007)
  * Thiffeault [2012] J.-L. Thiffeault, Using multiscale norms to quantify mixing and transport, Nonlinearity 25, R1 (2012)
  * Gao  _et al._ [2020] X. Gao, I. Bermejo-Moreno, and J. Larsson, Parametric numerical study of passive scalar mixing in shock turbulence interaction, J. Fluid Mech. 895, A21 (2020)
  * Aref [1984] H. Aref, Stirring by chaotic advection, J. Fluid Mech. 143, 1 (1984)
  * Meena and Taira [2021] M. Meena and K. Taira, Identifying vortical network connectors for turbulent flow modification, J. Fluid Mech. 915 (2021)
  * Chien  _et al._ [1986] W.-L. Chien, H. Rising, and J. Ottino, Laminar mixing and chaotic mixing in several cavity flows, J. Fluid Mech. 170, 355 (1986)
  * Kadoch  _et al._ [2020] B. Kadoch, W. Bos, and K. Schneider, Efficiency of laminar and turbulent mixing in wall-bounded flows, Phys. Rev. E 101, 043104 (2020)
  * Balogh  _et al._ [2005] A. Balogh, O. Aamo, and M. Krstic, Optimal mixing enhancement in 3-d pipe flow, IEEE Trans. Control Sys. Tech. 13, 27 (2005)
  * Gubanov and Cortelezzi [2010] O. Gubanov and L. Cortelezzi, Towards the design of an optimal mixer, J. Fluid Mech. 651, 27 (2010)
  * Lin  _et al._ [2011] Z. Lin, J.-L. Thiffeault, and C. Doering, Optimal stirring strategies for passive scalar mixing, J. Fluid Mech. 675, 465 (2011)
  * Liu [2008] W. Liu, Mixing enhancement by optimal flow advection, SIAM J. Control and Opt. 47, 624 (2008)
  * Foures  _et al._ [2012] D. Foures, C. Caulfield, and P. Schmid, Variational framework for flow optimization using seminorm constraints, Phys. Rev. E 86, 026306 (2012)
  * Foures  _et al._ [2013] D. Foures, C. Caulfield, and P. Schmid, Localized flow structures using $\infty$-norm optimization, J. Fluid Mech. 729, 672 (2013)
  * Foures  _et al._ [2014] D. Foures, C. Caulfield, and P. Schmid, Optimal mixing in plane Poiseuille flow, J. Fluid Mech. 748, 241 (2014)
  * Marcotte and Caulfield [2018] F. Marcotte and C. Caulfield, Optimal mixing in two-dimensional stratified plane Poiseuille flow at finite Péclet and Richardson numbers, J. Fluid Mech. 853, 359 (2018)
  * Vermach and Caulfield [2018] L. Vermach and C. Caulfield, Optimal mixing in three-dimensional plane Poiseuille flow at high Péclet number, J. Fluid Mech. 850, 875 (2018)
  * Wu  _et al._ [2018] Z. Wu, D. Fan, Y. Zhou, R. Li, and B. Noack, Jet mixing optimization using machine learning control, Exp. in Fluids 59, 1 (2018)
  * Hilgers and Boersma [2001] A. Hilgers and B. Boersma, Optimization of turbulent jet mixing, Fluid Dyn. Res. 29, 345 (2001)
  * Alexias and Giannakoglou [2020] P. Alexias and K. C. Giannakoglou, Shape optimization of a two-fluid mixing device using continuous adjoint, Fluids 5, 11 (2020)
  * Eggl and Schmid [2018] M. Eggl and P. Schmid, A gradient-based framework for maximizing mixing in binary fluids, J. Comp. Phys. 368, 131 (2018)
  * Eggl and Schmid [2020a] M. Eggl and P. Schmid, Mixing enhancement in binary fluids using optimised stirring strategies, J. Fluid Mech. 899 (2020a)
  * Angot  _et al._ [1999] P. Angot, C.-H. Bruneau, and P. Fabrie, A penalization method to take into account obstacles in incompressible viscous flows, Num. Math. 81, 497 (1999)
  * Eggl and Schmid [2020b] M. Eggl and P. Schmid, Shape optimization of stirring rods for mixing binary fluids, IMA J. Appl. Math. (2020b)
  * Kolomenskiy and Schneider [2009] D. Kolomenskiy and K. Schneider, A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles, J. Comp. Phys. 228, 5687 (2009)
  * Engels  _et al._ [2015] T. Engels, D. Kolomenskiy, K. Schneider, and J. Sesterhenn, FLUSI: a novel parallel simulation tool for flapping insect flight using a Fourier method with volume penalization, SIAM J. Sci. Comp. 38, S03 (2015)
  * Griewank and Walther [2000] A. Griewank and A. Walther, Algorithm 799: revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation, ACM Trans. Math. Software (TOMS) 26, 19 (2000)
  * Blonigan and Wang [2012] P. Blonigan and Q. Wang, Least-squares shadowing for chaotic nonlinear dynamical systems, J. Comp. Phys. 34, 1 (2012)
  * Blonigan [2017] P. Blonigan, Adjoint sensitivity analysis of chaotic dynamical systems with non-intrusive least squares shadowing, J. Comp. Phys. 348, 803 (2017)
  * Chhabra and Richardson [1999] R. Chhabra and J. Richardson,  _Non-Newtonian flow in the process industries: fundamentals and engineering applications_ (Butterworth-Heinemann, 1999)
  * Peng  _et al._ [2014] Y. Peng, B. Lv, J. Yuan, H. Ji, L. Sun, and C. Dong, Application and prospect of the non-newtonian fluid in industrial field, in  _Materials Science Forum_ , Vol. 770 (Trans Tech Publ, 2014) pp. 396–401