Adjoint-based interfacial control of viscous drops

Abstract We develop a continuous adjoint formulation for controlling the deformation of clean, neutrally buoyant droplets in Stokes flow. The focus is on surface tension-driven flows where the interface is deformed by the local fluid velocity. We apply results from shape optimization to rigorously derive the optimality conditions for a range of interfacial problems. In the cases of interest, we make use of boundary integral methods as a natural choice for the numerical discretization of the flow variables. In the static case, our methodology is tested on a tracking-type cost functional and corresponds to classic shape optimization problems. We show agreement with black-box finite difference-based gradients and accurate minimization of the cost functionals. Finally, we demonstrate the methodology on the control of unsteady droplet deformation through external forcing.


Introduction
In this article we develop and apply an adjoint-based approach for controlling incompressible two-phase flows with sharp interfaces governed by the Stokes equations. The main model we consider is that of where α is an indicator function. However, recent work on similar systems with discontinuous coefficients, such as Pantz (2005) for a heat equation and Allaire, Jouve & van Goethem (2011) for linear elasticity, has shown that the bulk variables are not differentiable with respect to the domain, i.e. not shape differentiable. The present manuscript considers the restrictions (u ± , p ± ) separately and develops an adjoint method for the coupled system, instead of a single-equation system like the one-fluid model. The novelty of our contribution lies in extending and verifying the methodology proposed by Pantz (2005) to models similar to the above surface tension-driven two-phase flow. In particular, we derive the adjoint equations and adjoint-based gradient for static and quasi-static two-phase Stokes systems, which have been successfully used in many studies of droplet stability and breakup (see Stone & Leal 1989;Pozrikidis 1990;Hou, Lowengrub & Shelley 1994;Stone 1994;Lac & Homsy 2007;Ojala & Tornberg 2015) and provide a practical set-up for advancing the control of multiphase flows. Furthermore, the use of shape calculus allows for deriving an adjoint-based gradient independent of the representation of the interface. From the continuous adjoint perspective, the smoothness of the resulting flow and control parameters is also very important. For example, smooth solutions to the Stokes problem require a smooth curvature, but this is not a priori guaranteed by the optimization procedure. For the forward problem, several time-varying interfacial models have been recently studied in Pruss & Simonett (2016) and shown to be well-posed, paving the way to similar results for the adjoint problem. The smoothness requirement on the interface and curvature also poses problems numerically. For example, the widely used volume of fluid (known as VOF) method requires substantial effort to provide an accurate representation of the curvature (see Popinet 2009). Our numerical discretization relies on boundary element methods (see Pozrikidis 1992) to provide accurate representations of all quantities of interest. In particular, since the interface is expressed explicitly, we can compute accurate geometric quantities. The use of boundary element methods in shape optimization problems also has a long history, as found in classic monographs such as Aliabadi (2002) and in recent applications from Alouges, Desimone & Heltai (2011), Bonnet, Liu & Veerapaneni (2020) and others.
The current work can also provide a starting point for the adjoint-based optimization of higher Reynolds number multiphase flows governed by the incompressible Navier-Stokes equations. However, the flow nonlinearity and interface breakup will require changes to the formulation. Nonetheless, the present work can be viewed as a verification formulation as well as a predictive approach for droplets in the Stokes regime. Note that optimal control of multiphase Navier-Stokes systems has been approached in the past with various models. For example, in Feppon et al. (2019) the authors derive a coupled thermal fluid-structure model and in Deng, Liu & Wu (2017) and Garcke et al. (2018) a phase field model was used to describe the interface. Phase-field models have seen the most development in this area, due in part to the fact that solutions are quickly varying but continuous and allow for a standard approach through adjoint-based methods. In contrast, the current method requires the more complex shape calculus to handle the sharp discontinuities explicitly, but the resulting gradients are not quickly varying and the optimization is shown to be well behaved. Other multiphase applications include the optimal control of a two-phase Stefan problem with level set methods in Bernauer & Herzog (2011), free-surface Stokes problems in Repke, Marheineke & Pinnau (2011), Palacios, Alonso & Jameson (2012) and Sellier (2016), etc. To the authors' knowledge, no previous work has incorporated sharp interfaces and surface tension in a fluid optimal control problem.
The paper is organized as follows. We begin by introducing the required notation, fluid model and notions of boundary perturbations in § 2. In § 3, we define the optimization problems of interest and derive the adjoint equations and the gradients required for the optimization using gradient descent-type methods. The discretization details are presented in § 4. The discretization and numerical results are focused on a single droplet under axisymmetric assumptions. The presented results are then numerically verified in § 5 for several test problems. Finally, concluding remarks are given in § 6.

Two-phase Stokes system
We are interested in the quasi-static evolution of multiple clean droplets in an incompressible flow. The droplets experience no phase change, but are driven by a constant surface tension at the fluid-fluid interface. Such a system is governed by the two-phase Stokes equations, which have been described in detail in classic monographs (see Ladyzhenskaya 1963;Pozrikidis 1992;. The notation and main auxiliary results used below are presented in appendix A.
We consider a three-dimensional domain Ω − ⊂ R d , d = 3, and its complement Ω + = R d \Ω − to denote the union of the droplets' interiors and the ambient fluid, respectively. The boundary between the two domains is Σ Ω − ∩Ω + . In the following, we assume Σ is a finite set of disjoint, closed, bounded and orientable surfaces. Furthermore, we require that Σ is of class C 2 to allow for the definition of a unique curvature at every point.
The droplets and the ambient fluid have dynamic viscosities μ − and μ + , respectively. The flow in each phase is described by the velocity fields u ± : Ω ± → R d and the pressures p ± : Ω ± → R. We also define the Cauchy stress tensor and the rate-of-strain tensor by σ ± −p ± I + 2μ ± ε ± , ε ± 1 2 (∇u ± + (∇u ± ) T ).
(2.1) Following these definitions, the static incompressible two-phase Stokes equations are given by The boundary conditions at the interface between the fluids are given as a set of jump conditions on Σ. Namely, the velocity must satisfy a no-slip condition and the surface stresses are related to the additive curvature of the interface by the Young-Laplace law. We write where f ± n · σ ± is the surface traction, γ is the surface tension coefficient and n is the exterior normal vector to Ω − . The jump operator is defined as It is convenient to non-dimensionalize the equations to reduce the number of parameters in the problem. We choose the following non-dimensionalizations, common in the literature: where U is a characteristic velocity magnitude and R is a characteristic droplet radius. The resulting non-dimensional parameters governing the systems are the viscosity ratio λ and the capillary number Ca, defined as In the remainder, we will continue by using the non-dimensional equations. We will write the non-dimensional stress jump condition (2.3) as follows: In order to evolve the shape of the droplets in time, we use a quasi-static approach in which we first compute the velocity using (2.2) and then displace the interface using an ordinary differential equation. The evolution of the interface is described by where X is a parameterization of the interface and V is a given source term. Typical examples of the source term are V (u, X ) (u · n)n + (I − n ⊗ n)w, (2.9) where w can be chosen to be 0, u or a different tangential correction to u. The choice of w is left open, since only the normal component of the velocity field governs the deformation of the interface. The tangential component can be artificially imposed for numerical reasons, e.g. to cluster points in regions of high curvature; for example, see the works of Hou et al. (1994) and Ojala & Tornberg (2015).

Shape derivatives
Interfacial control in multiphase flows with sharp interfaces is inherently a form of shape optimization. With this in mind, we will review some of the concepts behind shape optimization and perturbations with respect to the domain. We consider here Hadamard's method for boundary variations, which has been detailed rigorously in works such as Moubachir & Zolésio (2006), Allaire & Schoenauer (2007), Delfour & Zolésio (2001) and, from the view of differential geometry, in Walker (2015).
The method is based on defining so-called perturbations of the identity by where Id is the identity mapping,X ∈ W 2,∞ (R d ; R d ) is a perturbation and > 0 is a real constant. If is sufficiently small, then X is a diffeomorphism on R d (Allaire & Schoenauer 2007, lemma 6.13). As such, by Hadamard's method we can identify every admissible shape Ω with a perturbation X . Therefore, the method allows, among other things, defining differentiation of shapes by differentiation in the Banach space W 2,∞ (R d ; R d ) in terms of well known concepts, such as the Gâteaux and Fréchet derivatives.
In fact, using this simple definition, we can define the Fréchet derivatives of relevant functionals of the type We refer to Allaire & Schoenauer (2007) and Walker (2015) for proofs of these results and additional shape derivatives of quantities of interest, such as the normal and the additive curvature. The lemma below reproduces the main derivatives of functionals used in the current work.
LEMMA 2.1 ( Walker 2015, lemma 5.7). LetX ∈ W 2,∞ (R d ; R d ) be a sufficiently small perturbation, Ω ⊂ R d and Σ ⊂ R d a d − 1 manifold without boundary. Then, we have that the shape derivatives of the functionals (2.11) are Extensions to moving domains Ω(t) and Σ(t) are presented in Moubachir & Zolésio (2006). An important result in shape optimization is the Hadamard-Zolésio structure theorem from Delfour & Zolésio (2001, theorem 3.6, chapter 9), which guarantees that we can always express the first-order shape derivative of functionals (2.11) as boundary integrals that only depend on the normal component ofX . Restricting the shape gradient to the boundary of a domain Ω is of special interest to us, since it greatly simplifies the numerical aspects. In particular, it allows restricting the whole problem to a problem on the interface Σ between the droplets and the surrounding fluid.
Finally, we will informally consider the shape differentiability of surface functionals that depend on the problem variables p ± , u ± and X . The main issue, also discussed in Pantz (2005) and Allaire et al. (2011), is that u as a function on Ω is not shape differentiable, but u ± as restrictions to Ω ± are shape differentiable.
Recall that, for Ca < ∞, there is always a non-zero jump in pressure We can see from lemma 2.1 that the integrand must be differentiable up to the boundary to allow defining ∇j · n. In the following, we will restrict to differentiable functionals of the type i.e. that only depend on the normal component of the velocity field and the surface parameterization, both of which have continuous normal gradients. In the case of volume integrals, we would not need further restrictions on the cost functionals we consider.

Optimization problem
We now define two constrained optimization problems that we will attempt to characterize. In particular, we will use the Lagrange multiplier theory to determine adjoint equations and derive an expression for first-order variations of the cost functionals of interest. Firstly, we define a classic shape optimization problem where U 1 is the set of all admissible parameterizations of interfaces Σ ⊂ R d and E 1 corresponds to the constraints (2.2) and (2.3). In general, we require that the fluid-fluid interfaces remain of class C 2 .
Remark 3.1. Note that the parameterizations of a surface Σ are not unique. To rigorously define the set U 1 , we would be required to define an appropriate notion of equivalence classes (Delfour & Zolésio 2001).
Secondly, we consider the optimization with respect to a parameter of the problem, namely the capillary number Ca. This problem is defined by where U 2 R + and E 2 encodes the constraints (2.2), (2.3) and (2.8). At least conceptually, we can see that the second problem is a superset of the first one, since it simply includes an additional constraint. For this reason, we will focus our attention on the first problem, with extensions to the second problem as necessary.
The two optimization problems use tracking-type cost functionals restricted to the interface Σ. They are where u d and X d represent a desired surface normal velocity field and geometry. We differentiate between J k , which is a function of the control only, and J k , which is a function of the control and the state variables (as required). The equivalence between the two, through the corresponding constraints E k , is given by the implicit function theorem under some assumptions. The main theorem that is used in the theory of equality constrained optimization is stated below.
A proof of the Lagrange multiplier theorem above can be found in Luenberger (1997). For a control problem, we have that the variable x ( y, u), where y are the state variables and u is the control. In this case, the statement that (x 0 , z * 0 ) is a stationary point of the Lagrangian L can be split into the following equations: which are valid for all directions h in the necessary spaces. In the following we will show the specific Lagrangian functionals for our problems of interest and the equations satisfied by the adjoint variables.
Remark. Rigorous proofs of existence and uniqueness of solutions to (P1) and (P2) are out of scope for this work. We refer to the seminal work of Ladyzhenskaya (1963) on smooth domains for results related to the static problem (2.2) and (2.3). Extensions to the unsteady problem can be found in Denisova & Solonnikov (1991) and subsequent work.
3.1. Static problem We start by looking at the optimization problem (P1) with the static constraints (2.2) and (2.3). The resulting Lagrangian can be written as where z * are the adjoint variables and w ± (u ± , p ± ) are the state variables on Ω ± , such that w (w + , w − ). Expanding the constraints, we write where f λ denotes a weighted average, defined as The weak form of the Stokes equations, (2.2), can be obtained by integration by parts and then applying the following equality: Remark 3.3. The main subtlety in our derivation comes from the choice in constructing the Lagrangian. In classic single-phase works, such as Bewley, Moin & Temam (2001) and Wei & Freund (2006), the strong form of the equations is used with an integral over the whole domain. In the context of multiphase flows, such an approach would be possible using the one-fluid model, where (u, p) are considered as bulk variables with the jump conditions carried over as singular source terms. However, such a choice results in an incorrect shape derivative. As discussed in Pantz (2005), the reason for this is that the state variables are not shape differentiable over Ω. Therefore, equations with discontinuous coefficients and, in particular, multiphase flows must be treated with care. The correct construction involves considering the restrictions to Ω ± , i.e. (u ± , p ± ) as above, and enforcing values that are well-defined on the interface. For example, an alternative expansion of would lead to ill-defined results. The choice of a weak form in (3.6) is not formally required, but allows us to directly apply the results from lemma 2.1 and is consistent with the methods used in the shape optimization community, e.g. in Allaire & Schoenauer (2007).
We start by stating the equations satisfied by the adjoint variables in the following result.
LEMMA 3.4. At a stationary point of the Lagrangian (3.6), the adjoint variables w * ± satisfy the following equations: (3.11) The adjoint Cauchy stress tensor, rate-of-strain tensor and traction vector are in complete analogy to the definitions of the primal Stokes problem given in § 2. For the cost functional Proof . The adjoint equations are easily obtained by differentiating the Lagrangian with respect to the state variables w ± . We have that by algebraic manipulations. This is the weak formulation of the adjoint equations from lemma 3.4 with test functions (ũ ± ,p ± ).
We are now in a position to give the definition of the shape derivative of the cost functional from (P1). It is stated as follows.
THEOREM 3.5. We assume that (u ± , p ± ) and J 1 are shape differentiable. Then, the shape gradient of the cost functional (P1) is For the cost functional (3.1), we have that

17)
where j u · n − u d . The operators κ * and n * correspond to the adjoint operators of the Eulerian shape derivative of the additive curvature and normal vector, respectively. They are defined in appendix B.
Proof . A detailed derivation of the shape gradient is given in appendix C. We note that both the derivation of the shape gradient in theorem 3.5 and the adjoint equations in lemma 3.4 are purely formal; they assume that the state variables have shape derivatives of the required regularity.

3.2.
Quasi-static problem The extension from (P1) to the unsteady problem (P2) is straightforward from the point of view of the Lagrangian formalism. However, introducing the unsteady element complicates the theoretical analysis of the problem, since it is not clear if the solution will maintain the required smoothness. Similar issues appear at a numerical level, as we will see in subsequent sections. Furthermore, the shape derivative formalism from § 2.2 needs to be extended. Intuitively, we now want to perturb the time-dependent domains Ω ± (t). In the (d + 1)-dimensional (t, x) space, one can think of a fixed domain Ω as a cylinder, while a moving domain Ω(t) will be a general tubular manifold. Domain perturbations of this kind were made rigorous in Moubachir & Zolésio (2006), to which we refer to any additional details.
As before, the Lagrangian functional associated with (P2) is given by 18) or, explicitly, LEMMA 3.6. We assume that w ± and J 2 are shape differentiable. Furthermore, the source term V in (2.8) is also assumed to be shape differentiable. Then, at a stationary point of the Lagrangian (3.19), the adjoint variables w * ± satisfy the Stokes equations (2.2) with the jump conditions where X * is given by lemma 3.7. In the case of the cost functional from (3.1), we have that Also, using the forcing term described in (2.9), we have that Proof . The proof is completely equivalent to that of lemma 3.4.
LEMMA 3.7. Under the assumptions of lemma 3.6, we have that at a stationary point of the Lagrangian (3.19), the adjoint variable X * satisfies the following evolution equation: where C * and S * are defined in lemma 3.4 and V * is the adjoint operator corresponding to (2.8). In the case of the cost functional (3.1), we have that Similarly, for the choice of forcing term described in (2.9), we have that Proof . The derivation of the above evolution equation is presented in appendix C.
Determining the shape gradient of the cost functional (3.1) is a direct application of lemma 2.1.
THEOREM 3.8. The gradient of the cost functional for (P2) is given by Proof . The proof of theorem 3.8 follows from the definition of the Lagrangian given in (3.19). Note that the gradient is simply a classic derivative in R in this case.
This completes the statement of the two optimization problems we have proposed. In both cases, we have access to first-order variations that can be used as part of gradient descent algorithms to perform the optimization.

Numerical methods
In this section, we will detail the numerical methods used to solve the Stokes system (2.2) and (2.3) and the adjoint Stokes system described in lemma 3.4. Since the Stokes equations are symmetric, we require a single solver in both cases. We will also look into discretizations of the shape gradient in theorem 3.5, the interface evolution equation (2.8) and the corresponding adjoint evolution equation from lemma 3.7. The numerical methods described here have been implemented in an open-source code that can be found in Fikl (2020).

Boundary integral equations
An efficient method of solving the Stokes equation is based on boundary integral equations (Pozrikidis 1992). In our case, this is particularly useful because the entire problem reduces to the interface if the cost functional is expressed on Σ alone, as is the case for the choices presented in (3.1).
Boundary integral methods are based on Green's function fundamental solutions to the Stokes problem. Using the fundamental solutions as kernels in a linear integral operator, there are many possible representations of the two-phase Stokes problem. A common representation is based on the Lorentz reciprocal theorem and has been used in studies of interfacial dynamics (see Pozrikidis 1990;Stone 1994;Lac & Homsy 2007). However, this representation is less suitable for our case, since the pressure and stress kernels are hypersingular. Hypersingular integrals are generally harder to analyse and accurately evaluate numerically.
A second representation for two-phase flows with a stress discontinuity has been discussed in Pozrikidis (1990) and is described in detail in Pozrikidis (1992, chapter 5.3). It is a so-called single-layer potential representation, with the velocity field described by, where q : Σ → R d and G ij are usually referred to as the density and the Stokeslet, (4.4). This representation is complemented by definitions of the pressure and stress, for x / ∈ Σ, where both have been appropriately non-dimensionalized, as described in § 2.1. The usual summation convention over repeated indices has been used in the definitions. The fundamental solutions for the three main variables in Stokes flow are given by where r x − y and r r . To solve for the density we use the provided boundary and jump conditions, which, however, require limits as x → Σ of the layer potentials above. These limits are generally not trivial, since the kernels in question are all singular at Σ. The jumps in the velocity, pressure and traction are well known and given in Pozrikidis (1992, § 4.1). Using the surface limits of the traction and the jump (2.3), we obtain the following equation (Pozrikidis 1992, (5.3.9)): In the non-homogeneous far-field boundary conditions case, additional terms can be added, as described in Pozrikidis (1992). Equation (4.5) is a Fredholm integral equation of the second kind. The linear operator acting on the density q is generally well-conditioned. In fact, it matches the conditioning of the underlying physical problem, which is well-conditioned for 0 < λ < ∞. This allows for efficient results by classic iterative linear solvers such as GMRES.
The layer potentials presented so far are sufficient to solve the forward Stokes equations (2.2), (2.3) and (2.8). However, in the case of the adjoint equations, we further require knowledge of the normal and tangential velocity gradients. They can be obtained by differentiating the velocity (4.1) and are given, for x / ∈ Σ, by The kernel for the velocity gradient is given by The surface limits of the new layer-potential operators are not standard, but they can be easily derived from the main variables and the definition of the stress tensor (2.1).
THEOREM 4.1. The jump conditions for the normal and tangential velocity gradient layer potentials, defined in (4.6)-(4.7), are given by where p.v. denotes the Cauchy principal value interpretation of the singular integral.
We must also obtain formulae for the stress and strain-rate tensors that appear in theorem 3.5. The stress tensor can most easily be obtained algebraically from the traction, the pressure and the components of the velocity gradient from (4.6) and (4.7), by means of (2.1). Then, the strain rate tensor can be directly obtained from (2.1). Explicit formulae for the stress components in both two and three dimensions are given in Aliabadi (2002, § 2.5.1).
Finally, for the numerical implementation, we will make an axisymmetric assumption on the problem set-up, with coordinates x = (x, ρ, φ) (see figure 1). With this assumption, we must also change the meaning of the fundamental solutions we have used thus far. For example, in the three-dimensional case, G ij represents the interactions between a source point y and a target point x of strength q. However, in the axisymmetric setting, we must consider the interactions between a target point x and a ring of source points y for fixed (x 0 , ρ 0 ). This is achieved by analytically integrating the kernels in the azimuthal direction φ, with the aid of elliptic integrals. Several axisymmetric kernels have been derived in Pozrikidis (1992). We give in appendix D a short description of all the kernels required in this work.

Boundary element methods
The numerical discretization of the layer potentials from § 4.1 has been performed using a standard collocation method (also known as a boundary element method), similar to the methods described in Kress (1989, chapter 13.3). We will only present the details for a single droplet, as seen in figure 1. In an axisymmetric setting, all the integrals over the surface Σ can be reduced to integrals over the top half-plane and a one-dimensional interface Γ {x ∈ Σ : φ = 0}. We refer to Kress (1989), Aliabadi (2002) or Pozrikidis (1992) for additional details on the construction of one-dimensional collocation methods and only briefly describe the chosen notation for this work.
We define a uniform one-dimensional computational grid with parameter ξ ∈ [0, 1/2] discretized into M elements Γ m [ξ m , ξ m+1 ), for m ∈ {0, . . . , M}. On each element Γ m , we define a set of N + 1 uniformly distributed collocation points ξ m,n , for a total of N c (M × N) + 1 unique points. For a high-order finite element type discretization, we use the standard nodal Lagrange basis functions φ n . Finally, we use the standard Gauss-Legendre quadrature rule on regular elements and specialized methods for elements containing singularities.

Singularity handling and singular quadrature rules
The main difficulty in implementing a boundary element method lies in the handling of the singularities for each kernel. As seen in the previous sections, the adjoint problem requires all the quantities present in the Stokes problem with several different singularity types. The surface limit of the velocity layer potential (4.1) is weakly singular, with a log-type singularity, and the remaining layer potentials, e.g. (4.3), are defined by means of the Cauchy principal value and are strongly singular, with a r −1 -type singularity.
There are multiple methods used to handle singular integrals, such as changes of variables, singularity subtraction, regularization or generalized quadrature rules (see Aliabadi 2002, chapter 11). We focus on making use of generalized quadrature rules and analytical techniques for accurate numerical integration. For all singularity types we rely on an element subdivision method, as seen in figure 2.
This method proceeds by isolating the element containing the currently evaluated target point x i and splitting it into two elements with specially constructed quadrature rules on each. This allows reducing the problem to that of endpoint singularities, and minimizing the required quadrature rules. This is important because in the case of generalized quadrature rules, their construction depends on the location of the singularity, e.g. Carley (2006) and Kolm & Rokhlin (2001). There are two special cases that need to be handled: element endpoints consider the two existing adjacent elements as the subdivision and the two points at the poles (ρ = 0) only make use of the interior element for a one-sided evaluation.

Singularity subtraction
The pressure kernel (4.2), the traction kernel (4.3) and the normal velocity gradient kernel (4.6) can all be handled by the method of singularity subtraction. The singularity subtraction for the traction is described in Pozrikidis (1990). In the case of the pressure kernel, we make use of the well known identities Σ p j (x, y)n j ( y) dS = 4π and Σ ε ijk p j (x, y)n k (x) dS = 0, (4.10) where ε ijk is the Levi-Civita symbol forming the cross product. We can express the pressure layer potential (4.2), for x ∈ Σ, as (4.11) The new kernelp j is given bỹ (4.12) The axisymmetric form of this kernel can be found in appendix D. Finally, the layer-potential representation of the normal velocity gradient (4.6) can be desingularized by noting that and using the formulae derived for the pressure and traction layer potentials. Since the integrals are no longer singular, we use the standard Gauss-Legendre quadrature rules on each subdivision. The subdivisions are still required in this case, even though the singularity is removed, because the resulting integrand is generally not smooth (only continuous) and can degrade the accuracy of the solution.

Weakly singular integrals
For the weakly singular integrals in our problem, such as the velocity single-layer potential representation (4.1), we make use of the known Alpert quadrature rules. Alpert quadrature rules are a set of corrected trapezoidal rules that achieve high accuracy for log singularities. These quadrature rules can be constructed for endpoint singularities, as described in Alpert (1999), and require no additional changes to the kernels.

Strongly singular integrals
The remaining layer potential representation that requires special handling is that of the tangential velocity gradient, (4.7). This kernel is strongly singular and does not have known solutions for singularity subtraction. Therefore, we require a set of generalized quadrature rules for Cauchy principal value integrals. We have made use of the quadrature rules described in Carley (2006), based on the work of Kolm & Rokhlin (2001). However, we cannot make use of the subdivision method in this case. Applying the quadrature rules relies on a change of variables from the reference element [−1, 1] to the surface element, which cannot be naively applied for strongly singular integrals with endpoint singularities (Monegato 1994). Instead, we construct unique quadrature rules for each interior collocation point in the reference element. As before, the element endpoints are handled by considering the union of the adjacent elements with a singularity at the origin ξ = 0. The two poles need to be handled separately, since they are by definition one-sided and do not allow a change of variables. Therefore, we are forced to construct quadrature rules directly on the elements Γ 1 and Γ M . This is done by extending the results from Carley (2006) to arbitrary [a, b] intervals. Unlike the previous cases, we require N + 1 unique quadrature rules to handle the tangential velocity layer potential (4.7).

Geometry representation
The geometry is represented spectrally by Fourier modes, for high accuracy. Together with the discretization of the layer potentials from § 4.2, the result is a superparametric collocation method.
To compute the tangent vector, the normal vector and the curvature, we make use of fast Fourier transform. The fast Fourier transform is computed on the uniform grid formed by the collocation points, to avoid issues with non-uniform spacing. Since we only have access to the top half-plane, we mirror the geometry to form a closed curve on which the Fourier transform can be directly applied. Finally, all the quantities are interpolated to the required quadrature points using the Fourier coefficients computed at the uniform collocation points.

Evolution equation
The interface evolution (2.8) is integrated using the classic third-order strong-stabilitypreserving Runge-Kutta method At each stage of the Runge-Kutta method, we solve the discrete Fredholm integral, equation (4.5), and compute the velocity with (4.1). The time step used is based on restrictions by the capillary number from Denner & van Wachem (2015) and is given by is the smallest panel size in physical space. To maintain a smooth surface, we also perform smoothing by the penalized least squares method described in Eilers (2003) at regular time intervals. In the studied cases, the interface deformation is sufficiently small and does not require advanced mesh adaptation techniques. We refer to Pozrikidis (1990) and Ojala & Tornberg (2015) for a discussion on existing methods. 4.6. Discrete shape gradient The discretization of the shape gradient in theorem 3.5 will be performed using finite element-based methods. The collocation points and Lagrange polynomial basis from the boundary element method are used for this purpose. The inner product used in the discretization is given by (4.17) Discretely, this inner product is described using the mass matrix as follows: where the mass matrix M is given by for Gauss-Legendre quadrature nodes and weights (η p , w p ). We obtain weak forms of the adjoint normal operator n * and the adjoint curvature operator κ * from theorem 3.5 by a direct application of integration by parts. This helps avoid even higher-order derivatives of the geometry and reduces the required smoothness. Taking the results from theorem 3.5, multiplying by the basis functions and integrating by parts as required, we obtain a linear system of the type where g is the gradient and b is the right-hand side. To enforce a smooth boundary without cusps at the poles, we additionally require that dg dξ (0) = dg dξ (1/2) = 0. (4.21) In our work, we use Lagrange multipliers to impose the boundary condition and augment the mass matrix accordingly. Note, however, that this condition should be automatically satisfied if the flow quantities are similarly smooth in theorem 3.5.

Discrete adjoint quasi-static equation
The adjoint evolution equation in lemma 3.7 is also discretized by a finite element method. Finite element methods on moving meshes have been described in Dziuk & Elliott (2007) and we employ the same ideas to discretize (3.23). By integrating and using the surface Reynolds transport theorem from A.3, we obtain where V * now denotes the combined right-hand side of the evolution equation in lemma 3.7. As the terms are the same as those present in the shape gradient from § 4.6, the same steps are used to obtain an appropriate weak form. To maintain dual consistency at the time advancement level, we use the discrete adjoint of (4.14), namely whereV * M −1 V * . As the adjoint to (4.14), this method maintains the same stability region, but an analysis of its order conditions reveals that it is only second-order accurate. The gradient from theorem 3.8 is integrated in time using the same method, by noting that solving dg dt = − 1 Ca 2 Σ(t) κu * · n dS (4.24) with g(T) = 0, results in ∇J 2 = −g(0).

Optimization algorithm
The optimization method used in this work is based on gradient descent and requires first-order shape derivatives. For the problem (P1) they are given in theorem 3.5 and for problem (P2) they are given in theorem 3.8. The results provided in the respective theorems refer to the normal component of the gradient. Therefore, a descent direction for the directional derivative for a positive constant α. A similar reasoning holds for (P2), where the problem is greatly simplified by the fact that ∇ Ca J 2 ∈ R. This leads us to the gradient descent method described in algorithm 1.
In algorithm 1, we have directly used the gradient as a descent direction, resulting in a steepest descent method. This choice is not always optimal or can lead to a jamming phenomenon. A survey of gradient descent methods is given in Hager & Zhang (2006). In our work, we have used a conjugate gradient method based on the Fletcher-Reeves update with a step size α (n) chosen to satisfy the Armijo condition, but not the full Wolfe conditions.
It is well known in shape optimization problems, like (P1), that the surface shape gradient can have less regularity than the original solution. Intuitively, we can see this from the appearance of the normal velocity gradient or the derivatives of curvature in the shape gradient from theorem 3.5. This lack of smoothness can result in diverging results when the gradient is used repeatedly in a gradient descent method. To counter some of these issues we perform additional post-processing steps on the geometry X (n+1) after the update in algorithm 1. They are as follows.
(1) The geometry is smoothed by a penalized least squares method, described in Eilers (2003), where the penalization parameter is chosen at runtime. (2) Since X (n) describes only the top half-plane of the axisymmetric geometry, we must maintain ρ (n+1) > 0. This is ensured by a choosing a descent step α (n+1) that is sufficiently small. (3) We also remove any translations by enforcing that the on-axis points satisfy (4.27) With these additional steps, we have observed that the optimization procedure is sufficiently robust. In particular, accuracy in the computation of the required layer potentials in § 4.1 is maintained and the axisymmetric assumptions are enforced throughout. For more general applications, it may be desired to add additional regularization terms to the cost functional. For example, we can consider to penalize solutions with a non-smooth curvature. We can also use an H 1 inner product, which would require solving an additional Laplace-Beltrami equation on the interface.
Remark 4.2. The proposed algorithm for the shape optimization problem, (P1), does not take into account the structure of shape spaces (Michor & Mumford 2006). Since shape spaces are not flat, a general update formula can be written as where T is a parallel transport operator. Taking the additional structure into account can lead to better behaved algorithms, as seen in Ring & Wirth (2012) and Schulz (2014 The optimization problem defined by (P1) is simpler in this regard, since it takes place in R + and relies on the smoothness of the evolution equation (2.8). Here, we can use the standard steepest descent method with no additional post processing.

Numerical results
In this section, we present our results for the optimization of the two problems we have set out in § 3, namely the shape optimization problem, (P1), and the parameter optimization problem, (P2). All the presented results will be for a single axisymmetric droplet. The formulation of the adjoint equation and gradients from § 3 are general and can be applied to one-, two-or three-dimensional systems. The numerical methods described in § 4 can also be generalized to multiple droplets, with care needed for near-singular evaluations of the layer potentials.
In the following, we will present results for an (N = 4)-order collocation scheme. We have chosen this order as a baseline because the adjoint-based gradient from theorem 3.5 contains a second derivative of the velocity field (in the form of κ * ), so the result will always be of order N − 1 due to the weak formulation. We maintain the same discretization parameters for all the tests we present, as detailed in table 1. The different types of quadrature rules are given in § 4.3.
The errors are computed with the aid of the discrete L 2 norm, defined as where M is the mass matrix. The reported errors are relative, i.e.
where g is a known exact solution andĝ is the approximation. In computing orders of convergence, we use the largest element size as an abscissa, i.e.

Adjoint gradient
We analyse here the validity of the adjoint-based gradients from theorem 3.5 and theorem 3.8. In both cases, the gradients will be validated by comparison against finite different approximations. Further validation is presented in the next section by finding optimal solutions to known problems.

Test 1: static problem gradient validation
We start by validating the adjoint gradient from theorem 3.5 by a standard first-order finite difference approximation where (·, ·) is the inner product (4.17). To ensure that the perturbed interface X = X + i h i remains smooth and satisfies the axisymmetric assumptions for every i, we choose To compare against the adjoint gradient, we also dot the adjoint gradient with the bump functions in the same way, i.e. (g AD ) i = (g, h i ) with g given by (4.20).
We present results for panel sizes M ∈ {8, 16, 24, 32}, finite difference step sizes ∈ {10 −d : d ∈ [[1, 7]]} and a finite difference bump thickness σ 2 = 10 −2 . For the test set-up, we use u d ≡ 0 in (P1) with Ca = 1 and λ = 0.2 in an extensional flow where the geometry is given by an ellipse, for (a, b) = (1, 2), parameterized as x(ξ ) = r(ξ )(cos(2πξ), sin(2πξ )), The relative errors and gradients can be seen in figure 3. We expect that for a fixed and sufficiently fine mesh, we retrieve the first-order convergence in . We also present the exact values of the errors for = 10 −6 in table 2. For a fixed finite difference step size, we retrieve the expected discretization order (N − 1) or 3. We can also see that the errors plateau at ≈ , when the finite difference errors dominates.

Test 2: quasi-static problem gradient validation
To validate the adjoint gradient of the quasi-static problem, as given by theorem 3.8, we perform a standard finite difference method, since Ca ∈ R. We have that For this problem, we take Ca = 0.05, λ = 1 and the far-field (5.6a,b). The interface equation, (2.8), is evolved to T = 1 with Δτ = 10 −3 with the initial condition X (0) being a sphere of radius R = 1. The desired geometry X d is obtained by evolving the interface with the same set-up for Ca d = 0.1. We have not applied additional smoothing to this case, as we have found it skews the finite difference comparison. Choosing a small capillary number and a small time step instead ensures stability for the duration of the simulation.
In figure 4a, we present the convergence for ∈ {10 −d : d ∈ [[1, 7]]}. We can see that the error for a fixed discretization M = 32 behaves as expected, with an initial first-order accuracy. However, in the case from figure 4(b), we can see that the discretization error behaves considerably better than expected achieving close to fourth-order convergence. This is an artefact of the fact that λ = 1, resulting in smaller errors as in table 2, and that the deformation in the drop is small for Ca = 0.05. We present in table 3, the values for the case where Ca = 0.1 and Ca d = 0.05 and obtain values closer to the expected third order.

Optimization
Finally, we will look at finding close to optimal solutions to the optimization problems (P1) and (P2). These test cases use all the machinery described in § 4. In particular, we  will allow post processing of the interface during its evolution and during the gradient descent (as described in § 4.8).

Test 3: shape optimization
The first test we will perform is a classic shape optimization problem based on (P1). For this problem, we use a high viscosity ratio λ = 48.02 and capillary number of Ca = 1, corresponding roughly to a water droplet suspended in air. The far-field boundary conditions for the velocity and pressure are given by The initial shape in the optimization is a sphere of radius R = 1 and the desired normal velocity field u d in (3.1) is given by the solution on a geometry from Veerapaneni et al.
, namely r(ξ ) (a cos 2 (2π ξ)) 2 + (b sin(2π ξ)) 2 ) + cos 2 (2πk ξ), with a = 1, b = 1 and k = 1. The initial geometries and velocity fields are shown in figure 5. The optimization is performed using a standard conjugate gradient descent, as described in algorithm 1 from § 4.8, using the Fletcher-Reeves update formula from Hager & Zhang (2006). We choose a regularization parameter of 10 −2 (Eilers 2003) and the smoothing is performed at every step of the optimization. Finally, the discretization continues to follow the parameters described in table 1 with M = 64 elements. The optimal solutions are shown in figure 5. We can see that the algorithm reaches both the desired geometry and the desired normal velocity field. In figure 6 we have the per-iteration values of the normalized cost functional and gradient, defined aŝ The optimization algorithm was terminated after 128 iterations. We can see that the cost functional itself was decreased by approximately six orders of magnitude and plateaus when approaching the optimum. We have found that, with the chosen smoothing, the optimization is sufficiently robust and optimal solutions can be retrieved, with reasonable choices of initial conditions. Of course, this is a highly nonlinear optimization problem, so we do expect the existence of local minima, which a gradient descent algorithm cannot overcome.

Test 4: matching droplet evolution
Next, we look at a quasi-static optimization problem. The physical set-up we will be considering is Taylor's 'four roller' apparatus, described and analysed in Taylor (1934) (see figure 7). The same deformation problem of a droplet in extensional flow was studied in Stone & Leal (1989). The problem is one of boundary control, where the boundary conditions at 'infinity' are given by the (dimensional) velocity field  with C > 0, the strain rate, as a control parameter. In the non-dimensional equations, this amounts to using the capillary number as a control variable. The viscosity ratio is λ = 48.02. This is a high viscosity ratio, similar to Taylor's droplet formed from a tar-pitch mixture. The optimization is set-up as follows. The initial guess is given by Ca (0) = 0.1. The desired geometry in (3.1) is obtained by evolving (2.8) with a capillary number Ca d = 0.05. This allows achieving J 2 ≈ 0. For the optimization, we use the steepest descent choice of the descent direction, since the control parameter is a scalar. The optimization is run until convergence.
Finally, to evolve the interface, we use the third-order strong-stability-preserving Runge-Kutta method described in § 4.5, with an initial condition X (0), a sphere of radius R = 1. The simulation is run to a final time T = 1 and a time step of Δτ = 10 −3 , with a coarse mesh of M = 16 elements and the parameters from table 1. We also include the tangential forcing described in Loewenberg & Hinch (1996) and a smoothing filter from Eilers (2003), with a filter strength of 10 −4 . The results of the optimization can be seen in figure 8. As before, we report the cost and gradient normalized by their initial values. We can see that the cost functional converges to machine precision in only 14 conjugate gradient iterations.

Test 5: propulsion by Marangoni effects
Finally, we present a quasi-static problem where the control is a prescribed variable surface tension. Such a variation in surface tension can be achieved by a varying surfactant concentration, temperature gradients or electric fields. The first two cases have been studied in Muradoglu & Tryggvason (2008), for example, to which we refer for further details. We will attempt to find an optimal surfactant distribution on the droplet surface to maximize its propulsion in an otherwise quiescent flow. Due to the variable surface tension, the jump in traction becomes (Prosperetti & Tryggvason 2009) where we take a non-dimensional surface tension that is a perturbation of the constant case γ 1 + γ, (5.15) for a small . The control is taken to beγ in the cost functional (5.16) where x d is a desired centroid and x c (t) is the centroid at time t, as given by 1 2 x 2 n dS. (5.17) As this is a different optimization problem, we must expand on the results from § 3 to allow for the proposed changes. First, we have that the shape gradient and adjoint-based gradient of the cost functional are To finalize the definitions from lemma 3.7, the source term arising from the jump conditions becomes where j γ κn + ∇ Σ γ . An important difference in this case is that the last term now involves an averaging of the normal adjoint velocity gradient, since j also has a tangential component. All the terms can be obtained by following the steps described in appendix C.
For the fluid problem, we choose a quiescent flow u ∞ = 0 and p ∞ = 0 with an initial condition X (0) given as a sphere of radius R = 1. The parameters are given by Ca 0 = 0.05, λ = 1 and = −0.15. The equations are evolved to a final time of T = 1 with a Courant number of θ = 0.75 in (4.15). We have found that the resulting time steps ensure stability for all the configurations found during the optimization.
The optimization is started with an initial guess of (5.20) which is symmetric around θ = π/2, as seen in figure 9(b). Due to this symmetry, the initial guess does not transport the droplet, but simply changes its shape, while remaining centred at x c (0) = 0. The desired centroid is taken to be at a distance of 10 times the radius of the initial sphere, i.e. we have x d = (10R, 0). For the descent direction, we have used the simpler steepest descent. The cost functional and the controlγ during the optimization can be seen in figure 9 and figure 10. As before we can see that the cost correctly converges to the desired centroid location with a final value that approaches machine precision. The optimization was stopped when a step size that satisfied the Armijo condition could not be found, due to the small value of the cost functional. Additional results for a non-symmetric initial distribution ofγ have shown a clear tendency to evolve the droplets in the direction of larger surface tension, as also shown in Muradoglu & Tryggvason (2008) for a more general case. This is a well known result for variable surface tension problems.

Conclusions
In the current work, we have presented a derivation of the adjoint equations and adjoint-based gradient for two multiphase flow problems with a sharp interface. On a theoretical level, the main difficulty in handling sharp interfaces lies in the fact that the problem variables are not shape differentiable over the whole domain and must be considered separately over each fluid domain. We extend preliminary results from scalar problems with discontinuous coefficients (Pantz 2005) to the vector Stokes problem. This allows us to derive a general set of adjoint equations for static and quasi-static cases.
The optimization problems we have investigated naturally reduced to the fluid-fluid interface. However, the resulting gradients are highly non-trivial functions of the state and adjoint variables, as well as their derivatives (see § 3). This makes a numerical approximation of the shape gradients difficult. To work around these difficulties, we use boundary integral methods, which allow high-order representations of boundary restrictions for all the quantities of interest. As we have reported in § 5, using boundary integral methods we are able to numerically verify the adjoint-based gradients to at least single precision. We present two optimization problems -a static shape optimization problem and a quasi-static interface tracking optimization -that also reach their known minima to orders of 10 −8 or less. A non-trivial application to variable surface tension problems is also given, showing how the current results can be extended to more complex scenarios.
The numerical results we have presented are axisymmetric and reduce to one-dimensional interfaces. However, the theoretical results concerning the adjoint equations and adjoint-based gradients are independent of the spatial dimension. A fully three-dimensional implementation would allow for a much wider variety of applications. We believe the current work sufficiently motivates the continuous equations and prepares for such future applications. In particular, we anticipate the formulation and results to be useful in verifying forthcoming finite-inertia multiphase flow solvers with adjoint capabilities and in extensions to droplet control with electrostatic or other body forces.
Even so, many open questions remain. A more rigorous footing is required regarding the existence and uniqueness of optimal controls for the presented optimization problems. In general, we have not shown that the adjoint equations and adjoint-based gradients have the required regularity. Extensions to different classes of cost functionals are also required in practice, e.g. problems of drag minimization are very common. However, their shape differentiability in the multiphase case is rarely obvious. Finally, the numerical methods required to accurately represent shape gradients, such as the ones presented in theorem 3.5, are not trivial to implement in the general three-dimensional case. A possible route forward in this regard would require deriving an equivalent volumetric representation of the shape gradient (Giacomini, Pantz & Trabelsi 2017;Feppon et al. 2019), which requires less regularity.
Funding. This work was sponsored, in part, by the Office of Naval Research (ONR) as part of the Multidisciplinary University Research Initiatives (MURI) Program, under grant number N00014-16-1-2617.
Declaration of interests. The authors report no conflict of interest.

Appendix A. Notation
We give here a brief description of the notation and main auxiliary results used throughout the paper. We start by defining some functional spaces that appear.
We denote by C k (Ω), with Ω ⊂ R d , the space of all k-times continuously differentiable functions equipped with the uniform convergence norm. We also define W s,p (Ω) as the standard Sobolev spaces of real order s with derivatives in L p (Ω), as described in Temam (2001, chapter 1), Delfour & Zolésio (2001) and the references therein. We denote by H s (Ω) W s,2 (Ω) the corresponding family of Hilbert spaces. If a mapping f is a function of both time and space, we make use of the common notation for Bochner spaces, i.e. we say that f ∈ X ([0, T]; Y) if f (·, x) ∈ X([0, T]), for fixed x, and f (t, ·) ∈ Y, for fixed t. A similar notation will be employed for vector quantities, where we say that v ∈ X(Ω; Let X be a Banach space, then we denote its topological dual by X * . If A : X → Y is a linear operator, where X and Y are Banach spaces, then we define its adjoint by the linear operator A * : Y * → X * that satisfies for all x ∈ X and y * ∈ Y * , where ·, · X×X * denotes the dual mapping. In the case of Hilbert spaces, such as H k (Ω) W 2,k (Ω), we can identify the dual with the space itself and the dual mapping with the chosen inner product. We have the following definition of the Fréchet derivative from Luenberger (1997

then f is said to be Fréchet differentiable at x and Df (x)[h] is its Fréchet derivative in the direction h.
If the space X is the direct sum of spaces X n , i.e. X = n X n , then we define the following notation for the partial derivative ∂f ∂x n [h n ] D(x) [(0, . . . , h n , . . . , 0)].
If the space Y = R, f is usually called a functional. In this case, we can also define the gradient of f , if it exists, as the point ∇f (x) ∈ X * that satisfies for the chosen dual mapping on X. If X is a Hilbert space, the existence of the gradient is guaranteed by the Riesz representation theorem (Luenberger 1997, § 5.3). If X R d , we employ the usual notions of tensor calculus on Cartesian spaces. To avoid ambiguities, we will define our tensor derivatives using the common summation convention over repeated indices. We have that, for a vector u and a tensor σ , are the vector gradient and tensor divergence, respectively. We also denote the left and right multiplications between a vector and a tensor as (n · σ ) j n i σ ij and (σ n) i σ ij n j . (A6a,b) Finally, we use the single dot notation for the Frobenius inner product of two second-order tensors, i.e. A · B A ij B ij . Of special interest in our case are notions of tangential calculus (Delfour & Zolésio 2001). The divergence theorem and Reynolds transport theorem are well known on a volume Ω and the following theorems provide extensions to moving surfaces. THEOREM A.2 (Surface divergence theorem). Let v ∈ C 1 (Σ; R d ), where Σ is a smooth manifold of dimension d − 1. Then where κ denotes the additive curvature (twice the mean curvature) and the surface divergence is defined in Delfour & Zolésio (2001, § 5.4). The conormal is defined as μ τ × n, where τ is the boundary tangent positively oriented with respect to n.  (e x · f ± )ρ ds.
We compare the approximation of the above integral with the analytic results C D,+ = 2π 2 + 3λ 1 + λ , We choose the same discretization parameters as before, with a viscosity ratio of λ = 5 and capillary number of Ca = 0.01. The resulting relative errors can be observed in table 5.

E.3. Test 3: droplet deformation
Finally, we perform a test for the evolution of the droplet described in (2.8). This test is based on the work of Stone & Leal (1989), where we will recreate figure 4 in their paper.  The test case tracks the deformation of the droplet by looking at where L and B are the distances from the right-most and top-most points on the droplet to the centre of mass (here the origin). Steady state is determined based on the magnitude of the normal velocity, i.e. the simulation is stopped if for = 10 −4 . The far-field conditions are given by the extensional flow (5.6a,b). The methods used in our implementation differ significantly from the ones presented in Stone & Leal (1989). We use a third-order Runge-Kutta time integrator and the tangential corrections proposed in Loewenberg & Hinch (1996) to cluster points towards in regions of higher curvature. The time step is computed at each iteration according to (4.15) with a Courant number of θ = 0.1. We use a coarse grid with only M = 16 elements and the same number of collocation and quadrature points per element as reported in table 1.
The droplet deformation is shown in figure 12. We can see that our results match those obtained by Stone & Leal (1989) and also agree very well with the asymptotic expansions for small Ca. Finally, in figure 13, we report the relative volume loss In Stone & Leal (1989), the authors report a loss of 1 % to 10 % of the volume of the droplet over the span of 600-3000 time steps. In our implementation, we can see a relative volume loss of 10 −5 after 200-3000 iterations.