{'fraction_non_alphanumeric': 0.026995122168519056, 'fraction_numerical': 0.006904692963895138, 'mean_word_length': 5.102759276879163, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mixing enhancement in binary fluids using optimised stirring strategies

Mixing of binary fluids by moving stirrers is a commonplace process in many industrial applications, where even modest improvements in mixing efficiency could translate into considerable power savings or enhanced product quality. We propose a gradient-based nonlinear optimization scheme to minimize the mix-norm of a passive scalar. The velocities of two cylindrical stirrers, moving on concentric circular paths inside a circular container, represent the control variables, and an iterative direct-adjoint algorithm is employed to arrive at enhanced mixing results. The associated stirring protocol is characterized by a complex interplay of vortical structures, generated and promoted by the stirrers' action. Full convergence of the optimization process requires constraints that penalize the acceleration of the moving bodies. Under these conditions, considerable mixing enhancement can be accomplished, even though an optimum cannot be guaranteed due to the non-convex nature of the optimization problem. Various challenges and extensions of our approach are discussed.


Introduction
The mixing of binary fluids -the process by which a heterogeneous mixture of two miscible fluids is manipulated into a homogeneous blend of uniform composition -is at the core of a great many industrial and technological applications. The food and beverage industry, as well as the consumer product industry abound with examples where multiple fluid components are mixed into a final product. Adhesives, sealants, cosmetics, inks and paints all consist of multiple ingredients that need to be mixed into their final state during a complex industrial process. Efficiency and consistency are paramount in maintaining a quality product that is cost-effective to manufacture. Some of the strictest tolerances in mixing quality can be found in the pharmaceutical industry where medication has to be mixed into precise doses. But also chemical engineering processes, such as polymer production, rely on accurate mixing to facilitate the proper chemical reactions and to reduce undesirable by-products (for an overview of theoretical and practical aspects of mixing, see Paul et al. (2003) or Uhl (2012)).
Mixing processes can be induced actively or passively. The active strategy is commonly based on a stirrer system, where paddles or rods agitate the binary mixture, induce vortical fluid structures and ultimately blend the initial ingredients. The geometry, path and speed of the stirrers have a great influence on the effectiveness and efficiency of the mixing process and are the subject of mixing optimisation. Passive systems, on the other hand, possess no moving parts, but instead rely on a complex baffle system inside an inflow-outflow device that mixes initially separated fluid components.
In this article, we will concentrate on an active stirrer system and develop a mathematical and computational framework for the formulation and solution of a constrained optimisation problem that yields favorable stirrer protocols for enhanced mixing results in binary fluid systems. Constraints stem from speed and path restrictions on the stirrers: stirrer systems are subject to mechanical and material limitations, and paddles or rods often cannot accelerate or change directions at will or too abruptly. In addition, while a significant part of industrial mixing processes involve non-Newtonian fluids, we will, for simplicity's sake, focus on Newtonian fluids. Furthermore, we will concentrate on inertial, laminar mixing. The inertial aspect of this parameter regime, described by a Reynolds number above the Stokes-flow regime, guarantees a rich and varied control space, taking advantage of advective, unsteady and diffusive processes, while the laminar aspect avoids divergences of the direct-adjoint optimisation scheme due to the existence of positive Lyapunov exponents linked to turbulent fluid motion. Despite these restrictions, a great many mixing processes fall into our chosen parameter regime.
Research in mixing has a long and remarkable history, covering theoretical aspects as well as technological applications. A large body of literature has been devoted to mixing in simplified fluid models, for example neglecting viscosity, surface tension, density differences or fluid inertia. The primary mechanism has been identified as streamline stretching (Spencer & Wiley 1951) where the interface between two fluids is repeatedly distorted and redistributed into the bulk of the mixing volume. Among these simplifications, Stokes mixing, i.e., the mixing of highly viscous fluids where inertial effects can be neglected, has arguably received the most attention. This tendency has been further fueled by the rise of micro-mixers where multiple fluid components are injected into a micro-device and extracted affter the mixing process is completed (see, e.g., Orsi et al. (2013); Galletti et al. (2015)).
More mathematical investigations studied the breakdown in scales and the statistical properties of the observed cascade of fluid filaments. Iterated maps have often been used to determine measures that describe the pertinent scale dynamics or to design optimal mixing results in these measures (Mathew et al. 2007;Gubanov & Cortelezzi 2010;Lin et al. 2011;Finn & Thiffeault 2011). Many of these findings can be found in Sturman et al. (2006). Of particular interest was the rise of chaotic mixing motion from pure advection, even for laminar flows (Aref 1984;Ottino 1989;Liu 2008).
Rather recently, the inertial, but laminar mixing regime has been explored using advances in optimisation techniques. These studies build on the definition of proper mixing measures (Mathew et al. 2005;Thiffeault 2012) and break with the focus on hydrodynamic instabilities (Balogh et al. 2005) to increase mixing. Using wall-mounted blowing/suction control in a channel, improved mixing could be accomplished by directly targeting a mixing measure, rather than a flow instability (Foures et al. 2014). Further studies (Vermach & Caulfield 2018;Marcotte & Caulfield 2018) have extended this approach to higher dimensions and stratified flows.
The present article will remain in the inertial, but laminar regime and accomplish mixing of a binary fluid by embedded stirrers. These stirrers are constrained to specific paths, but can move along them in a manner that enhances mixing over a user-specified time horizon. The mixedness of the binary fluid is quantified using mix-norms (Mathew et al. 2007;Thiffeault 2012) for a passive scalar; an objective based on this measure is then optimised by a nonlinear, gradient-based scheme, which in turn provides the associated stirring protocol.

Governing equations
The focus of this study is the mixing process of a binary, miscible and Newtonian fluid by multiple circular stirrers on prescribed paths, and its optimisation by manipulating the 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 θ. We have with u as the velocity vector, p as the pressure field, and θ as a passive scalar (ranging from zero in one fluid to one in the other). The governing equations have been expressed in non-dimensional form using a characteristic length L 0 and velocity scale u 0 . This choice introduces the Reynolds number Re and the Péclet number Pe, to express the kinematic viscosity and the diffusion coefficient of the mixing fluid in non-dimensional form. Furthermore, the system of equations (2.1) contains terms that model the embedded stirrers via a Brinkmann penalisation approach (see Angot et al. 1999). The multiple solid stirrers, indexed by the subscript k, are characterised by their velocity u s,k , and are taken as circular in cross-section. The masks χ k describing the embedded solid bodies equal one for points occupied by the k-th stirrer and zero outside of it. The mask χ accounts for the overall geometry, such as the domain boundaries. The constant C η ensures the rapid relaxation of the fluid variables towards the respective values imposed by the stirrers or the geometry. The above formalism allows the efficient treatment of objects moving through a background grid on which the motion of the surrounding fluid is describes.
Details of this approach and its numerical implementation can be found in Eggl & Schmid (2018). The setup shown in (2.1) imposes no-slip velocity boundary conditions on the stirrers and Neumann conditions for the passive scalar on the solid bodies.

Measuring mixedness
In anticipation of our stated goal of enhancing mixing efficiency, we have to introduce a measure that quantifies the degree of mixedness of a particular flow state. This measure shall be based solely on the passive scalar field θ.
In general, mixing is defined as the reduction of inhomogeneities of a given indicator field (Paul et al. 2003), which still leaves open a precise mathematical definition to be used in our case. Several norms of the passive scalar θ that attempt to mathematically define the measure of mixedness have been proposed and used in the past (Mathew et al. 2005), among them the variance or the more complex negative-index and fractional-index Sobolev norms (Thiffeault 2012;Foures et al. 2014). As the choice of norm may influence the outcome of the optimisation, but will not affect the design of our computational optimisation platform, we will focus on the Sobolev norm of negative fractional index of the passive scalar θ. A measure of this general type downplays the role of small scales and instead directs mixing efforts towards larger fluid elements. It attains higher values for an unmixed field (with high levels of inhomogeneities) and decreases as the scalar field becomes more mixed. Mathematically, the mixedness measure is given as with Ω denoting our computational domain, and |Ω| representing its size (volume or area). In the above definition, we have assumed, without loss of generality, a zero mean of the passive scalar field θ. Throughout this paper we will be optimising with respect to this quantity, but we stress again that other norms can be employed without conceptual changes in the optimisation procedures. The fractional exponent of 2/3 can be justified using arguments from optimal transport and ergodic theory. Examples of previous studies using mix-norm optimisation employed −1/2 (Foures et al. 2014) or −1 (Lin et al. 2011).

Mixing protocol
As a first attempt at optimising the mixing of a binary fluid, we will concentrate on a stationary circular vessel with two stirrers on circular paths or distinct radii (see figure 1 for a sketch of this configuration). The stirrers have a circular cross-section, and their velocities along their respective paths are undetermined and subject to optimisation and constraints. The circular path is conveniently defined in polar coordinates, while we formulate the remaining equations in Cartesian coordinates, and we thus introduce the vector-valued function l, which transforms between the two coordinate systems according to (l 1 (φ), l 2 (φ)) = (− sin(φ), cos(φ)) with φ as the angle traversed along the path of the circle. The parameterisation of the velocity of the k-th stirrer thus becomes where ϕ is the sum of the angles travelled along the path, i.e., ϕ k (t) = t 0 ω k (s)ds. (2.4) Following the notation of the governing equations, we define ω k (t) as the rotational speed of the k-th solid about the centre of the vessel, r k denotes the distance from the same centre, Ω k stands for the rotation about the stirrer's centre and a k represents the stirrer's (vectorial) radius. For simplicity, we take Ω k = 0; the alternative choice Ω k = ω k yielded very similar results in enhancing mixing efficiency.

Constrained optimisation
We can then state the optimisation problems as finding a time-dependent velocity protocol ω k (t) for each of the two stirrers such that the mix-norm of the passive scalar is minimised over a prescribed time horizon. This minimum has to be achieved while satisfying the governing equations and respecting constraints and bounds on the stirrer velocities. Mathematically we have and u s,lower u s,k u s,upper k = 1, 2, (2.5d) and a s,lower a s,k a s,upper k = 1, 2. (2.5e) The constraints on the stirrer strategy are threefold: the first constraint limits the L 2norm of u s , i.e., the kinetic energy of stirrers' motion along their paths, expended over the time horizon T to a maximum value of E 0 ; the second and third impose upper and lower bounds directly on the stirrers' velocities and accelerations, respectively. All restrictions could conceivably stem from mechanical limitations of the mixing apparatus. In our study, we will consider the constraints successively in order to determine the influence they impose on the optimisation results.

Computational framework
The implementation of the above optimisation problem requires the discretisation of the governing equations and the reformulation of the constrained problem (2.5) in terms of an unconstrained one.

Numerical scheme for the governing equations
Starting point for the numerical treatment of mixing enhancement is the open-source FLUSI software (Engels et al. 2015), in which the governing equations are discretised on a Cartesian, double-periodic domain for the two-dimensional case. This formulation allows the application of Fourier-spectral techniques to represent the spatial derivatives. The outer perimeter of the mixing vessel (with radius R) and the two stirrers on their respective circular paths (with radius r 1 and r 2 ) are described by a Brinkman penalisation technique as shown in (2.1). The original software has been augmented by the passive scalar field and embedded into a gradient-based optimisation formalism.
The Fourier-spectral discretisation allows the replacement of spatial derivatives with a multiplication by components of a wavenumber vector k = (k 1 , k 2 ) T , with the indices 1,2 indicating the two coordinate directions. Mathematically, we introduce this replacement where we introduced the Hadamard (element-wise) product • (see Horn & Johnson (2012)) and assumed implicit (Einstein) summation over identical indices. Particular care has to be exercised when evaluating the nonlinear terms, as aliasing errors can lead to inaccuracies and numerical instabilities. A low-pass Hou-Li filter (Hou & Li 2007) has been applied to avoid scattering of unresolved, small scales onto resolved, large scales. In addition, P3DFFT (Pekurovsky 2012), a highly efficient, parallel Fouriertransform library, is used to ensure scaling on parallel computer architectures.
Finally, the representation of solid bodies on an underlying Cartesian grid calls for a transfer of geometric information onto the background mesh. This transfer is accomplished by mollified delta-functions, smoothing the otherwise discontinuous mask onto the grid and thus avoiding numerical inaccuracies and instabilities (Kolomenskiy & Schneider 2009).

From constrained to unconstrained optimisation
A common reformulation of the constrained optimisation problem (2.5) as an unconstrained problem introduces Lagrange multipliers (or adjoint variables) for the dependent variables of equations (2.1): the adjoint velocity will enforce the momentum equation, the adjoint pressure the divergence condition, and the adjoint passive scalar the transport equation for θ. The augmented Lagrangian -consisting of the cost functional and the scalar product of the adjoint variables and the governing equations -then needs to minimised. A system of equations, referred to as the KKT-system, can then be derived by setting to zero the first variation of the augmented Lagrangian with respect to direct (original) and adjoint variables. The first variation with respect to the adjoint variables recovers the original set of governing equations. The first variation with respect to the original variables produces, after considerable algebra, a set of equations governing the adjoint variables. The first variation with respect to the control variables (in our case, the velocity strategies ω 1,2 (t)) will furnish the gradients of the cost functional with respect to ω 1,2 (t) which will be used to enhance the mixedness of our fluid system via improved stirring strategies.

Adjoint equations
Denoting the adjoint variables (velocity, pressure, passive scalar) by u † , p † and θ † , the governing equations for their evolution, in semi-discretized form, read It is important to note that the above adjoint equations are linear in the adjoint variables, but are dependent on the direct variables u i . Moreover, it should become apparent that the adjoint equations have to solved backwards in time, from t = T to t = 0. The optimality conditions, stemming from the first variation with respect to the control variables, result in the adjoint rotational velocity along the circular paths given by This value provides the gradient information in our iterative optimisation scheme and, together with a line-search routine, updates the current stirring protocol to a more effective one. The optimisation terminates when no more progress can be made, and the magnitude of the cost functional gradient drops below a prescribed threshold.

Implementing additional constraints
Additional constraints that need to be enforced are incorporated into the gradientbased optimisation routine. This is accomplished by projections and thresholding. In this case, the gradient -computed from the adjoint equations and the optimality condition, without imposed constraints -is projected and properly curtailed to comply with energy constraints and velocity bounds. Details of the numerical implementation of this procedure can be found in Eggl & Schmid (2018).

Summary of optimisation procedure
The full optimisation scheme then proceeds along the following lines. Starting with an initial guess of the stirring protocols ω 1,2 (t), we solve the governing equations (2.1) forward in time over a chosen time horizon [0, T ]. In a second step, the adjoint set of equations (3.2) are solved, starting with the proper terminal condition, backwards in time from t = T to t = 0. The adjoint variables are then used to evaluate the optimality condition and retrieve the gradient of the cost functional with respect to ω 1,2 . This gradient is then furnished to a standard optimisation routine (such as steepest descent or conjugate gradients) which, together with a line-search routine, produces a new and improved stirring protocal. In this last step, all constraints imposed on the stirring functions will be imposed by the aforementioned projections and thresholding. After this step, the next iteration is started. The optimisation terminates when no further progress can be made, within the constraints imposed on the system.
It is worth mentioning that additional complications arise from the fact that the governing equations (2.1) are nonlinear and, as a consequence, the adjoint equations (3.2), while linear in the adjoint variables, depend on the direct variables u, θ. This dependency requires the storage of direct fields during the simulation of (2.1) and their injection into (3.2), in reverse order, during the integration of the adjoint system. For efficiency reasons, this exchange between direct and adjoint simulations is handled by checkpointing, where we trade excessive storage requirements for an increased simulation time. The revolve library (Griewank & Walther 2000) accomplishes this task in an optimal manner.

Test cases and results
We follow the setup shown in figure 1 with two circular cylinders of radius a 1,2 = 1, moving on two concentric circular paths of radius r 1 = 3.5 and r 2 = 1.5 and embedded in a circular (stationary) vessel of radius R = 5. The Reynolds number and Péclet number are chosen as Re = P e = 1000.
A further parameter in the optimisation concerns the time horizons over which (i) the control strategy is applied, and (ii) over which gradient information is gathered. The former time interval determines the window given to the stirrers to be active mixers; after this window is passed, the motion of the stirrers will stop, and only the remaining inertia of the fluid and diffusion will contribute to further mixing. The latter time interval determines the amount of information extracted from the evolution process that is used to compute an enhanced stirring protocol (applied over the former time window). The control horizon may be chosen shorter than the information (predictive) horizon: in this case, a time-compressed strategy will be employed that accounts for, and optimises over, a more expansive time window. In our case, we will juxtapose a short-term strategy with T control = 1 and a longer-time strategy with T control = 8 and assess the optimised strategies in either case. Both protocols, however, have access to information over a temporal interval of T info = 8. Finally, the simulations have been continued over T sim = 32 units to track the further development of the instigated mixing processes; rest inertia and diffusion will remain the only mechanisms during this stage. A summary of our choice of parameters is sketched in figure 2.
Before proceeding to the various optimisation studies, it is instructive to reflect upon possible mixing mechanisms given the setup in figure 1. The most obvious strategy for mixing a binary fluid consists of a plunging motion, where the cylinders push through the initial interface, distort it and drag fluid one into regions occupied by fluid two, and vice versa (figures 3a-c). This type of strategy is nearly exclusively implemented in industrial mixers. Despite its omnipresence in applications, alternative strategies are often equally or more effective, foremost among them vortex shedding due to unsteady and abrupt motion of the stirrers, affably denoted as the vortex cannon strategy (figures 3d-f). In this case, the stirrer generates a sequence of startup and stopping vortices by rapid oscillations or abrupt directional changes along the circular paths. The shed vortices then act as effective autonomous mixers that, once they reach the initial or distorted interface, further deform the passive scalar field and locally (and globally) reduce the mix-norm. In this manner, a single stirrer can clone 'fluid stirrers' (shed vortices) and thus multiply its mixing effectiveness. In a further possible strategy, vortices can be generated in the fluid that collide with each other and thus generate filaments, which are then subjected to more rapid diffusion and homogenisation (figures 3g-i). Of course, this vortex collision strategy is strongly dependent on the initial condition of the passive scalar -and for this reason, less transferable to a general, realistic mixing strategy -, nonetheless, within our computational framework, it is a viable and pervasive strategy utilised by our directadjoint algorithm. A far more transferable mixing strategy is the collision of vortical structures with the outer wall (figures 3j-l) whereby a large fluid element is broken up into smaller elements which further interact with other vortices and are subject to increased diffusion due to the breakdown in scales. Finally, the embedded physical stirrers can themselves interplay with the vortical structures they generate, acting as obstructions in the path of vortices (figures 3m-o). A collision between a stirrer and a vortex will split the vortex and yield smaller scales, hence contributing to a decrease in the mix-norm. This final strategy will continue to cause a moderate breakdown in scales, even after the control window has closed and no more stirring motion is allowed.
Given these five fundamental strategies, illustrated in figure 3 with samples from our simulations, the direct-adjoint looping algorithm will select from and combine these options into a coherent strategy, given the chosen parameters and user-specified constraints.

Overview of test cases
We will consider six cases, grouped into three examples. Each example consists of a short-time strategy with a rather limited control horizon of T control = 1 and a longtime strategy with a more generous horizon of T control = 8. These two T control -settings will impose noticeable constraints on the choice of strategies, the interplay of dynamic processes and the feasibility of the final protocol. The three examples further distinguish themselves by the number of external constraints: starting with pure energy constraints, via energy and velocity constraints, to energy, velocity and acceleration constraints. Along this course of action, algorithmic requisites and physical requirements will be encountered and discussed.
Convergence of the iterative scheme is principally governed by constraints imposed on the optimisation problem. First, the nonlinear nature of the governing equations precludes a guarantee to converge towards a global minimum; only a local minimum can be expected. More importantly, additional constraints on the stirrers, such as energy, velocity or acceleration bounds, can convert a semi-norm to a full-norm optimisation problem. Semi-norm optimisation problems (Foures et al. 2012;Blumenthal et al. 2017) are 'open-ended' in the sense that the stirrer velocities increase without bounds, while the mixing measure steadily improves. In this case, the iterative optimisation scheme terminates when the adjoint variables, due to excessive direct velocities, no longer furnish useful gradient information for further improvement. In other words, the optimisation comes to an end when the signal-to-noise ratio for the gradient information drops to a value near unity, even before the cost-functional gradient attains a small value. In the full-norm case, convergence is achieved when the cost-functional gradient falls below a small, user-specified threshold, indicating that a local minimum has been reached. In brief, iterations are halted when either the cost-functional gradient falls below a small threshold or the signal-to-noise ratio of the adjoint gradient information reaches unitywhichever scenario comes first.
Since mixing enhancement based on complex stirring strategies is a highly dynamic process -based on a rapid sequence of abundant vortical features -, a set of static snapshots cannot do justice to the intricacies of an optimised mixing protocol. For this reason, we urge the reader to turn to the animations in the supplemental material.

Cases 1 and 2: optimisation under energy constraints
The first two cases follow a common procedure whereby the cost functional (mixnorm of the passive scalar) of the constrained optimisation is minimised, subject to a penalisation of the control energy that accomplishes this minimum. Since the stirrers' kinetic energy is a measure of effort that goes into the mixing process, we add a corresponding term to the pure mix-norm cost functional. As a consequence, the energy expended by the stirrers is bounded to a user-specified value.    Figure 4 displays the results of our optimisation, visualised by iso-contours of the passive scalar θ at selected time instances. The control horizon is T control = 1. We observed that the optimisation does not utilise the 'plunging' option, as the stirrers remain nearly at their initial position. Instead, the entire energy available to the stirrers is used up in a rapid start-and-stop motion which initially causes multiple small-scale shed vortices that distort the plane interface, collide into each other and the stirrers, and merge into larger-scale vortex structures which eventually achieve good mixing. It is important to stress that for the calculation of this short-time mixing strategy, information about the full dynamics up to T info = 8 has been incorporated into the optimisation. In other words, the consequences of the limited stirring protocol up to T info = 8 are known to the optimisation, and adjustments to the control strategy can be made that affect the vortex dynamics beyond its active control window. The evolution of the passive scalar between T info = 8 and T sim = 32, however, is neither designed nor recognized by the optimzation algorithm; it simply plays out according to the action taken during the control and optimisation windows. We include this further evolution to underline our choice of the mix-norm as the mixedness measure, whose optimisation produces the small scales that are subsequently diffused during this 'cool-off' window.
We conclude that the absence of any plunging option points at the suboptimality of this particular strategy in achieving an enhanced mixing process. It is thus not pursued as a viable option by the direct-adjoint optimisation technique. In the interpretation of these results, it may be tempting to conclude that a different initial placement of the cylinders -closer to the initial interface -would have resulted in strategies that included plunging. However, a simulation of the same case (not shown), with the two cylinders starting immersed in the initial interface, came to the same conclusion: while, by design, there is a small amount of plunging in this case, the vast majority of the mixnorm reduction has been accomplished by the shedding of start-and-stop vortices by a vigorous oscillatory motion of either cylinder and a subsequent collision of the generated vortices. The utilization of the stirrers' energy to shed small "vortical stirrers" is a better strategy than the distortion of the interface by simply moving through it with the stirrers.
The later part of the stirring strategy includes vortex collisions (see figure 4b), obstruction by the stirrers (see figures 4c,d) and collision with the outer wall (see figure 4e) to yield a well-mixed state at the end of the simulation horizon (figure 4f).
Increasing the control horizon from T control = 1 to T control = 8 leads to similar conclusions, even though the stirring action by the cylinders is less abrupt and jarring. Still, the bulk of the mixing action is achieved by shedding start-and-stop vortices which collide with themselves, secondary vortices and the wall to produce a mixed state in the end. Again, the absence of plunging is noteworthy. This is even more remarkable, as the increased control time horizon would certainly allow the stirrers to approach and reach the interface; yet, they remain close to their initial position.
In both cases, the strategy found by the direct-adjoint optimisation technique will yield increasingly larger velocities, as long as the integrated energy is constant. Eventually, the energy expenditure becomes more and more localized in time, with the stirrers barely moving. This optimisation route is a logical consequence of our current setup. It is closely connected to the semi-norm problem (see Foures et al. (2012); Blumenthal et al. (2017)): the mix-norm only contains the passive scalar θ, but does not account for the other dynamic variable, the velocities, in the optimisation. The energy of the stirrers is not sufficient to arrive at realistic stirring protocols that could be implemented in an experimental or industrial setting. To ensure applicability of our stirrer strategies to real-life settings, additional constraints are required.

Cases 3 and 4: optimisation under energy and velocity constraints
Following the findings of the previous section, in the next examples we limit the velocity of the stirrers to avoid excessive values of u s,i . This capping of the velocity is implemented by a projection of the raw cost-functional gradient onto control strategies that satisfy the given constraints (for details of this projection technique see Foures et al. (2014)). The resulting limit on the stirrer velocities provides a longer time window (up to the control horizon T control ) over which the specified energy can be expended. As a consequence, an extended and smoother movement of the stirrers is expected.
For the shorter control horizon T control = 1, figure 6 shows the outcome of our optimisation. The top stirrer starts by an oscillatory motion, creating start-and-stop vortices. The capping of its velocity, however, keeps the shed vortices within bounds; nonetheless, the optimality of the "vortex cannon" strategy can still be exploited. Both stirrers then move closer to the (already distorted) interface. But rather than plunging through it, they abruptly stop short of it and let the overtaking stop-vortices carry out the distortion of the interface and the subsequent mixing. Again, the optimisation algorithm selects the mixing by shed vortices over the plunging of the stirrers through the interface. The remaining mixing process is characterized by vortex collision (see figures 6a,d), collision with the wall (see figure 6c) and stirrer obstruction (see figures 6c,d).
Extending the horizon T control over which control is applied results in a change of strategy (see figure 7). The top stirrer now plunges through the interface -but not before stopping and starting on its circular path towards it. This uneven motion creates more vortical structures in the stirrer's wake that add to the sole plunging action of the stirrer itself. The result is a far more distorted interface (and consequently a lower mix-norm) than would be generated by a simple traversal. At the end of the motion, a back-and-forth motion is performed to generate, within the chosen energy and velocity constraints, additional shed vortices that further interact with the interface and other vortical elements. The second stirrer does not follow the strategy of the first. It engages in an oscillatory motion along its circular path and generates, as before, the resulting start-and-stop vortices that distort the interface and interact with the other vortices inside the container. Again, obstruction by the cylinders (see figures 7b,d) and vortex and wall collisions (see figures 7c,d) contribute to the continued mixing.
In both cases, a gentler stirring strategy is observed. However, the problem of converging towards a realistic mixing protocol has not been solved completely. While we (a) explicitly avoid highly localized action of the stirrers with excessive velocities, we now tend towards favoring strategies with excessive acceleration. In other words, within our efforts to limit the total expended energy of the stirrers while capping their velocities, the optimisation algorithms tends towards strategies that are characterized by large accelerations (high velocity gradients). This should not come as a surprise as the strength of shed vortices from the stirrers' unsteady motion is proportional to their acceleration. Our imposed constraints do account for energy and velocities, but not velocity gradients, of the stirrers. As a consequence, we can seed our binary mixture with vortical elements of nearly unlimited strength. Again, this divergence is related to the above-mentioned semi-norm problem: the velocity field of the binary mixture is not accounted for in the mix-norm, and thus the optimisation scheme can achieve high-energy fluid states by highly accelerating stirrers (even though the stirrers' energy and velocities are capped). To limit the velocity of the fluid, we have to limit the acceleration of the stirrers. Again, additional constraints are necessary.

Cases 5 and 6: optimisation under energy, velocity and acceleration constraints
For accomplishing enhanced mixing in binary fluids, the direct-adjoint optimisation technique makes heavy use of an acceleration-based strategy: shed vortices generated by the abrupt motion of the stirring cylinders are injected into both fluids, and their interactions with the interface, themselves and the container wall yield a low mix-norm. A limit on this acceleration will result in a limit on the velocities in either fluid component and thus provide the necessary restriction for a successful semi-norm optimisation. To this end, we augment our optimisation scheme by additional terms accounting for the stirrers' acceleration. This type of penalisation is common in deblurring of images where strong gradients are detected and encouraged. In our case, additional projections are used to enforce the acceleration constraints.
For the short-term control with T control = 1 (see figure 8), the optimal strategy now includes a plunging of the first cylinder through the interface, while the second cylinder continues in a straight manner towards the interface but stops short of it. The wake A longer control horizon of T control = 8 yields a more varied stirring protocol. The first cylinder makes a farther excursion, plunging through the interface (not without stopping to generate additional shed vortices close to the interface) before stopping close to the interface and shedding two stop vortices. The second cylinder first approaches the interface, ejects a stop vortex before reversing and stopping short of the interface with another stop vortex. The generated structures interact with themselves and the wall to break down the binary fluid into a homogeneous mixture, although of less homogeneity (larger mix-norm) than for the short-term strategy. This reduced homogeneity can be attributed more to the restricted velocity range (which is due to the constraint of an equal energy budget across both time horizons) than the larger time-horizon strategy. The lower velocity maximum, combined with the limitation on acceleration, impedes the same amount of vortex shedding than for the shorter time-horizon case. For this reason, mixing cannot be as efficient.
At no point during either optimisation has the energy, velocity or acceleration of the stirrers exceeded the specified limits. As a consequence, these latter strategies are amenable to implementation in an experimental or industrial setting.

Summary, conclusions and remaining challenges
A direct-adjoint optimisation methodology has been applied to the problem of mixing of binary fluids. A circular configuration with two embedded stirrers on circular paths has been chosen, and the velocities of these two stirrers over a user-specified time interval have been determined in an attempt to enhance the homogeneity of the binary mixture. The gradient-based optimisation is effective in finding stirring protocols that yield enhanced mixing results, but convergence towards an optimum crucially relies on imposing proper constraints on the iterative algorithm. Since the mixing efficiency is based on only one dynamic variable, the passive scalar θ, but disregards the velocity fields, additional external constraints have to be imposed to properly define a feasible optimum. These constraints have to enforce limitations on the encountered fluid velocities, which are forced by accelerating the stirrers. The accelerations, in turn, inject vortical structures into both fluids by unsteady, Stokes-layer-type shedding of vortices. Thus, by restricting the stirrers' maximum acceleration, we arrive at a properly stated optimisation problem and a convergent direct-adjoint algorithm.
Under these conditions, the optimal strategy that improves on the mixedness of the binary fluid utilises a combination of prototypical mixing techniques, consisting of plunging, unsteady vortex shedding, collisions between vortices and the wall, and obstructions by the stirrers. While unsteady vortex shedding is the key strategy for the unconstrained (or insufficiently constrained) case, a more balanced protocol ensues when excessive accelerations are increasingly penalised. Nonetheless, a rather counterintuitive optimal mixing strategy has been determined for short and longer time-horizons. As a general tendency, shorter control windows reach a lower mix-norm state, as the stirrer motion is more vigorous over a more limited horizon.
In conclusion, the above direct-adjoint approach to mixing of binary fluids -when combined with an efficient spectral simulation scheme, a Brinkman-type penalisation to accommodate moving bodies, and a systematic checkpointing technology -has proven an effective and robust tool to design stirrer strategies for the enhancement of mixing. After proper limitations on the stirrers' acceleration have been taken into account and a proper measure of mixedness has been defined, stirring strategies exploiting the full range of fluid processes induced by fluid-structure interactions are found, that suggest realisable modifications to commonly employed stirrer-induced mixing methods for industrial applications. The accomplished increase in mixing efficiency is summarised T control iterations θ mix,t=1 θ mix,t=8 θ mix,t=32 AccPen in table 1 which lists the mix-norm values for t = 1, t = 8 and t = 32. With an imposed acceleration penalisation, the short control horizon (T control = 1) gives markedly better results than a longer one (T control = 8), even though the difference is less pronounced after rest-inertia and diffusion set in. The number of iterations taken by the direct-adjoint algorithm is displayed as well; longer time horizons typically converge faster, owing to a less abrupt stirrer protocol. The short-horizon case with only energy penalisation (case 1, above) is included for comparison. Despite a successful increase of mixing efficiency via uncommon and unexpected strategies and despite corroborating and supporting the chosen computational approach, the present study also raises a number of challenges and shortcomings. Foremost among them is the fact that gradient-based optimisation with nonlinear partial differential equations as constraints can only assure a local optimum; a globally optimal solution requires an additional, and often prohibitive, methodology, such as simulated annealing or other variations of the same concept. While a global optimum may certainly be desirable, we point out that the improvements in mixing efficiency from a local solution (as shown in this study) would already have a respectable impact on mixing results due to its omnipresence in many industrial settings. In this sense, improvements suggested by locally optimal solutions would be most welcome.
A further challenge consists of additional constraint handling, imposed by mechanical restrictions on the stirrer motion. Penalisation methods or auxiliary projections imposed on the gradient information are conceivable to address this issue.
The length of the optimisation window T info places a bound on the overall algorithm. The adjoint part of the simulations computes the sensitivity of the output functional (mix norm) with respect to our control parameters (stirrer velocities along their paths). For increasingly large optimisation windows these sensitivities diverge due to the quasichaotic behavior of the direct problem. As a result, meaningful sensitivity to aid our optimisation will get overwhelmed by general sensitivity due to chaotic motion, and if meaningful sensitivity is lost, the optimisation algorithm will stagnate or even diverge. An advance to higher Reynolds or Péclet numbers will encounter similar issues. While techniques to overcome this predicament are currently being developed (Blonigan & Wang 2012), their cost-efficient application to complex systems, such as mixing, is still an open problem.
However, within the constraints of this study, with a physical problem this rich in possibilities and with a computational approach to match, there is an abundance of extensions and opportunities. Besides obvious explorations of other parameter combinations, the optimisation of the stirrers' shape is certainly within the capabilities of the computational framework; a preliminary study in this direction can be found in Eggl & Schmid (2019) using cycloids and trochoids as cross-sectional stirrer geometries. The path of the stirrers (in our study, concentric circles) can also be optimised; a collisionavoiding constraint may pose an additional challenge in this case. Furthermore, since the wall constitutes an important component in the breakup of vortices, an optimisation of wall motion or wall corrugation may further facilitate a more rapid breakdown in scales. With a view towards industrial applications, a non-Newtonian fluid model may be implemented. Finally, injection mixing (where the unmixed fluids are introduced into a mixing device or passive baffle system, and extracted when fully mixed) could be treated within the same direct-adjoint framework. These possible extensions, some of which will be reported in future efforts, attest to the flexibility and efficacy of the computational setup; marked enhancements in mixing are expected in the above cases.
Declaration of Interests. The authors report no conflict of interest.