**Shenan Grossberg 1, Daniel S. Jarman <sup>2</sup> and Gavin R. Tabor 1,\***


Received: 24 January 2020; Accepted: 07 March 2020; Published: 11 March 2020

**Abstract:** The continuous adjoint approach is a technique for calculating the sensitivity of a flow to changes in input parameters, most commonly changes of geometry. Here we present for the first time the mathematical derivation of the adjoint system for multiphase flow modeled by the commonly used drift flux equations, together with the adjoint boundary conditions necessary to solve a generic multiphase flow problem. The objective function is defined for such a system, and specific examples derived for commonly used settling velocity formulations such as the Takacs and Dahl models. We also discuss the use of these equations for a complete optimisation process.

**Keywords:** adjoint optimization; multiphase flow; computational fluid dynamics

#### **1. Introduction**

The adjoint method is currently attracting significant interest as an optimization process in CFD. The objective of the adjoint approach is to calculate the sensitivity of the flow solution with respect to changes in the input parameters, most commonly changes in the geometry. This can then in principle be used as the basis for an iterative optimization algorithm based on gradient information (the sensitivities) which can optimize the design with many fewer function evaluations than would be the case for non-gradient-based approaches (such as genetic algorithms). Calculating the sensitivities requires differentiating the governing equations with respect to the changes of the input parameters, and since the governing equations for fluid flow are the Navier-Stokes equations (or equations derived from these such as the Reynolds Averaged Navier Stokes equations), this is understandably very challenging. There are two main approaches; the discrete adjoint approach, and the continuous adjoint approach. In the discrete adjoint approach, the sensitivity matrix is calculated numerically by evaluating the system for small changes in the inputs and applying standard finite difference methods. In the continuous adjoint approach, the sensitivities are calculated mathematically using lagrange multipliers. This is more elegant and provides an implementation which is easier to code, requires fewer evaluations and can be made numerically consistent with the evaluation of the original equation set. However it does require significant mathematical analysis in advance, and if the problem formulation changes (different equations, boundary conditions etc) this has to be repeated. Examples of the application of the continuous adjoint method for single phase flow can be found in a range of areas [1,2] such as automotive [3–5], aerospace [6,7] and turbomachinery [8–10], and implementations of the equations can be found in general purpose CFD codes such as STAR-CCM, ANSYS Fluent [11] and Engys Helyx [4]. However the equations are complex to develop and application to multiphase systems is only just starting [12]. In many cases, even just the evaluation of the sensitivities is valuable, as they can be used to indicate possible changes to the design engineers. Beyond this the sensitivities can also be used as the basis for an optimization loop [2]. This of course necessitates the morphing of

the geometry through techniques such as volumetric B-splines [13] or Radial Basis Functions [14] and consequent updating of the mesh [15].

Multiphase flow is the simultaneous flow of two or more immiscible phases in a system. In dispersed multiphase systems, one or more of the phases exists as fluid particles small enough not to be resolved in the simulation; examples include gas bubbles in water, emulsions (liquid droplets in another immiscible liquid) and actual solid particles in gas or liquid. A wide variety of different mathematical models have been derived over the years to describe dispersed multiphase flow, including mixture models, lagrangian particle tracking, and eulerian n-fluid models [16,17]. Which is used depends on the exact physics of the problem, as well as factors such as available computing resources and desired accuracy. In many physical systems, the density ratio between the two phases is low, generally less than 2:1, and the drag force between them is high. Therefore, to a good approximation, the two phases can be considered to respond to pressure gradients as a single phase. Additionally, the slip (drift) between the phases is primarily due to the gravitational settling of the dispersed phase. This might adequately describe solid particles in water or an emulsion of immiscible liquids, and in these cases a commonly used mathematical model is the drift flux model. Hence it is this set of equations we have decided to focus on.

In the drift flux model, the two phases are treated as one: the momentum and continuity equations for both phases are summed to create a mixture-momentum and mixture-continuity equation, and the transport of the dispersed phase is modelled using a drift equation. The three equations, collectively called the drift flux equations, are listed below:

$$\begin{split} \frac{\partial \rho\_m \mathbf{v}\_m}{\partial t} + (\mathbf{v}\_m \cdot \nabla)(\rho\_m \mathbf{v}\_m) &= -\nabla(\rho\_m p\_m) + \nabla \cdot \left(2\mu\_m \mathbf{D}(\mathbf{v}\_m)\right) \\ &- \nabla \cdot \left(\frac{\kappa}{1-\alpha} \frac{\rho\_c \rho\_d}{\rho\_m} \mathbf{v}\_{dj} \mathbf{v}\_{dj}\right) + \rho\_m \mathbf{g} + F\_\prime \end{split} \tag{1a}$$

$$\frac{\partial \rho\_m}{\partial t} + \nabla \cdot (\rho\_m \mathbf{v}\_m) = 0,\tag{1b}$$

$$\frac{\partial a}{\partial t} + \nabla \cdot (a\mathbf{v}\_m) = -\nabla \cdot \left(\frac{a\rho\_\mathbf{c}}{\rho\_m}\mathbf{v}\_{dj}\right) + \nabla \cdot K \nabla a\_\mathbf{v} \tag{1c}$$

where:


In summing the momentum equations, not only have the number of equations been reduced from four to three, but the inter-phase momentum transfer terms have also been eliminated which were numerically unstable [18]. Hence, a far more robust equation set has been produced and the computational resources required to solve the system have been reduced. This also makes it a very appropriate basis from which to develop an adjoint formulation suitable for applying to dispersed

multiphase flows in this regime. This is the challenge of the current paper. We focus in particular on wall-bounded or ducted flows, in which there is no contribution to the objective function from the interior of the domain, in other words, the performance of the system is entirely governed by the boundary properties.

The paper is organized as follows. The optimization problem is stated in Section 2 and the adjoint equations for the drift flux model are derived for the general case in Section 2. These equations are then applied to the specific case of ducted or wall-bounded flows in Section 3, with the objective function for this case being specified in Section 4, and different settling velocities in Section 5. Finally, the conclusions follow in Section 6.

#### **2. The Optimization Problem**

If the performance of a device is measured by an objective function, *J*, and the residuals of the primal (flow) equations are given by R, the optimisation problem can be stated as,

$$\text{optimise } f(\mathbf{x}, \mathbf{y}) \text{ subject to } \mathcal{R}(\mathbf{x}, \mathbf{y}) = \mathbf{0},\tag{2}$$

where x are the design parameters and y are the primal variables [19]. It can then be formulated as,

$$
\mathcal{L}\mathcal{L} = \mathcal{J} + \int\_{\Omega} \lambda \,\mathcal{R} \,d\Omega,\tag{3}
$$

where L is the Lagrange function, *λ* are the Lagrange multipliers (also referred to as the adjoint variables) and Ω is the flow domain. In this case, the primal equations are the steady state drift flux equations, with the capillary force taken to be zero [18] and a Darcy term included in the mixture-momentum equation. They are rearranged in terms of their residuals, R = (*R*1, *<sup>R</sup>*2, *<sup>R</sup>*3, *<sup>R</sup>*4, *<sup>R</sup>*5)*T*, as follows:

$$\begin{aligned} (R\_1, R\_2, R\_3)^T &= (\mathbf{v}\_m \cdot \nabla)(\rho\_m \mathbf{v}\_m) + \nabla(\rho\_m p\_m) - \nabla \cdot \left(2\mu\_m \mathbf{D}(\mathbf{v}\_m)\right) \\ &+ \nabla \cdot \left(\frac{a}{1-a} \frac{\rho\_c \rho\_d}{\rho\_m} \mathbf{v}\_{dj} \mathbf{v}\_{dj}\right) - \rho\_m \mathbf{g} + \mathbb{N}\rho\_m \mathbf{v}\_{m\prime} \end{aligned} \tag{4a}$$

$$\mathcal{R}\_4 = -\nabla \cdot (\rho\_m \mathbf{v}\_m)\_\prime \tag{4b}$$

$$R\_5 = \nabla \cdot (a\mathbf{v}\_{\rm int}) + \nabla \cdot \left(\frac{a\rho\_c}{\rho\_{\rm m}}\mathbf{v}\_{d\dot{\jmath}}\right) - \nabla \cdot K\nabla a\_\prime \tag{4c}$$

where ℵ is the porosity, associated with the Darcy term. The variation of the Lagrange function with respect to the primal variables, (**v***m*, *pm*, *α*), and the design parameter, ℵ, is,

$$
\delta \lrcorner \mathcal{C} = \delta\_{\mathbf{v}\_m} \mathcal{C} + \delta\_{p\_m} \mathcal{C} + \delta\_{\mathbf{a}} \mathcal{C} + \delta\_{\mathbb{N}} \mathcal{C} \,,\tag{5}
$$

where, for example, *δα* L = L (*α* + *δα*) − L (*α*). We choose the adjoint variables, (**u**, *q*, *β*) = (*u*1, *u*2, *u*3, *q*, *β*), so that the variation with respect to the primal variables vanishes, i.e.,

$$
\delta\_{\mathbf{v}\_{m}} \mathcal{X}^{\rho} + \delta\_{p\_{m}} \mathcal{X}^{\rho} + \delta\_{n} \mathcal{X}^{\rho} = 0,\tag{6}
$$

and the Lagrange function now varies only with respect to the design parameter,

$$
\delta \mathcal{L}^{\rho} = \delta\_{\mathbb{N}} \mathcal{L}^{\rho} = \delta\_{\mathbb{N}} \mathcal{I} + \int\_{\Omega} (\mathbf{u}, q, \boldsymbol{\rho}, \boldsymbol{\beta}) \delta\_{\mathbb{N}} \, \mathbb{R} \, \, d\Omega. \tag{7}
$$

#### *Derivation of the Adjoint Drift Flux Equations*

The adjoint drift flux equations are derived by substituting Equation (3) into Equation (6), giving,

$$\begin{aligned} \delta\_{\mathbf{v}\_{n}} \| + \delta\_{p\_{n}} \| + \delta\_{\mathbf{n}} \| \\ &+ \int\_{\Omega} (\mathbf{u}, q, \boldsymbol{\rho}, \boldsymbol{\xi}) \delta\_{\mathbf{v}\_{n}} \mathcal{R} \, d\Omega + \int\_{\Omega} (\mathbf{u}, q, \boldsymbol{\rho}, \boldsymbol{\xi}) \delta\_{p\_{n}} \mathcal{R} \, d\Omega + \int\_{\Omega} (\mathbf{u}, q, \boldsymbol{\rho}, \boldsymbol{\xi}) \delta\_{\mathbf{n}} \mathcal{R} \, d\Omega = 0, \end{aligned} \tag{8}$$

which can be expanded to,

$$\begin{aligned} \delta \boldsymbol{\delta}\_{\mathbf{w}} \boldsymbol{f} + \delta \boldsymbol{\delta}\_{\mathbf{p}} \boldsymbol{f} + \delta \boldsymbol{\delta}\_{\mathbf{f}} &+ \int\_{\Omega} d\Omega \, \mathbf{u} \cdot \delta \mathbf{v}\_{\mathbf{n}} (\mathcal{R}\_{1}, \mathcal{R}\_{2}, \mathcal{R}\_{3})^{T} + \int\_{\Omega} d\Omega \, q \delta \boldsymbol{\delta}\_{\mathbf{u}} \mathcal{R}\_{4} + \int\_{\Omega} d\Omega \, \boldsymbol{\delta} \delta \boldsymbol{\delta}\_{\mathbf{v}} \mathcal{R}\_{5} \\ &+ \int\_{\Omega} d\Omega \, \mathbf{u} \cdot \delta \boldsymbol{\delta}\_{\mathbf{f}} (\mathcal{R}\_{1}, \mathcal{R}\_{2}, \mathcal{R}\_{3})^{T} + \int\_{\Omega} d\Omega \, q \delta\_{p\_{\mathbf{n}}} \mathcal{R}\_{4} + \int\_{\Omega} d\Omega \, \boldsymbol{\delta} \delta \boldsymbol{\delta}\_{p\_{\mathbf{n}}} \mathcal{R}\_{5} \\ &+ \int\_{\Omega} d\Omega \, \mathbf{u} \cdot \delta\_{\mathbf{a}} (\mathcal{R}\_{1}, \mathcal{R}\_{2}, \mathcal{R}\_{3})^{T} + \int\_{\Omega} d\Omega \, q \delta\_{\mathbf{a}} \mathcal{R}\_{4} + \int\_{\Omega} d\Omega \, \boldsymbol{\delta} \delta \boldsymbol{\delta}\_{\mathbf{a}} \mathcal{R}\_{5} = 0. \end{aligned}$$

The variation of R with respect to the primal variables can be determined as:

$$\begin{split} \delta\_{\mathbf{v}\_{m}} (R\_{1}, R\_{2}, R\_{3})^{T} &= (\delta \mathbf{v}\_{m} \cdot \nabla) (\rho\_{m} \mathbf{v}\_{m}) + (\mathbf{v}\_{m} \cdot \nabla)(\rho\_{m} \delta \mathbf{v}\_{m}) - \nabla \cdot \left( 2 \mu\_{m} \mathbf{D} (\delta \mathbf{v}\_{m}) \right) \\ &- \nabla \cdot \left( 2 \delta\_{\mathbf{v}\_{m}} \mu\_{d} \mathbf{D} (\mathbf{v}\_{m}) \right) + \aleph \rho\_{m} \delta \mathbf{v}\_{m} \end{split} \tag{10a}$$

$$
\delta\_{\mathbf{v}\_m} \mathcal{R}\_{\mathbf{4}} = -\nabla \cdot (\rho\_m \delta \mathbf{v}\_m),
\tag{10b}
$$

$$
\delta\_{\mathbf{v}\_m} R\_5 = \nabla \cdot (a \delta \mathbf{v}\_m),
\tag{10c}
$$

$$\delta\_{p\_m}(R\_1, R\_2, R\_3)^T = \nabla(\rho\_m \delta p\_m),\tag{10d}$$

$$
\delta\_{p\_m} R\_4 = 0,
$$

$$
\delta\_{p\_m} R\_{\overline{5}} = 0,
$$

$$\begin{split} \delta\_{\mathsf{a}} (R\_1, R\_2, R\_3)^T &= (\rho\_d - \rho\_c) \left( (\mathbf{v}\_m \cdot \nabla)(\delta \mathbf{a} \mathbf{v}\_m) + \nabla(\delta \mathbf{a} p\_m) + \delta \mathbf{a} (\mathbb{N} \mathbf{v}\_m - \mathbf{g}) \right) \\ &- \nabla \cdot \left( 2 \delta\_{\mathsf{a}} \mu\_d \mathbf{D} (\mathbf{v}\_m) \right) + \nabla \cdot \delta\_{\mathsf{a}} (a \rho\_d \mathbf{v}\_d \mathbf{v}\_{di}), \end{split} \tag{10g}$$

$$-\nabla \cdot \left(2\delta\_{\mathfrak{a}}\mu\_{d}\mathbf{D}(\mathbf{v}\_{\mathfrak{m}})\right) + \nabla \cdot \delta\_{\mathfrak{a}}(\mathsf{a}\rho\_{d}\mathbf{v}\_{d\bar{j}}\mathbf{v}\_{d\bar{j}})\_{\mathfrak{l}} \tag{10g}$$

$$
\delta\_{\mathsf{A}} R\_{\mathsf{A}} = - (\rho\_d - \rho\_{\mathsf{c}}) \nabla \cdot (\delta a \mathbf{v}\_{\mathsf{m}}),
\tag{10h}
$$

$$\begin{split} \delta\_{a}R\_{5} &= \nabla \cdot (\delta a \mathbf{v}\_{m}) + \nabla \cdot \delta\_{a} (a \mathbf{v}\_{dj}) - \nabla \cdot \left(\frac{\mu\_{m}^{t}}{\rho\_{c}} \nabla \delta a\right) \\ &+ \frac{\mu\_{m}^{t}}{\rho\_{c}} \left(\frac{\rho\_{d}}{\rho\_{c}} - 1\right) \nabla \cdot (\delta a \nabla a) + \frac{\mu\_{m}^{t}}{\rho\_{c}} \left(\frac{\rho\_{d}}{\rho\_{c}} - 1\right) \nabla \cdot (a \nabla \delta a). \end{split} \tag{10i}$$

Derivation of Equations (10a), (10g) and (10i) can be found in Appendies A–C, respectively, where the variation of *μ<sup>t</sup> <sup>m</sup>* has been neglected. This is correct only for laminar flow regimes. For turbulent flows, neglecting this variation constitutes a common approximation, known as *frozen turbulence* [19]. This may introduce errors into the optimisation [20], although there are cases in the literature where the frozen turbulence assumption can be demonstrated to be acceptable [21].

With these variations, Equation (9) now reads,

*δ***v***<sup>m</sup> J* + *δpm J* + *δα J* + Ω *d*Ω **u** · (*δ***v***<sup>m</sup>* · ∇)(*ρm***v***m*)+(**v***<sup>m</sup>* · ∇)(*ρmδ***v***m*) −∇·  2*μm*D(*δ***v***m*) −∇·  2*δ***v***<sup>m</sup> μd*D(**v***m*) + ℵ*ρmδ***v***<sup>m</sup>* − Ω *<sup>d</sup>*<sup>Ω</sup> *<sup>q</sup>*∇ · (*ρmδ***v***m*) + Ω *d*Ω *β*∇ · (*αδ***v***m*) + Ω *d*Ω **u** · ∇(*ρmδpm*)+(*ρ<sup>d</sup>* − *ρc*) Ω *d*Ω **u** · (**v***<sup>m</sup>* · ∇)(*δα***v***m*) + ∇(*δαpm*) + *δα*(ℵ**v***<sup>m</sup>* − **g**) − Ω *d*Ω **u** · ∇ ·  2*δαμd*D(**v***m*) + Ω *d*Ω **u** · ∇ · *δα*(*αρd***v***dj***v***dj*) − (*ρ<sup>d</sup>* − *ρc*) Ω *<sup>d</sup>*<sup>Ω</sup> *<sup>q</sup>*∇ · (*δα***v***m*) + Ω *d*Ω *β* ∇ · (*δα***v***m*) + ∇ · *δα*(*α***v***dj*) −∇· *μ<sup>t</sup> m ρc* ∇*δα* + *ρ<sup>d</sup> ρc* − 1 ∇ · *μ<sup>t</sup> m ρc δα*∇*α* + *ρ<sup>d</sup> ρc* − 1 ∇ · *μ<sup>t</sup> m ρc <sup>α</sup>*∇*δα* <sup>=</sup> 0. (11)

Decomposing the objective function into contributions from the boundary, Γ, and interior, Ω, of the domain,

$$J = \int\_{\Gamma} J\_{\Gamma} \, d\Gamma + \int\_{\Omega} J\_{\Omega} \, d\Omega\_{\prime} \tag{12}$$

Equation (11) can be reformulated as,

$$\begin{split} \int\_{\Gamma} \mathrm{d}\Gamma \left( \mathbf{n} (\mathbf{u} \cdot \rho\_{\mathbf{w}} \mathbf{w}\_{n}) + \mathbf{u} (\rho\_{\mathbf{w}} \mathbf{w}\_{n} \cdot \mathbf{n}) + 2 \mu\_{\mathbf{m}} \mathbf{n} \cdot \mathbf{D} (\mathbf{u} - \mathbf{q} \cdot \mathbf{n} + \mathbf{a} \beta \mathbf{n} + \beta \mathbf{n} \frac{\delta \mathbf{v}}{\delta \mathbf{v}\_{n}}) \cdot \delta \mathbf{v}\_{m} \right. \\ \quad - \int\_{\Gamma} \mathrm{d}\Gamma \left( 2 \mu\_{\mathbf{m}} \mathbf{n} \cdot \mathbf{D} (\delta \mathbf{v}\_{m}) \cdot \mathbf{n} + 2 \delta\_{\mathbf{w}} \mu\_{\mathbf{m}} \mathbf{n} \cdot \mathbf{D} (\mathbf{w}\_{n}) \cdot \mathbf{u} - 2 \delta\_{\mathbf{v} \mathbf{n}} \mu\_{\mathbf{m}} \mathbf{u} \cdot \mathbf{D} (\mathbf{u} \cdot \mathbf{v}\_{m}) \right. \\ \left. + \int\_{\Gamma} \mathrm{d}\Omega \left( -\nabla \mathbf{u} \cdot (\rho\_{\mathbf{w}} \mathbf{v}\_{n}) - (\rho\_{\mathbf{m}} \mathbf{v}\_{n} \cdot \nabla) \mathbf{u} - \nabla \cdot (2 \mu\_{\mathbf{m}} \mathbf{D} (\mathbf{u})) \right. \\ \left. \quad + \mathrm{N} \rho\_{\mathbf{m}} \mathbf{u} + \mu\_{\mathbf{m}} \nabla q - a \, \nabla \delta \frac{\partial \mathbf{v}\_{n}}{\partial \mu\_{\mathbf{m}}} \right. \\ \left. + \int\_{\Gamma} \mathrm{d}\Gamma \left( \rho\_{\mathbf{m}} \mathbf{u} \cdot \mathbf{n} + \frac{\delta \mathbf{J} \nu$$

*Fluids* **2020**, *5*, 31

Derivation of Equation (13) can be found in Appendix D. In order to satisfy Equation (13) in general, the integrals must vanish individually. The adjoint drift flux equations are deduced from the integrals over the interior of the domain:

$$\begin{split}-\nabla\mathbf{u}\cdot\left(\rho\_{\text{m}}\mathbf{v}\_{\text{m}}\right)-\left(\rho\_{\text{m}}\mathbf{v}\_{\text{m}}\cdot\nabla\right)\mathbf{u}-\nabla\cdot\left(2\mu\_{\text{m}}\mathbf{D}(\mathbf{u})\right) &= -\rho\_{\text{m}}\nabla q + a\nabla\beta - \aleph\_{\mathcal{P}m}\mathbf{u} - \frac{\partial I\_{\Omega}}{\partial \mathbf{v}\_{\text{m}}} \\ &+ \nabla\cdot\left(2\frac{\partial \mu\_{d}}{\partial \mathbf{v}\_{\text{m}}}\mathbf{D}(\mathbf{u})\right)\cdot\mathbf{v}\_{\text{m}} \\ &- \dots \end{split} \tag{14a}$$

$$\nabla \cdot (\rho\_m \mathbf{u}) = \frac{\partial f\_{\Omega}}{\partial p\_m},\tag{14b}$$

$$-\left(\mathbf{v}\_m + \frac{\partial}{\partial a}(a\mathbf{v}\_{dj}) + \frac{\mu\_m^t}{\rho\_c}\left(\frac{\rho\_d}{\rho\_c} - 1\right)\nabla a\right) \cdot \nabla \beta = \nabla \cdot \frac{\mu\_m^t}{\rho\_c}\left(1 - a\left(\frac{\rho\_d}{\rho\_c} - 1\right)\right)\nabla \beta$$

$$+ S\_1 + S\_2 - \frac{\partial f\_{\Omega}}{\partial a},\tag{14c}$$

where:

$$S\_1 = (\rho\_d - \rho\_{\mathbb{C}}) \left( (\mathbf{v}\_{\mathcal{W}} \cdot \nabla)(\mathbf{u} \cdot \mathbf{v}\_{\mathcal{W}} - q) - \mathbf{u} \cdot (\mathbb{N} \mathbf{v}\_{\mathcal{W}} - \mathbf{g}) \right) + \nabla \cdot \left( 2 \frac{\partial \mu\_d}{\partial \mathbf{a}} \mathbf{D}(\mathbf{u}) \right) \cdot \mathbf{v}\_{\mathcal{W}} \tag{15a}$$

$$S\_2 = \left( (\rho\_d - \rho\_c) p\_m + \frac{\partial}{\partial \mathbf{a}} (\mathbf{a} \rho\_d \mathbf{v}\_{d\uparrow} \mathbf{v}\_{d\uparrow}) \right) \nabla \cdot \mathbf{u}\_{\prime} \tag{15b}$$

and the boundary conditions for the adjoint variables are deduced from the surface integrals:

$$\int\_{\Gamma} d\Gamma \left( \mathbf{n} (\mathbf{u} \cdot \rho\_m \mathbf{v}\_m) + \rho\_m \mathbf{v}\_m \cdot \mathbf{n} \mathbf{u} + 2\mu\_m \mathbf{n} \cdot \mathbf{D}(\mathbf{u}) - q\rho\_m \mathbf{n} + a\beta \mathbf{n} + \frac{\partial f\_{\Gamma}}{\partial \mathbf{v}\_m} \right. $$

$$+ 2 \frac{\partial \mu\_d}{\partial \mathbf{v}\_m} \mathbf{n} \cdot \left( \mathbf{D}(\mathbf{v}\_m) \cdot \mathbf{u} - \mathbf{D}(\mathbf{u}) \cdot \mathbf{v}\_m \right) \Big) \cdot \delta \mathbf{v}\_m - \int\_{\Gamma} d\Gamma 2\mu\_m \mathbf{n} \cdot \mathbf{D}(\delta \mathbf{v}\_m) \cdot \mathbf{u} = 0,\tag{16a}$$

$$\int\_{\Gamma} d\Gamma \left( \rho\_m \mu\_\mathbf{n} + \frac{\partial f\_\Gamma}{\partial p\_m} \right) \delta p\_\mathbf{n} = 0,\tag{16b}$$

$$\begin{split} \int\_{\Gamma} d\Gamma \left( \left( \mathbf{v}\_{m} \cdot \mathbf{n} + \frac{\partial}{\partial a} (a \mathbf{v}\_{dj}) \cdot \mathbf{n} + \frac{\mu\_{m}^{t}}{\rho\_{c}} \left( \frac{\rho\_{d}}{\rho\_{c}} - 1 \right) (\mathbf{n} \cdot \nabla) a \right) \beta \right. \\ \left. + \frac{\mu\_{m}^{t}}{\rho\_{c}} \left( 1 - a \left( \frac{\rho\_{d}}{\rho\_{c}} - 1 \right) \right) (\mathbf{n} \cdot \nabla) \beta + \mathbf{C}\_{1} + \mathbf{C}\_{2} + \frac{\partial f\_{\Gamma}}{\partial a} \right) \delta a \\ - \int\_{\Gamma} d\Gamma \beta \frac{\mu\_{m}^{t}}{\rho\_{c}} \left( 1 - a \left( \frac{\rho\_{d}}{\rho\_{c}} - 1 \right) \right) (\mathbf{n} \cdot \nabla) \delta a = 0, \end{split} \tag{16c}$$

where:

$$\mathbf{C}\_{1} = \left( (\rho\_{d} - \rho\_{c})(\mathbf{u} \cdot \mathbf{v}\_{m} - q)\mathbf{v}\_{m} + 2\frac{\partial \mu\_{d}}{\partial a} \left( \mathbf{D}(\mathbf{u}) \cdot \mathbf{v}\_{\text{ff}} \right) \right) \cdot \mathbf{n},\tag{17a}$$

$$\mathcal{C}\_2 = \left( \left( (\rho\_d - \rho\_c) p\_m + \frac{\partial}{\partial a} (a \rho\_d \mathbf{v}\_{dj} \mathbf{v}\_{dj}) \right) \mathbf{u} - 2 \frac{\partial \mu\_d}{\partial a} \left( \mathcal{D} (\mathbf{v}\_m) \cdot \mathbf{u} \right) \right) \cdot \mathbf{n} \tag{17b}$$

and *u*<sup>n</sup> = **u** · **n** is the normal component of the adjoint velocity. This is the general form of the adjoint equation system for the steady state drift flux equations with Darcy porosity term and frozen turbulence.

#### **3. Application to Wall Bounded Flows**

Thus far in the paper we have presented the optimisation problem in as generic a way as possible. To proceed further with the derivation we now need to derive expressions for the boundary conditions, objective function and slip velocity. We will examine these for the case of wall-bounded or ducted

flows, for which there is no contribution to the objective function from the interior of the domain. So, in the cases where the objective function only involves integrals over the surface of the flow domain rather than over its interior, the adjoint equations reduce to:

$$\begin{split}-\nabla\mathbf{u}\cdot(\rho\_{m}\mathbf{v}\_{m})-(\rho\_{m}\mathbf{v}\_{m}\cdot\nabla)\mathbf{u}-\nabla\cdot\left(2\mu\_{m}\mathbf{D}(\mathbf{u})\right)&=-\rho\_{m}\nabla q+a\nabla\beta-\mathbb{N}\rho\_{m}\mathbf{u} \\ &+\nabla\cdot\left(2\frac{\partial\mu\_{d}}{\partial\mathbf{v}\_{m}}\mathbf{D}(\mathbf{u})\right)\cdot\mathbf{v}\_{m},\end{split}\tag{18a}$$

$$\nabla \cdot (\rho\_m \mathbf{u}) = 0,\tag{18b}$$

$$-\left(\mathbf{v}\_{\mathrm{m}} + \frac{\partial}{\partial \boldsymbol{a}} (\boldsymbol{a} \mathbf{v}\_{d\boldsymbol{j}}) + \frac{\mu\_{\mathrm{m}}^{\mathrm{f}}}{\rho\_{\mathrm{c}}} \left(\frac{\rho\_{\mathrm{d}}}{\rho\_{\mathrm{c}}} - 1\right) \nabla \boldsymbol{a}\right) \cdot \nabla \boldsymbol{\beta} = \nabla \cdot \frac{\mu\_{\mathrm{m}}^{\mathrm{f}}}{\rho\_{\mathrm{c}}} \left(1 - \boldsymbol{a} \left(\frac{\rho\_{\mathrm{d}}}{\rho\_{\mathrm{c}}} - 1\right)\right) \nabla \boldsymbol{\beta} \tag{18c}$$
 
$$+ S\_{1}. \tag{18c}$$

These equations no longer depend on the objective function, so when switching from one optimisation objective to another, they remain unchanged and only the boundary conditions have to be adapted to the specific objective function. Note that as a result of Equation (18b), ∇ · **u** = 0 [22] and, therefore, *S*<sup>2</sup> = 0.

For the adjoint boundary conditions, the terms in Equation (16a) involving *μ<sup>m</sup>* can be rewritten as,

$$\int\_{\Gamma} d\Gamma \, 2\mu\_{\mathfrak{m}} \mathbf{n} \cdot \left( \mathrm{D}(\mathbf{u}) \cdot \delta \mathbf{v}\_{\mathfrak{m}} - \mathrm{D}(\delta \mathbf{v}\_{\mathfrak{m}}) \cdot \mathbf{u} \right) \\ \quad = \int\_{\Gamma} d\Gamma \, \mu\_{\mathfrak{m}} \left( (\mathbf{n} \cdot \nabla) \mathbf{u} \cdot \delta \mathbf{v}\_{\mathfrak{m}} - (\mathbf{n} \cdot \nabla) \delta \mathbf{v}\_{\mathfrak{m}} \cdot \mathbf{u} \right) \tag{19}$$

(Ref. [19]) and therefore the adjoint boundary conditions, Equation (16), reduce to:

$$\begin{split} \int\_{\Gamma} d\Gamma \Big( \mathbf{n} \left( \mathbf{u} \cdot \rho\_{m} \mathbf{v}\_{m} \right) + \rho\_{m} \mathbf{v}\_{m} \cdot \mathbf{n} \mathbf{u} + \mu\_{m} (\mathbf{n} \cdot \nabla) \mathbf{u} - q \rho\_{m} \mathbf{n} + a \beta \mathbf{n} + \frac{\partial f\_{\Gamma}}{\partial \mathbf{v}\_{m}} \\ + 2 \frac{\partial \mu\_{d}}{\partial \mathbf{v}\_{m}} \mathbf{n} \cdot \Big( \mathbf{D} (\mathbf{v}\_{m}) \cdot \mathbf{u} - \mathbf{D} (\mathbf{u}) \cdot \mathbf{v}\_{\text{m}} \Big) \Big) \cdot \delta \mathbf{v}\_{m} - \int\_{\Gamma} d\Gamma \mu\_{\text{fl}} (\mathbf{n} \cdot \nabla) \delta \mathbf{v}\_{m} \cdot \mathbf{u} = 0, \end{split} \tag{20a}$$

 Γ

$$d\Gamma \left( \rho\_m u\_\mathbf{n} + \frac{\partial f\_\Gamma}{\partial p\_m} \right) \delta p\_\mathbf{m} = 0,\tag{20b}$$

$$\begin{split} \int\_{\Gamma} d\Gamma \left( \left( \mathbf{v}\_{m} \cdot \mathbf{n} + \frac{\partial}{\partial a} (a\mathbf{v}\_{dj}) \cdot \mathbf{n} + \frac{\mu\_{m}^{t}}{\rho\_{c}} \left( \frac{\rho\_{d}}{\rho\_{c}} - 1 \right) (\mathbf{n} \cdot \nabla) a \right) \beta \right. \\ \left. + \frac{\mu\_{m}^{t}}{\rho\_{c}} \left( 1 - a \left( \frac{\rho\_{d}}{\rho\_{c}} - 1 \right) \right) (\mathbf{n} \cdot \nabla) \beta + \mathbf{C}\_{1} + \mathbf{C}\_{2} + \frac{\partial I\_{\Gamma}}{\partial a} \right) \delta a \\ - \int\_{\Gamma} d\Gamma \beta \frac{\mu\_{m}^{t}}{\rho\_{c}} \left( 1 - a \left( \frac{\rho\_{d}}{\rho\_{c}} - 1 \right) \right) (\mathbf{n} \cdot \nabla) \delta a = 0. \end{split} \tag{20c}$$

In order to determine the boundary conditions of the adjoint variables, the boundary conditions imposed on the primal variables are listed in Table 1. We will derive expressions for the three main boundary conditions.

**Table 1.** Primal boundary conditions.


#### *3.1. Adjoint Boundary Conditions at the Inlet*

At an inlet, the primal velocity and dispersed-phase volume fraction are usually fixed, so,

$$
\delta \mathbf{v}\_m = 0 \text{ and } \delta \mathbf{a} = 0. \tag{21}
$$

The first integrals in Equations (20a) and (20c) therefore go to zero and Equation (20) reduces to:

$$\int\_{\Gamma} d\Gamma \mu\_{\rm{m}} (\mathbf{n} \cdot \nabla) \delta \mathbf{v}\_{\rm{m}} \cdot \mathbf{u} = 0,\tag{22a}$$

$$\int\_{\Gamma} d\Gamma \left( \rho\_m u\_\mathbf{n} + \frac{\partial f\_\Gamma}{\partial p\_m} \right) \delta p\_\mathbf{m} = 0,\tag{22b}$$

$$\int\_{\Gamma} d\Gamma \beta \frac{\mu\_m^t}{\rho\_c} \left(1 - a\left(\frac{\rho\_d}{\rho\_c} - 1\right)\right) (\mathbf{n} \cdot \nabla) \delta \mathbf{a} = 0. \tag{22c}$$

When both fluids are incompressible, ∇ · **v***<sup>m</sup>* = 0 [22], and as *δ***v***m*<sup>t</sup> = 0 along the inlet, (**n** · ∇)*δ***v***<sup>m</sup>* = (**n** · ∇)*δ***v***m*<sup>t</sup> [19], where **v***m*<sup>t</sup> is the tangential component of the mixture velocity. Hence, Equation (22) reduces to:

$$\int\_{\Gamma} d\Gamma \mu\_{\text{tt}} (\mathbf{n} \cdot \nabla) \delta \mathbf{v}\_{\text{int}} \cdot \mathbf{u}\_{\text{t}} = 0,\tag{23a}$$

$$\int\_{\Gamma} d\Gamma \left( \rho\_m u\_\mathbf{n} + \frac{\partial f\_\Gamma}{\partial p\_m} \right) \delta p\_m = 0,\tag{23b}$$

$$\int\_{\Gamma} d\Gamma \beta \frac{\mu\_m^t}{\rho\_c} \left(1 - \alpha \left(\frac{\rho\_d}{\rho\_c} - 1\right)\right) (\mathbf{n} \cdot \nabla) \delta \mathbf{a} = 0,\tag{23c}$$

where **u**<sup>t</sup> is the tangential component of the adjoint velocity, from which we deduce the boundary conditions for the adjoint variables at the inlet to be:

$$\mathbf{u}\_{\mathbb{H}} = \mathbf{0},\tag{24a}$$

$$\mu\_{\rm n} = -\frac{1}{\rho\_m} \frac{\partial f\_{\Gamma}}{\partial p\_m} \,, \tag{24b}$$

$$
\beta = 0 \iff \mu\_m^t \neq 0. \tag{24c}
$$

Note that these derivations do not impose a condition for *q*. Since *q* enters the adjoint drift flux equations in a manner similar to the way *pm* enters the primal drift flux equations, the zero gradient boundary condition of *pm* at the inlet is applied to *q* as well,

$$(\mathbf{n} \cdot \nabla)\boldsymbol{q} = 0.\tag{25}$$

#### *3.2. Adjoint Boundary Conditions at the Wall*

At a wall, typical primal conditions are zero velocity and zero gradient of the dispersed-phase volume fraction. Therefore, we have,

$$
\delta \mathbf{v}\_m = 0,\\
\delta \mathbf{v}\_m = 0 \text{ and } (\mathbf{n} \cdot \nabla) \delta \mathbf{a} = 0. \tag{26}
$$

*Fluids* **2020**, *5*, 31

The first integral in Equation (20a) and the second integral in Equation (20c) therefore go to zero and the terms in the first integral in Equation (20c), containing **v***m*, go to zero. Equation (20) therefore reduces to:

$$\int\_{\Gamma} d\Gamma \mu\_m (\mathbf{n} \cdot \nabla) \delta \mathbf{v}\_m \cdot \mathbf{u} = 0,\tag{27a}$$

$$\int\_{\Gamma} d\Gamma \left(\rho\_m u\_\mathbf{n} + \frac{\partial f\_\Gamma}{\partial p\_m}\right) \delta p\_\mathbf{n} = 0,\tag{27b}$$

$$\begin{split} &\int\_{\Gamma} d\Gamma \left( \left( \frac{\partial}{\partial a} (a \mathbf{v}\_{dj}) \cdot \mathbf{n} + \frac{\mu\_m^t}{\rho\_c} \left( \frac{\rho\_d}{\rho\_c} - 1 \right) (\mathbf{n} \cdot \nabla) a \right) \beta \right. \\ &\left. + \frac{\mu\_m^t}{\rho\_c} \left( 1 - a \left( \frac{\rho\_d}{\rho\_c} - 1 \right) \right) (\mathbf{n} \cdot \nabla) \beta + \mathbb{C}\_2 + \frac{\partial f\_{\Gamma}}{\partial a} \right) \delta a = 0. \end{split} \tag{27c}$$

As at the inlet, the primal velocity does not diverge and *δ***v***m*<sup>t</sup> = 0 along the wall, so Equation (27) reduces to:

$$\int\_{\Gamma} d\Gamma \mu\_{\rm m}(\mathbf{n} \cdot \nabla) \delta \mathbf{v}\_{\rm mt} \cdot \mathbf{u}\_{\rm t} = 0,\tag{28a}$$

$$\int\_{\Gamma} d\Gamma \left(\rho\_m \mu\_\mathbf{n} + \frac{\partial f\_\Gamma}{\partial p\_m}\right) \delta p\_m = 0,\tag{28b}$$

$$\begin{split} &\int\_{\Gamma} d\Gamma \left( \left( \frac{\partial}{\partial \boldsymbol{\alpha}} (a \mathbf{v}\_{dj}) \cdot \mathbf{n} + \frac{\mu\_m^t}{\rho\_c} \left( \frac{\rho\_d}{\rho\_c} - 1 \right) (\mathbf{n} \cdot \nabla) \boldsymbol{\alpha} \right) \boldsymbol{\beta} \right. \\ &\left. + \frac{\mu\_m^t}{\rho\_c} \left( 1 - \boldsymbol{\alpha} \left( \frac{\rho\_d}{\rho\_c} - 1 \right) \right) (\mathbf{n} \cdot \nabla) \boldsymbol{\beta} + \mathbf{C}\_2 + \frac{\partial f\_{\Gamma}}{\partial \boldsymbol{\alpha}} \right) \boldsymbol{\delta} \boldsymbol{a} = \boldsymbol{0}, \end{split} \tag{28c}$$

from which we deduce the boundary conditions for the adjoint variables at the wall to be:

$$\mathbf{u}\_{\mathbf{t}} = \mathbf{0},\tag{29a}$$

$$u\_{\mathbf{n}} = -\frac{1}{\rho\_m} \frac{\partial f\_{\Gamma}}{\partial p\_m},\tag{29b}$$

$$\begin{split} \left( \frac{\partial}{\partial \boldsymbol{\alpha}} (\boldsymbol{a} \mathbf{v}\_{dj}) \cdot \mathbf{n} + \frac{\mu\_{\rm m}^{\rm t}}{\rho\_{\rm c}} \left( \frac{\rho\_{\rm d}}{\rho\_{\rm c}} - 1 \right) (\mathbf{n} \cdot \nabla) \boldsymbol{\alpha} \right) \boldsymbol{\beta} \\ + \frac{\mu\_{\rm m}^{\rm t}}{\rho\_{\rm c}} \left( 1 - \boldsymbol{a} \left( \frac{\rho\_{\rm d}}{\rho\_{\rm c}} - 1 \right) \right) (\mathbf{n} \cdot \nabla) \boldsymbol{\beta} = -\mathsf{C}\_{2} - \frac{\partial f\_{\Gamma}}{\partial \boldsymbol{\alpha}}. \end{split} \tag{29c}$$

Equation (29c) is used to determine *β* and, as at the inlet, Equation (25) applies.

## *3.3. Adjoint Boundary Conditions at the Outlet*

At an outlet, typical primal conditions are zero pressure and zero gradient of velocity and dispersed-phase volume fraction. Therefore, we have,

$$
\delta \rho\_m = 0,\\
\text{ (n} \cdot \nabla) \delta \mathbf{v}\_m = 0 \text{ and } (\mathbf{n} \cdot \nabla) \delta \mathbf{a} = 0. \tag{30}
$$

The second integral in Equations (20a) and (20c) therefore goes to zero and, with *δpm* = 0, Equation (20b) is identically fulfilled. The remaining terms in Equation (20) are the first integrals in Equations (20a) and (20c), which can be made to go to zero by enforcing the integrands to vanish:

$$\begin{split} \mathbf{n}(\mathbf{u}\cdot\rho\_{m}\mathbf{v}\_{m}) + \rho\_{m}\mathbf{v}\_{m}\cdot\mathbf{n}\mathbf{u} + \mu\_{m}(\mathbf{n}\cdot\nabla)\mathbf{u} - q\rho\_{m}\mathbf{n} + a\beta\mathbf{n} + \frac{\partial f\_{\Gamma}}{\partial \mathbf{v}\_{m}} \\ -2\frac{\partial \mu\_{d}}{\partial \mathbf{v}\_{m}}\mathbf{n}\cdot\mathbf{D}(\mathbf{u})\cdot\mathbf{v}\_{m} = 0, \\ \left(\mathbf{v}\_{m}\cdot\mathbf{n} + \frac{\partial}{\partial \alpha}(a\mathbf{v}\_{dj})\cdot\mathbf{n} + \frac{\mu\_{m}^{t}}{\rho\_{c}}\left(\frac{\rho\_{d}}{\rho\_{c}} - 1\right)(\mathbf{n}\cdot\nabla)a\right)\beta \\ + \frac{\mu\_{m}^{t}}{\rho\_{c}}\left(1 - \alpha\left(\frac{\rho\_{d}}{\rho\_{c}} - 1\right)\right)(\mathbf{n}\cdot\nabla)\beta + \mathbf{C}\_{1} + \mathbf{C}\_{2} + \frac{\partial f\_{\Gamma}}{\partial a} = 0. \end{split} \tag{31b}$$

Note that the term containing D(**v***m*) = 0, because (**n** · ∇)**v***<sup>m</sup>* = 0. Decomposing Equation (31a) into its normal and tangential components yields:

$$\rho\_m \mathbf{u} \cdot \mathbf{v}\_m + \rho\_m u\_n \mathbf{v}\_m \cdot \mathbf{n} + \mu\_m (\mathbf{n} \cdot \nabla) u\_n - \rho\_m q + a\beta + \frac{\partial f\_\Gamma}{\partial \mathbf{v}\_m \cdot \mathbf{n}}$$

$$-2 \frac{\partial \mu\_d}{\partial \mathbf{v}\_m} \mathbf{n} \cdot \mathbf{D}(\mathbf{u}) \cdot \mathbf{v}\_m \cdot \mathbf{n} = 0,\tag{32a}$$

$$
\rho\_m \mathbf{v}\_m \cdot \mathbf{n} \mathbf{u}\_t + \mu\_m (\mathbf{n} \cdot \nabla) \mathbf{u}\_t + \frac{\partial f\_\Gamma}{\partial \mathbf{v}\_{\rm mt}} - 2 \frac{\partial \mu\_d}{\partial \mathbf{v}\_m} \mathbf{n} \cdot \mathbf{D}(\mathbf{u}) \cdot \mathbf{v}\_{\rm mt} = 0. \tag{32b}
$$

Equations (31b), (32a) and (32b) are used to determine *β*, *q* and **u**t, respectively. Since *u*<sup>n</sup> is prescribed at the inlet, the adjoint continuity equation, Equation (18b), is used to calculate *u*<sup>n</sup> at the outlet, ΦΓ. The boundary conditions for the adjoint variables at the outlet are summarised as:

$$
\rho\_m \mathbf{v}\_m \cdot \mathbf{n} \mathbf{u}\_t + \mu\_m (\mathbf{n} \cdot \nabla) \mathbf{u}\_t - 2 \frac{\partial \mu\_d}{\partial \mathbf{v}\_m} \mathbf{n} \cdot \mathbf{D}(\mathbf{u}) \cdot \mathbf{v}\_{\text{mt}} = -\frac{\partial f\_\Gamma}{\partial \mathbf{v}\_{\text{mt}}},\tag{33a}
$$

*u*<sup>n</sup> = ΦΓ, (33b)

$$\begin{split} \left( \mathbf{v}\_{m} \cdot \mathbf{n} + \frac{\partial}{\partial \mathbf{a}} (a \mathbf{v}\_{dj}) \cdot \mathbf{n} + \frac{\mu\_{m}^{t}}{\rho\_{c}} \left( \frac{\rho\_{d}}{\rho\_{c}} - 1 \right) (\mathbf{n} \cdot \nabla) \mathbf{a} \right) \boldsymbol{\beta} \\ + \frac{\mu\_{m}^{t}}{\rho\_{c}} \left( 1 - a \left( \frac{\rho\_{d}}{\rho\_{c}} - 1 \right) \right) (\mathbf{n} \cdot \nabla) \boldsymbol{\beta} = -\mathbf{C}\_{1} - \mathbf{C}\_{2} - \frac{\partial \boldsymbol{J}\_{\Gamma}}{\partial \mathbf{a}}, \end{split} \tag{33c}$$

$$\begin{split} \mathbf{u} \cdot \mathbf{v}\_m + \boldsymbol{u}\_n \mathbf{v}\_m \cdot \mathbf{n} + \nu\_m (\mathbf{n} \cdot \nabla) u\_n + \frac{a\beta}{\rho\_m} + \frac{1}{\rho\_m} \frac{\partial f\_\Gamma}{\partial \mathbf{v}\_m \cdot \mathbf{n}} \\ -\frac{2}{\rho\_m} \frac{\partial \mu\_d}{\partial \mathbf{v}\_m} \mathbf{n} \cdot \mathbf{D}(\mathbf{u}) \cdot \mathbf{v}\_m \cdot \mathbf{n} = q\_\prime \end{split} \tag{33d}$$

where *<sup>ν</sup><sup>m</sup>* <sup>=</sup> *<sup>μ</sup><sup>m</sup> ρ<sup>m</sup>* is the mixture kinematic viscosity. A summary of the boundary conditions for the adjoint variables is presented in Table 2.

**Table 2.** Adjoint boundary conditions for ducted flows.


*Fluids* **2020**, *5*, 31

#### **4. Objective Function**

The objective function is related to the dispersed-phase mass-flow rate at the boundaries of the domain,

$$J\_{\Gamma} = \alpha \rho\_d \mathbf{v}\_d \cdot \mathbf{n}\_\prime \tag{34}$$

where *αρ<sup>d</sup>* is the dispersed-phase mass fraction and **v***<sup>d</sup>* · **n** is the dispersed-phase velocity normal to the boundary. Since the phase fraction at the inlet is specified, the objective function is defined as the mass-flow rate of solid at the outlet, and Equation (12) becomes,

$$J = \int\_{\Gamma\_o} J\_{\Gamma\_o} \, d\Gamma\_{o\prime} \tag{35}$$

where *o* refers to the outlet. The derivatives of the objective function, Equation (34), with respect to the primal variables are:

$$\frac{\partial \|\Gamma\|}{\partial \mathbf{v}\_{\text{ill}} \cdot \mathbf{n}} = \mathbf{a} \rho\_{d\prime} \tag{36a}$$

$$\frac{\partial \mathbf{J}\_{\rm int}}{\partial \mathbf{v}\_{\rm int}} = \mathbf{0},\tag{36b}$$

$$\frac{\partial f\_{\Gamma}}{\partial \alpha} = \rho\_d \left( \mathbf{v}\_m + \mathbf{v}\_{d\bar{j}} + \alpha \frac{\partial \mathbf{v}\_{d\bar{j}}}{\partial \alpha} \right) \cdot \mathbf{n}\_{\prime} \tag{36c}$$

$$\frac{\partial f\_{\Gamma}}{\partial p\_{m}} = 0.\tag{36d}$$

Derivation of Equation (36c) can be found in Appendix E. Using these derivatives, the adjoint boundary conditions at an inlet reduces to:

$$\mathbf{u}\_{\mathbf{t}} = \mathbf{0},\tag{37a}$$

$$
u\_{\mathbf{n}} = \mathbf{0},\tag{37b}$$

$$
\beta = 0 \iff \mu\_m^t \neq 0,\tag{37c}
$$

$$(\mathbf{n} \cdot \nabla)q = 0.\tag{37d}$$

At a wall, there is no contribution from the objective function, so:

$$\mathbf{u}\_{\mathbf{t}} = \mathbf{0},\tag{38a}$$

$$
\mu\_{\rm n} = 0,\tag{38b}
$$

 *∂ ∂α* (*α***v***dj*) · **<sup>n</sup>** <sup>+</sup> *<sup>μ</sup><sup>t</sup> m ρc ρ<sup>d</sup> ρc* − 1 (**n** · ∇)*α β*

$$+\frac{\mu\_m^t}{\rho\_c}\left(1-\alpha\left(\frac{\rho\_d}{\rho\_c}-1\right)\right)(\mathbf{n}\cdot\nabla)\beta=0,\tag{38c}$$

$$(\mathbf{n} \cdot \nabla)q = 0.\tag{38d}$$

*Fluids* **2020**, *5*, 31

Note that *C*<sup>2</sup> = 0 because **u** = 0. At an outlet, to satisfy the adjoint continuity equation, Equation (18b), *u*<sup>n</sup> = 0, so:

$$\mathbf{v}\_{\mathrm{m}} \cdot \mathbf{n} \mathbf{u}\_{\mathrm{t}} + \nu\_{\mathrm{m}} (\mathbf{n} \cdot \nabla) \mathbf{u}\_{\mathrm{t}} - \frac{2}{\rho\_{\mathrm{m}}} \frac{\partial \mu\_{d}}{\partial \mathbf{v}\_{\mathrm{m}}} \mathbf{n} \cdot \mathbf{D}(\mathbf{u}) \cdot \mathbf{v}\_{\mathrm{mt}} = 0,\tag{39a}$$

*u*<sup>n</sup> = 0, (39b)

$$\begin{split} \left(\mathbf{v}\_{\mathrm{m}} \cdot \mathbf{n} + \frac{\partial}{\partial \boldsymbol{\alpha}} (\boldsymbol{a} \mathbf{v}\_{dj}) \cdot \mathbf{n} + \frac{\mu\_{m}^{t}}{\rho\_{c}} \left(\frac{\rho\_{d}}{\rho\_{c}} - 1\right) (\mathbf{n} \cdot \nabla) \boldsymbol{\alpha}\right) \boldsymbol{\beta} \\ + \frac{\mu\_{m}^{t}}{r} \left(1 - \boldsymbol{a} \left(\frac{\rho\_{d}}{r} - 1\right)\right) (\mathbf{n} \cdot \nabla) \boldsymbol{\beta} = -\mathbb{C}\_{1} - \frac{\partial I\_{\Gamma}}{\mathbb{S}\_{1}}, \end{split} \tag{39c}$$

$$\left(\underset{\begin{subarray}{c}\mathbf{a}\end{subarray}}{+\frac{\mu\_{m}^{\mathrm{I}}}{\rho\_{c}}}\left(1-\mathfrak{a}\left(\frac{\rho\_{d}}{\rho\_{c}}-1\right)\right)(\mathbf{n}\cdot\nabla)\beta=-\mathrm{C}\_{1}-\frac{\hat{\mathfrak{d}}f\_{\Gamma}}{\hat{\alpha}\alpha},\tag{39c}$$

$$\mathbf{u} \cdot \mathbf{v}\_{\mathrm{m}} + \nu\_{\mathrm{m}} (\mathbf{n} \cdot \nabla) u\_{\mathrm{m}} + \frac{a\beta}{\rho\_{\mathrm{m}}} + \frac{1}{\rho\_{\mathrm{m}}} \frac{\partial l\_{\Gamma}}{\partial \mathbf{v}\_{\mathrm{m}} \cdot \mathbf{n}} - \frac{2}{\rho\_{\mathrm{m}}} \frac{\partial \mu\_{d}}{\partial \mathbf{v}\_{\mathrm{m}}} \mathbf{n} \cdot \mathbf{D}(\mathbf{u}) \cdot \mathbf{v}\_{\mathrm{m}} \cdot \mathbf{n} = q. \tag{39d}$$

Note that *C*<sup>2</sup> = 0 because *u*<sup>n</sup> = 0 and D(**v***m*) = 0 when (**n** · ∇)**v***<sup>m</sup>* = 0. A summary of the adjoint boundary conditions, using the objective function defined in Equation (34), is presented in Table 3.

**Table 3.** Adjoint boundary conditions, using objective function Equation (34).


#### **5. Settling Velocity**

The equations thus far have been derived for the most general case in which the settling (drift) velocity has not been specified. Of course the settling velocity is key to the behaviour of the drift flux model, and incorporates much of the physics of the multiphase system. Here we will derive the appropriate additional equations for two common settling velocity models, vis. the Dahl [23] and Takacs [24] models.

#### *5.1. Dahl Model*

In this formulation **v***dj* is modelled using,

$$\mathbf{v}\_{dj} = \mathbf{v}\_0 \mathbf{1} \mathbf{0}^{-ka},\tag{40}$$

where **v**<sup>0</sup> is the maximum theoretical settling velocity and *k* is a settling parameter, and its partial derivative with respect to *α* is given by,

$$\frac{\partial \mathbf{v}\_{d\bar{j}}}{\partial \alpha} = -k \ln 10 \,\mathbf{v}\_{d\bar{j}}.\tag{41}$$

#### *5.2. Takacs Model*

In this formulation **v***dj* is modelled using,

$$\mathbf{v}\_{d\boldsymbol{j}} = \mathbf{v}\_0 \left( e^{-a(\boldsymbol{a} - \boldsymbol{u}\_r)} - e^{-a\_1(\boldsymbol{a} - \boldsymbol{u}\_r)} \right),\tag{42a}$$

$$\mathbf{0} \lessapprox \mathbf{v}\_{d\dot{j}} \lessapprox \mathbf{v}\_{00\prime} \tag{42b}$$

where:


and its partial derivative with respect to *α* is given by,

$$\frac{\partial \mathbf{v}\_{dj}}{\partial a} = \mathbf{v}\_0 \left( -a e^{-a(a - a\_r)} + a\_1 e^{-a\_1(a - a\_r)} \right). \tag{43}$$

For both models the partial derivative of *α***v***dj* with respect to *α* is given by,

$$\frac{\partial}{\partial \alpha} (\alpha \mathbf{v}\_{d\dot{j}}) = \mathbf{v}\_{d\dot{j}} + \alpha \frac{\partial \mathbf{v}\_{d\dot{j}}}{\partial \alpha}. \tag{44}$$

#### **6. Conclusions**

In this paper we have derived, for the first time, the adjoint equations based on the drift flux model for dispersed multiphase flow. In addition to the adjoint drift flux equations themselves we have presented the adjoint boundary conditions for the common boundary conditions (Inlet, Outlet, Wall), as well as a treatment of the generic objective function, and specific formulations corresponding to the common settling velocity models proposed by Dahl [23] and Takacs [24]. From these elements a full adjoint set of equations can be derived for any specific ducted flow problem, and of course implemented in an appropriate numerical code. This of course also presents many, largely numerical and coding, challenges, and this will be the subject of subsequent papers.

**Author Contributions:** Formal analysis, S.G.; Funding acquisition, D.S.J.; Investigation, S.G. and G.R.T.; Methodology, D.S.J.; Project administration, G.R.T.; Writing–original draft, S.G.; Writing–review and editing, D.S.J. and G.R.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the STREAM CDT under grant number EP/L015412/1 and InnovateUK under KTP grant KTP010621.

**Conflicts of Interest:** There is no conflict of interest for Hydro International Ltd. in the publication of the paper. Hydro International holds interests in the implementation and application of the methods discussed in the paper for specific scenarios and geometries, but this work will remain private and unpublished.

#### **Appendix A. Derivation of Equation (10a)**

The variation of (*R*1, *R*2, *R*3)*<sup>T</sup>* with respect to **v***<sup>m</sup>* is calculated as,

$$\begin{split} \delta\_{\mathbf{v}\_{m}}(R\_{1},R\_{2},R\_{3})^{T} &= \delta\_{\mathbf{v}\_{m}} \Big( (\mathbf{v}\_{m}\cdot\nabla)(\rho\_{m}\mathbf{v}\_{m}) + \nabla(\rho\_{m}p\_{m}) - \nabla\cdot(2\mu\_{m}\mathbf{D}(\mathbf{v}\_{m}) \Big) \\ &+ \nabla\cdot\Big(\frac{a}{1-a}\frac{\rho\_{\ell}\rho\_{\ell}}{\rho\_{m}}\mathbf{v}\_{\ell j}\mathbf{v}\_{\ell j}\Big) - \rho\_{m}\mathbf{g} + \mathbb{N}\rho\_{m}\mathbf{v}\_{m} \Big) \\ &= (\delta\mathbf{v}\_{m}\cdot\nabla)(\rho\_{m}\mathbf{v}\_{m}) + (\mathbf{v}\_{m}\cdot\nabla)(\rho\_{m}\delta\mathbf{v}\_{m}) - \nabla\cdot\Big(2\mu\_{m}\mathbf{D}(\delta\mathbf{v}\_{m})\Big) \\ &- \nabla\cdot\Big(2\delta\_{\mathbf{v}\_{m}}\mu\_{m}\mathbf{D}(\mathbf{v}\_{m})\Big) + \mathbb{N}\rho\_{m}\delta\mathbf{v}\_{m}. \end{split} \tag{A1}$$

As stated above, *μ<sup>m</sup>* is defined as the sum of the continuum, dispersed-phase and mixture turbulent viscosities,

$$
\mu\_m = \mu\_c + \mu\_d + \mu\_{m'}^t \tag{A2}
$$

where *μ<sup>c</sup>* is constant, *μ<sup>d</sup>* is a function of **v***<sup>m</sup>* and *α*, and *μ<sup>t</sup> <sup>m</sup>* is obtained from turbulence modelling. Equation (A1) can now be rewritten as,

$$\begin{split} \delta\_{\mathbf{v}\_{m}} (R\_{1}, R\_{2}, R\_{3})^{T} &= (\delta \mathbf{v}\_{m} \cdot \nabla) (\rho\_{m} \mathbf{v}\_{m}) + (\mathbf{v}\_{m} \cdot \nabla)(\rho\_{m} \delta \mathbf{v}\_{m}) - \nabla \cdot \left( 2 \mu\_{m} \mathbf{D} (\delta \mathbf{v}\_{m}) \right) \\ &- \nabla \cdot \left( 2 \delta\_{\mathbf{v}\_{m}} \mu\_{d} \mathbf{D} (\mathbf{v}\_{m}) \right) + \aleph \rho\_{m} \delta \mathbf{v}\_{m} \end{split} \tag{A3}$$

where *δ***v***<sup>m</sup> μ<sup>t</sup> <sup>m</sup>* has been neglected.

#### **Appendix B. Derivation of Equation (10g)**

The variation of (*R*1, *R*2, *R*3)*<sup>T</sup>* with respect to *α* is calculated as,

$$\delta\_{\mathfrak{a}} (R\_1, R\_2, R\_3)^T = \delta\_{\mathfrak{a}} \left( (\mathbf{v}\_m \cdot \nabla)(\rho\_m \mathbf{v}\_m) + \nabla(\rho\_m p\_m) - \nabla \cdot (2\mu\_m \mathbf{D}(\mathbf{v}\_m)) \right)$$

$$\begin{split} &+ \nabla \cdot \left( \frac{\mathbf{a}}{1-\mathbf{a}} \frac{\rho\_{\mathfrak{f}} \varrho\_{\mathbf{f}}}{\rho\_m} \mathbf{v}\_{\operatorname{d}\mathbf{f}} \mathbf{v}\_{\operatorname{d}\mathbf{f}} \right) - \rho\_{\mathfrak{m}} \mathbf{g} + \aleph\_{\mathfrak{M}} \mathbf{v}\_{\mathfrak{m}} \Big) \\ &= \delta\_{\mathfrak{a}} \left( (\mathbf{v}\_m \cdot \nabla)(\rho\_m \mathbf{v}\_m) \right) + \delta\_{\mathfrak{a}} \nabla (\rho\_m p\_m) - \delta\_{\mathfrak{a}} \nabla \cdot \left( 2\mu\_m \mathbf{D}(\mathbf{v}\_m) \right) \\ &+ \delta\_{\mathfrak{a}} \nabla \cdot \mathbf{A} + \delta\_{\mathfrak{a}} \left( \rho\_m (\mathbf{N} \mathbf{v}\_m - \mathbf{g}) \right), \end{split} \tag{A4}$$

where

$$\mathbf{A} = \frac{\boldsymbol{\alpha}}{1 - \boldsymbol{\alpha}} \frac{\rho\_{\ell} \rho\_{d}}{\rho\_{m}} \mathbf{v}\_{dj} \mathbf{v}\_{dj} \tag{A5}$$

and

$$\begin{split} \rho\_m &= \mathfrak{a}\rho\_d + (1 - \mathfrak{a})\rho\_c \\ &= (1 - \mathfrak{a})\rho\_\mathbb{C} \left( 1 + \frac{\mathfrak{a}}{1 - \mathfrak{a}} \frac{\rho\_d}{\rho\_c} \right) . \end{split} \tag{A6}$$

Substituting Equation (A6) into Equation (A5) and rewriting the parentheses as binomial expansions,

$$\mathbf{A} = a(1 - \boldsymbol{\alpha})^{-1} \rho\_c \rho\_d \frac{(1 - \boldsymbol{\alpha})^{-1}}{\rho\_c} \left( 1 + \frac{\boldsymbol{\alpha}}{1 - \boldsymbol{\alpha}} \frac{\rho\_d}{\rho\_c} \right)^{-1} \mathbf{v}\_{dj} \mathbf{v}\_{dj}$$

$$= a \rho\_d \left( 1 + \boldsymbol{\alpha} \left( 2 - \frac{1}{1 - \boldsymbol{\alpha}} \frac{\rho\_d}{\rho\_c} \right) + \cdots \right) \mathbf{v}\_{dj} \mathbf{v}\_{dj}.\tag{A7}$$

As *α* 1 and *ρ<sup>d</sup>* ≈ 2*ρ<sup>c</sup>* =⇒ <sup>2</sup> <sup>−</sup> <sup>1</sup> 1−*α ρd ρc* <sup>&</sup>lt; 1 and ignoring terms containing squared and higher powers of *α*,

$$\mathbf{A} \approx \alpha \rho\_d \mathbf{v}\_{dj} \mathbf{v}\_{dj}.\tag{A8}$$

Substituting Equations (A8) and (A2) into Equation (A4),

$$\begin{split} \delta\_{a} (R\_{1}, R\_{2}, R\_{3})^{T} \approx & (\rho\_{d} - \rho\_{c}) \left( (\mathbf{v}\_{m} \cdot \nabla)(\delta a \mathbf{v}\_{m}) + \nabla(\delta a p\_{m}) + \delta a (\mathbb{N} \mathbf{v}\_{m} - \mathbf{g}) \right) \\ & - \nabla \cdot \left( 2 \delta\_{d} \mu\_{d} \mathbf{D} (\mathbf{v}\_{m}) \right) + \nabla \cdot \delta\_{a} (a \rho\_{d} \mathbf{v}\_{d\uparrow} \mathbf{v}\_{d\uparrow}), \end{split} \tag{A9}$$

where *δαμ<sup>t</sup> <sup>m</sup>* has been neglected.

## **Appendix C. Derivation of Equation (10i)**

Similarly, the variation of *R*<sup>5</sup> with respect to *α* is calculated as,

$$\delta\_a R \mathfrak{z} = \delta\_a \left( \nabla \cdot (a \mathbf{v}\_m) + \nabla \cdot \left( \frac{a \rho\_c}{\rho\_m} \mathbf{v}\_{dj} \right) - \nabla \cdot (K \nabla a) \right)$$

$$= \delta\_a \nabla \cdot (a \mathbf{v}\_m) + \delta\_a \nabla \cdot \mathbf{B} - \delta\_a \nabla \cdot (K \nabla a),\tag{A10}$$

where

$$\mathbf{B} = \frac{\kappa \rho\_c}{\rho\_m} \mathbf{v}\_{dj}.\tag{A11}$$

*Fluids* **2020**, *5*, 31

Substituting Equation (A6) into Equation (A11) and rewriting the parentheses as binomial expansions,

$$\begin{split} \mathbf{B} &= \kappa \rho\_c \frac{(1-a)^{-1}}{\rho\_c} \left( 1 + \frac{a}{1-a} \frac{\rho\_d}{\rho\_c} \right)^{-1} \mathbf{v}\_{dj} \\ &= \kappa \left( 1 + a \left( 1 - \frac{1}{1-a} \frac{\rho\_d}{\rho\_c} \right) + \cdots \right) \mathbf{v}\_{dj}. \end{split} \tag{A12}$$

As *α* 1 and *ρ<sup>d</sup>* ≈ 2*ρ<sup>c</sup>* =⇒ <sup>1</sup> <sup>−</sup> <sup>1</sup> 1−*α ρd ρc* <sup>&</sup>lt; 2 and ignoring terms containing squared and higher powers of *α*,

$$\mathbf{B} \approx \boldsymbol{\kappa} \mathbf{v}\_{dj}.\tag{A13}$$

As stated above, *K* is defined as the mixture eddy diffusivity,

$$\nu\_m^t = \frac{\mu\_m^t}{\rho\_m}.\tag{A14}$$

From Equation (A6), <sup>1</sup> *<sup>ρ</sup><sup>m</sup>* can be written as,

$$\begin{split} \frac{1}{\rho\_m} &= \frac{1}{\rho\_c} (1 - a)^{-1} \left( 1 + \frac{a}{1 - a} \frac{\rho\_d}{\rho\_c} \right)^{-1} \\ &= \frac{1}{\rho\_c} \left( 1 + a \left( 1 - \frac{\rho\_d}{\rho\_c} \right) \right), \end{split} \tag{A15}$$

ignoring terms containing squared and higher powers of *α*.

Substituting Equations (A14) and (A15) into the term containing *K* in Equation (A10),

$$\begin{split} \delta\_a \nabla \cdot (K \nabla a) &= \delta\_a \nabla \cdot \left( \frac{\mu\_m^t}{\rho\_m} \nabla a \right) \\ &= \delta\_a \nabla \cdot \left( \frac{\mu\_m^t}{\rho\_c} \left( 1 + a \left( 1 - \frac{\rho\_d}{\rho\_c} \right) \right) \nabla a \right) \\ &= \nabla \cdot \left( \frac{\mu\_m^t}{\rho\_c} \left( 1 + (a + \delta a) \left( 1 - \frac{\rho\_d}{\rho\_c} \right) \right) \nabla (a + \delta a) \right) \\ &\quad - \nabla \cdot \left( \frac{\mu\_m^t}{\rho\_c} \left( 1 + a \left( 1 - \frac{\rho\_d}{\rho\_c} \right) \right) \nabla a \right) \\ &= \nabla \cdot \left( \frac{\mu\_m^t}{\rho\_c} \delta a \left( 1 - \frac{\rho\_d}{\rho\_c} \right) \nabla a \right) + \nabla \cdot \left( \frac{\mu\_m^t}{\rho\_c} \nabla \delta a \right) \\ &\quad + \nabla \cdot \left( \frac{\mu\_m^t}{\rho\_c} a \left( 1 - \frac{\rho\_d}{\rho\_c} \right) \nabla \delta a \right), \end{split} \tag{A16}$$

ignoring the term containing *δα*∇*δα*, because when substitued into Equation (9) becomes terms containing squared powers of *δα*. Equation (A10) now becomes,

$$\begin{split} \delta\_{\mathsf{a}} R\_{\mathsf{F}} & \approx \nabla \cdot (\delta a \mathbf{v}\_{m}) + \nabla \cdot \delta\_{\mathsf{a}} (a \mathbf{v}\_{d\bar{j}}) - \nabla \cdot \left(\frac{\mu\_{m}^{t}}{\rho\_{\mathsf{c}}} \nabla \delta a\right) \\ & + \frac{\mu\_{m}^{t}}{\rho\_{\mathsf{c}}} \left(\frac{\rho\_{d}}{\rho\_{\mathsf{c}}} - 1\right) \nabla \cdot (\delta a \nabla a) + \frac{\mu\_{m}^{t}}{\rho\_{\mathsf{c}}} \left(\frac{\rho\_{d}}{\rho\_{\mathsf{c}}} - 1\right) \nabla \cdot (a \nabla \delta a), \end{split} \tag{A17}$$

where *δαμ<sup>t</sup> <sup>m</sup>* has been neglected.

#### **Appendix D. Derivation of Equation (13)**

Decomposing the objective function into contributions from the boundary and interior of the domain, according to Equation (12), the terms in Equation (11) can be written as follows. The variations of the objective function can be written as,

$$
\delta\boldsymbol{\delta}\_{\mathbf{v}}\boldsymbol{J} = \int\_{\Gamma} d\Gamma \, \frac{\partial \boldsymbol{f}\_{\Gamma}}{\partial \mathbf{v}\_{\boldsymbol{m}}} \cdot \delta\mathbf{v}\_{\boldsymbol{m}} + \int\_{\Omega} d\Omega \, \frac{\partial \boldsymbol{f}\_{\Omega}}{\partial \mathbf{v}\_{\boldsymbol{m}}} \cdot \delta\mathbf{v}\_{\boldsymbol{m}} \tag{A18}
$$

$$\delta\_{p\_m} J = \int\_{\Gamma} d\Gamma \, \frac{\partial f\_{\Gamma}}{\partial p\_m} \delta p\_m + \int\_{\Omega} d\Omega \, \frac{\partial f\_{\Omega}}{\partial p\_m} \delta p\_m \tag{A19}$$

and

$$
\delta\mathfrak{sl}f = \int\_{\Gamma} d\Gamma \, \frac{\partial f\_{\Gamma}}{\partial a} \delta a + \int\_{\Omega} d\Omega \, \frac{\partial f\_{\Omega}}{\partial a} \delta a. \tag{A20}
$$

Applying the product rule, divergence theorem and continuity equation, and using the Einstein notation for clarity, the terms containing **u**, **v***<sup>m</sup>* and ∇ can be written as,

$$\begin{split} \int\_{\Omega} d\Omega \,\mathbf{u} \cdot (\delta \mathbf{v}\_{\mathrm{m}} \cdot \nabla) (\rho\_{\mathrm{m}} \mathbf{v}\_{\mathrm{m}}) &= \int\_{\Omega} d\Omega \, u\_{k} \delta v\_{\mathrm{m}} \frac{\partial}{\partial x\_{i}} (\rho\_{\mathrm{m}} \mathbf{v}\_{\mathrm{mk}}) \\ &= \int\_{\Omega} d\Omega \, \frac{\partial}{\partial x\_{i}} (u\_{k} \rho\_{\mathrm{m}} \mathbf{v}\_{\mathrm{m}} \delta \mathbf{v}\_{\mathrm{m}}) - \int\_{\Omega} d\Omega \, \rho\_{\mathrm{m}} \mathbf{v}\_{\mathrm{mk}} \frac{\partial (u\_{k} \delta v\_{\mathrm{m}})}{\partial x\_{i}} \\ &= \int\_{\Gamma} d\Gamma \, \mathbf{n}\_{i} u\_{k} \rho\_{\mathrm{m}} \mathbf{v}\_{\mathrm{mk}} \delta v\_{\mathrm{m}i} - \int\_{\Omega} d\Omega \, \rho\_{\mathrm{m}} \mathbf{v}\_{\mathrm{mk}} \delta \mathbf{v}\_{\mathrm{m}} \frac{\partial u\_{k}}{\partial x\_{i}} \\ &\qquad - \int\_{\Omega} d\Omega \, \rho\_{\mathrm{m}} \mathbf{v}\_{\mathrm{mk}} u\_{k} \frac{\partial \delta v\_{\mathrm{m}}}{\partial x\_{i}} \\ &= \int\_{\Gamma} d\Gamma \, \mathbf{n} (\mathbf{u} \cdot \rho\_{\mathrm{m}} \mathbf{v}\_{\mathrm{m}}) \cdot \delta \mathbf{v}\_{\mathrm{m}} - \int\_{\Omega} d\Omega \, \nabla \mathbf{u} \cdot (\rho\_{\mathrm{m}} \mathbf{v}\_{\mathrm{m}}) \cdot \delta \mathbf{v}\_{\mathrm{m}}. \end{split} \tag{A21}$$

$$\begin{split} \int\_{\Omega} d\Omega \,\mathbf{u} \cdot (\mathbf{v}\_{m} \cdot \nabla) \rho\_{m} \delta \mathbf{v}\_{m} &= \int\_{\Omega} d\Omega \, u\_{k} v\_{mi} \frac{\partial}{\partial x\_{i}} (\rho\_{mk} \delta v\_{mk}) \\ &= \int\_{\Omega} d\Omega \, \frac{\partial}{\partial x\_{i}} (u\_{k} v\_{mi} \rho\_{mk} \delta v\_{mk}) - \int\_{\Omega} d\Omega \, \rho\_{mk} \delta v\_{mk} \frac{\partial}{\partial x\_{i}} (u\_{k} v\_{mi}) \\ &= \int\_{\Gamma} d\Gamma \, n\_{i} u\_{k} v\_{mi} \rho\_{mk} \delta v\_{mk} - \int\_{\Omega} d\Omega \, \rho\_{mk} \delta v\_{mk} v\_{mi} \frac{\partial u\_{k}}{\partial x\_{i}} \\ &\quad - \int\_{\Omega} d\Omega \, \rho\_{mk} \delta v\_{mk} u\_{k} \frac{\partial v\_{mi}}{\partial x\_{i}} \\ &= \int\_{\Gamma} d\Gamma \, \mathbf{u} \left( \rho\_{m} \mathbf{v}\_{m} \cdot \mathbf{n} \right) \cdot \delta \mathbf{v}\_{m} - \int\_{\Omega} d\Omega \, (\rho\_{m} \mathbf{v}\_{m} \cdot \nabla) \mathbf{u} \cdot \delta \mathbf{v}\_{m} \end{split} \tag{A22}$$

and

$$\begin{split} \int\_{\Omega} d\Omega \,\mathbf{u} \cdot (\mathbf{v}\_{\text{m}} \cdot \nabla) (\delta \mathbf{a} \mathbf{v}\_{\text{m}}) &= \int\_{\Omega} d\Omega \, u\_{k} v\_{\text{m}} \frac{\partial}{\partial x\_{i}} (\delta a\_{k} v\_{\text{m}}) \\ &= \int\_{\Omega} d\Omega \, \frac{\partial}{\partial x\_{i}} (u\_{k} v\_{\text{m}i} \delta a\_{k} v\_{\text{m}k}) - \int\_{\Omega} d\Omega \, \delta a\_{k} v\_{\text{m}k} \frac{\partial}{\partial x\_{i}} (u\_{k} v\_{\text{m}i}) \\ &= \int\_{\Gamma} d\Gamma \, \mathbf{n}\_{i} u\_{k} v\_{\text{m}i} \delta a\_{k} v\_{\text{m}k} - \int\_{\Omega} d\Omega \, \delta a\_{k} v\_{\text{m}k} \boldsymbol{v}\_{\text{m}i} \frac{\partial \boldsymbol{u}\_{k}}{\partial x\_{i}} \\ &\quad - \int\_{\Omega} d\Omega \, \delta a\_{k} v\_{\text{m}k} \frac{\partial \boldsymbol{v}\_{\text{m}i}}{\partial x\_{i}} \\ &= \int\_{\Gamma} d\Gamma \, \mathbf{u} (\mathbf{v}\_{\text{m}} \cdot \mathbf{n}) \cdot \mathbf{v}\_{\text{m}} \delta \mathbf{a} - \int\_{\Omega} d\Omega \, (\mathbf{v}\_{\text{m}} \cdot \nabla) \mathbf{u} \cdot \mathbf{v}\_{\text{m}} \delta \mathbf{a}. \end{split} \tag{A23}$$

*Fluids* **2020**, *5*, 31

Applying the tensor-vector identity [25], the divergence theorem and a property of the colon product, demonstrated below,

$$\begin{split} \nabla \mathbf{u} : \mathcal{D}(\delta \mathbf{v}\_{\mathrm{m}}) &= \nabla \mathbf{u} : \frac{1}{2} \left( \nabla \delta \mathbf{v}\_{\mathrm{m}} + (\nabla \delta \mathbf{v}\_{\mathrm{m}})^{T} \right) \\ &= \frac{1}{2} \left( \nabla \mathbf{u} : \nabla \delta \mathbf{v}\_{\mathrm{m}} + \nabla \mathbf{u} : (\nabla \delta \mathbf{v}\_{\mathrm{m}})^{T} \right) \\ &= \frac{1}{2} \left( \nabla \mathbf{u} : \nabla \delta \mathbf{v}\_{\mathrm{m}} + (\nabla \mathbf{u})^{T} : \nabla \delta \mathbf{v}\_{\mathrm{m}} \right) \\ &= \frac{1}{2} \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^{T} \right) : \nabla \delta \mathbf{v}\_{\mathrm{m}} \\ &= \mathbf{D}(\mathbf{u}) : \nabla \delta \mathbf{v}\_{\mathrm{m}}. \end{split} \tag{A24}$$

the term containing *μ<sup>m</sup>* can be written as,

$$\begin{split} \int\_{\Omega} d\Omega \,\mathbf{u} \cdot \nabla \cdot \left( 2\mu\_{m} \mathbf{D}(\delta \mathbf{v}\_{m}) \right) &= \int\_{\Omega} d\Omega \,\nabla \cdot \left( 2\mu\_{m} \mathbf{D}(\delta \mathbf{v}\_{m}) \cdot \mathbf{u} \right) - \int\_{\Omega} d\Omega \,\nabla \mathbf{u} : 2\mu\_{m} \mathbf{D}(\delta \mathbf{v}\_{m}) \\ &= \int\_{\Gamma} d\Gamma \, 2\mu\_{m} \mathbf{n} \cdot \mathbf{D}(\delta \mathbf{v}\_{m}) \cdot \mathbf{u} - \int\_{\Omega} d\Omega \, 2\mu\_{m} \mathbf{D}(\mathbf{u}) : \nabla \delta \mathbf{v}\_{m} \\ &= \int\_{\Gamma} d\Gamma \, 2\mu\_{m} \mathbf{n} \cdot \mathbf{D}(\delta \mathbf{v}\_{m}) \cdot \mathbf{u} - \int\_{\Omega} d\Omega \, \nabla \cdot \left( 2\mu\_{m} \mathbf{D}(\mathbf{u}) \cdot \delta \mathbf{v}\_{m} \right) \\ &\quad + \int\_{\Omega} d\Omega \, \nabla \cdot \left( 2\mu\_{m} \mathbf{D}(\mathbf{u}) \right) \cdot \delta \mathbf{v}\_{m} \\ &= \int\_{\Gamma} d\Gamma \, 2\mu\_{m} \mathbf{n} \cdot \mathbf{D}(\delta \mathbf{v}\_{m}) \cdot \mathbf{u} - \int\_{\Gamma} d\Gamma \, 2\mu\_{m} \mathbf{n} \cdot \mathbf{D}(\mathbf{u}) \cdot \delta \mathbf{v}\_{m} \\ &\quad + \int\_{\Omega} d\Omega \, \nabla \cdot \left( 2\mu\_{m} \mathbf{D}(\mathbf{u}) \right) \cdot \delta \mathbf{v}\_{m}. \end{split} \tag{A25}$$

Similarly, the terms containing *μ<sup>d</sup>* can be written as,

$$\begin{split} \int\_{\Omega} d\Omega \,\mathbf{u} \cdot \nabla \cdot \left( \delta\_{\mathbf{u}} \mu\_{d} \mathbf{D} (\mathbf{v}\_{m}) \right) &= \int\_{\Gamma} d\Gamma \,\delta\_{\mathbf{u}} \mu\_{d} \mathbf{n} \cdot \mathbf{D} (\mathbf{v}\_{m}) \cdot \mathbf{u} - \int\_{\Gamma} d\Gamma \,\delta\_{\mathbf{u}} \mu\_{d} \mathbf{n} \cdot \mathbf{D} (\mathbf{u}) \cdot \mathbf{v}\_{m} \\ &+ \int\_{\Omega} d\Omega \,\nabla \cdot \left( \delta\_{\mathbf{u}} \mu\_{d} \mathbf{D} (\mathbf{u}) \right) \cdot \mathbf{v}\_{m} \end{split} \tag{A26}$$

and

$$\begin{split} \int\_{\Omega} d\Omega \, \mathbf{u} \cdot \nabla \cdot \left( \delta\_{\mathbf{v}\_{m}} \mu\_{d} \mathrm{D}(\mathbf{v}\_{m}) \right) &= \int\_{\Gamma} d\Gamma \, \delta\_{\mathbf{v}\_{m}} \mu\_{d} \mathbf{n} \cdot \mathrm{D}(\mathbf{v}\_{m}) \cdot \mathbf{u} - \int\_{\Gamma} d\Gamma \, \delta\_{\mathbf{v}\_{m}} \mu\_{d} \mathbf{n} \cdot \mathrm{D}(\mathbf{u}) \cdot \mathbf{v}\_{m} \\ &+ \int\_{\Omega} d\Omega \, \nabla \cdot \left( \delta\_{\mathbf{v}\_{m}} \mu\_{d} \mathrm{D}(\mathbf{u}) \right) \cdot \mathbf{v}\_{m}. \end{split} \tag{A27}$$

Applying the product rule and divergence theorem, the remaining terms in Equation (11) can be written as,

$$\begin{split} \int\_{\Omega} d\Omega \, q \nabla \cdot (\rho\_{\text{m}} \delta \mathbf{v}\_{\text{m}}) &= \int\_{\Omega} d\Omega \, \nabla \cdot (q \rho\_{\text{m}} \delta \mathbf{v}\_{\text{m}}) - \int\_{\Omega} d\Omega \, \nabla q \cdot \rho\_{\text{m}} \delta \mathbf{v}\_{\text{m}} \\ &= \int\_{\Gamma} d\Gamma \, q \rho\_{\text{m}} \mathbf{n} \cdot \delta \mathbf{v}\_{\text{m}} - \int\_{\Omega} d\Omega \, \rho\_{\text{m}} \nabla q \cdot \delta \mathbf{v}\_{\text{m}} \end{split} \tag{A28}$$

$$\begin{split} \int\_{\Omega} d\Omega \,\beta \nabla \cdot (a \delta \mathbf{v}\_m) &= \int\_{\Omega} d\Omega \,\nabla \cdot (\beta a \delta \mathbf{v}\_m) - \int\_{\Omega} d\Omega \,\nabla \beta \cdot a \delta \mathbf{v}\_m \\ &= \int\_{\Gamma} d\Gamma \, a \beta \mathbf{n} \cdot \delta \mathbf{v}\_m - \int\_{\Omega} d\Omega \, a \nabla \beta \cdot \delta \mathbf{v}\_m \end{split} \tag{A29}$$

*Fluids* **2020**, *5*, 31

$$\begin{split} \int\_{\Omega} d\Omega \,\mathbf{u} \cdot \nabla (\rho\_m \delta p\_m) &= \int\_{\Omega} d\Omega \,\nabla \cdot (\mathbf{u} \rho\_m \delta p\_m) - \int\_{\Omega} d\Omega \,\nabla \cdot \mathbf{u} \rho\_m \delta p\_m \\ &= \int\_{\Gamma} d\Gamma \,\rho\_m \mathbf{u} \cdot \mathbf{n} \delta p\_m - \int\_{\Omega} d\Omega \,\nabla \cdot \rho\_m \mathbf{u} \delta p\_m \end{split} \tag{A30}$$

$$\begin{split} \int\_{\Omega} d\Omega \,\mathbf{u} \cdot \nabla (\delta a p\_m) &= \int\_{\Omega} d\Omega \,\nabla \cdot (\mathbf{u} \delta a p\_m) - \int\_{\Omega} d\Omega \,\nabla \cdot \mathbf{u} \delta a p\_m \\ &= \int\_{\Gamma} d\Gamma \,\mathbf{u} \cdot \mathbf{n} \delta a p\_m - \int\_{\Omega} d\Omega \,\nabla \cdot \mathbf{u} \delta a p\_m \end{split} \tag{A31}$$

$$\begin{split} \int\_{\Omega} d\Omega \,\mathbf{u} \cdot \nabla \cdot \boldsymbol{\delta}\_{a} (\boldsymbol{a}\rho\_{d} \mathbf{v}\_{d\boldsymbol{j}} \mathbf{v}\_{d\boldsymbol{j}}) &= \int\_{\Omega} d\Omega \,\nabla \cdot \left( \mathbf{u} \boldsymbol{\delta}\_{a} (\boldsymbol{a}\rho\_{d} \mathbf{v}\_{d\boldsymbol{j}} \mathbf{v}\_{d\boldsymbol{j}}) \right) \\ &- \int\_{\Omega} d\Omega \,\nabla \cdot \mathbf{u} \,\boldsymbol{\delta}\_{a} (\boldsymbol{a}\rho\_{d} \mathbf{v}\_{d\boldsymbol{j}} \mathbf{v}\_{d\boldsymbol{j}}) \\ &= \int\_{\Gamma} d\Gamma \,\mathbf{u} \cdot \mathbf{n} \boldsymbol{\delta}\_{a} (\boldsymbol{a}\rho\_{d} \mathbf{v}\_{d\boldsymbol{j}} \mathbf{v}\_{d\boldsymbol{j}}) \\ &- \int\_{\Omega} d\Omega \,\nabla \cdot \mathbf{u} \,\boldsymbol{\delta}\_{a} (\boldsymbol{a}\rho\_{d} \mathbf{v}\_{d\boldsymbol{j}} \mathbf{v}\_{d\boldsymbol{j}}), \end{split} \tag{A32}$$

$$\begin{split} \int\_{\Omega} d\Omega \, q \nabla \cdot (\delta a \mathbf{v}\_m) &= \int\_{\Omega} d\Omega \, \nabla \cdot (q \delta a \mathbf{v}\_m) - \int\_{\Omega} d\Omega \, \nabla q \cdot (\delta a \mathbf{v}\_m) \\ &= \int\_{\Gamma} d\Gamma \, q \mathbf{v}\_m \cdot \mathbf{n} \delta a - \int\_{\Omega} d\Omega \, (\mathbf{v}\_m \cdot \nabla) q \delta a, \end{split} \tag{A33}$$

$$\begin{split} \int\_{\Omega} d\Omega \,\beta \nabla \cdot (\delta a \mathbf{v}\_{m}) &= \int\_{\Omega} d\Omega \,\nabla \cdot (\beta \delta a \mathbf{v}\_{m}) - \int\_{\Omega} d\Omega \,\nabla \beta \cdot (\delta a \mathbf{v}\_{m}) \\ &= \int\_{\Gamma} d\Gamma \,\beta \mathbf{v}\_{m} \cdot \mathbf{n} \delta a - \int\_{\Omega} d\Omega \,(\mathbf{v}\_{m} \cdot \nabla) \beta \delta a, \end{split} \tag{A34}$$

$$\begin{split} \int\_{\Omega} d\Omega \,\boldsymbol{\beta} \nabla \cdot \boldsymbol{\delta}\_{\boldsymbol{a}} (\boldsymbol{a} \mathbf{v}\_{d\boldsymbol{j}}) &= \int\_{\Omega} d\Omega \,\nabla \cdot \left( \boldsymbol{\beta} \delta\_{\boldsymbol{a}} (\boldsymbol{a} \mathbf{v}\_{d\boldsymbol{j}}) \right) - \int\_{\Omega} d\Omega \,\nabla \beta \cdot \boldsymbol{\delta}\_{\boldsymbol{a}} (\boldsymbol{a} \mathbf{v}\_{d\boldsymbol{j}}) \\ &= \int\_{\Gamma} d\Gamma \,\boldsymbol{\beta} \delta\_{\boldsymbol{a}} (\boldsymbol{a} \mathbf{v}\_{d\boldsymbol{j}}) \cdot \mathbf{n} - \int\_{\Omega} d\Omega \,\boldsymbol{\delta}\_{\boldsymbol{a}} (\boldsymbol{a} \mathbf{v}\_{d\boldsymbol{j}}) \cdot \nabla \beta, \end{split} \tag{A35}$$

$$\begin{split} \int\_{\Omega} d\Omega \,\beta \frac{\mu\_{\text{in}}^{\text{in}}}{\rho\_{c}} \nabla \cdot \nabla \delta \boldsymbol{a} &= \frac{\mu\_{\text{in}}^{\text{in}}}{\rho\_{c}} \left( \int\_{\Omega} d\Omega \,\nabla \cdot (\beta \nabla \delta \boldsymbol{a}) - \int\_{\Omega} d\Omega \,\nabla \beta \cdot \nabla \delta \boldsymbol{a} \right) \\ &= \frac{\mu\_{\text{in}}^{\text{in}}}{\rho\_{c}} \left( \int\_{\Gamma} d\Gamma \,\beta \mathbf{n} \cdot \nabla \delta \boldsymbol{a} \right. \end{split} \tag{A36}$$
 
$$\begin{split} -\int\_{\Omega} d\Omega \,\nabla \cdot (\delta \boldsymbol{a} \nabla \beta) &+ \int\_{\Omega} d\Omega \,\delta \boldsymbol{a} \nabla \cdot \nabla \beta \\ &= \frac{\mu\_{\text{in}}^{\text{in}}}{\rho\_{c}} \left( \int\_{\Gamma} d\Gamma \,\beta (\mathbf{n} \cdot \nabla) \delta \boldsymbol{a} \right. \end{split} \tag{A36}$$
 
$$\begin{split} -\int\_{\Gamma} d\Gamma \,\delta \boldsymbol{a} (\mathbf{n} \cdot \nabla) \beta &+ \int\_{\Omega} d\Omega \,\delta \boldsymbol{a} \nabla \cdot \nabla \beta \end{split} \tag{A37}$$

$$\begin{split} \int\_{\Omega} d\Omega \,\beta \frac{\mu\_m^t}{\rho\_c} \left(\frac{\rho\_d}{\rho\_c} - 1\right) \nabla \cdot (\delta a \nabla a) &= \frac{\mu\_m^t}{\rho\_c} \left(\frac{\rho\_d}{\rho\_c} - 1\right) \left(\int\_{\Omega} d\Omega \,\nabla \cdot (\delta a \nabla a) \right. \\ &\quad \left. - \int\_{\Omega} d\Omega \,\nabla \beta \cdot (\delta a \nabla a) \right) \\ &= \frac{\mu\_m^t}{\rho\_c} \left(\frac{\rho\_d}{\rho\_c} - 1\right) \left(\int\_{\Gamma} d\Gamma \,\beta \delta a (\mathbf{n} \cdot \nabla) a \\ &\quad - \int\_{\Omega} d\Omega \,\delta a \nabla a \cdot \nabla \beta \right) \end{split} \tag{A37}$$

and

$$\begin{split} \int\_{\Omega} d\Omega \beta \frac{\mu\_{\text{m}}^{\text{f}}}{\rho\_{\text{c}}} \left(\frac{\rho\_{d}}{\rho\_{\text{c}}} - 1\right) \nabla \cdot (a \nabla \delta a) &= \frac{\mu\_{\text{m}}^{\text{f}}}{\rho\_{\text{c}}} \left(\frac{\rho\_{d}}{\rho\_{\text{c}}} - 1\right) \left(\int\_{\Omega} d\Omega \, \nabla \cdot (\beta a \nabla \delta a) \right. \\ &\quad \left. - \int\_{\Omega} d\Omega \, \nabla \beta \cdot (a \nabla \delta a) \right) \\ &= \frac{\mu\_{\text{m}}^{\text{f}}}{\rho\_{\text{c}}} \left(\frac{\rho\_{d}}{\rho\_{\text{c}}} - 1\right) \left(\int\_{\Gamma} d\Gamma \, \beta a \mathbf{n} \cdot \nabla \delta a \right. \\ &\quad \left. - \int\_{\Omega} d\Omega \, \nabla \cdot (a \delta a \nabla \beta) + \int\_{\Omega} d\Omega \delta a \nabla \cdot (a \nabla \beta) \right) \\ &= \frac{\mu\_{\text{m}}^{\text{f}}}{\rho\_{\text{c}}} \left(\frac{\rho\_{d}}{\rho\_{\text{c}}} - 1\right) \left(\int\_{\Gamma} d\Gamma \, a \beta (\mathbf{n} \cdot \nabla) \delta a \right. \\ &\quad \left. - \int\_{\Gamma} d\Gamma \, a \delta a (\mathbf{n} \cdot \nabla) \beta + \int\_{\Omega} d\Omega \delta a \nabla \cdot (a \nabla \beta) \right). \end{split} \tag{A}$$

Equation (11) can now be reformulated and rearranged as Equation (13).

#### **Appendix E. Derivation of Equation (36c)**

Decomposing the dispersed-phase velocity into the mixture and dispersed-phase diffusion velocities, Equation (34) can be rewritten as,

$$\begin{split} \|f\|\_{\Gamma} &= \mathfrak{a}\rho\_d \left(\mathbf{v}\_m + \mathbf{v}\_{dm}\right) \cdot \mathbf{n} \\ &= \mathfrak{a}\rho\_d \left(\mathbf{v}\_m + \frac{\rho\_c}{\rho\_m}\mathbf{v}\_{dj}\right) \cdot \mathbf{n}\_\prime \end{split} \tag{A39}$$

where **v***dj* is defined in terms of **v***dm*, the dispersed-phase velocity relative to the mixture velocity, as **<sup>v</sup>***dj* <sup>=</sup> *<sup>ρ</sup><sup>m</sup> <sup>ρ</sup><sup>c</sup>* **v***dm*. Substituting Equation (A6) into Equation (A39) and rewriting the parentheses as binomial expansions,

$$f\_{\Gamma} = a\rho\_d \left(\mathbf{v}\_m + (1 - a)^{-1} \left(1 + \frac{a}{1 - a} \frac{\rho\_d}{\rho\_c}\right)^{-1} \mathbf{v}\_{d\bar{\jmath}}\right) \cdot \mathbf{n} \tag{A40}$$

$$\mathbf{u} = a\rho\_d \left( \mathbf{v}\_m + (1 + a + \cdots) \left( 1 - \frac{a}{1 - a} \frac{\rho\_d}{\rho\_c} + \cdots \right) \mathbf{v}\_{d\bar{j}} \right) \cdot \mathbf{n} \tag{A41}$$

$$\mathbf{u} = a\rho\_d \left( \mathbf{v}\_m + \left( 1 + a\left( 1 - \frac{1}{1 - a} \frac{\rho\_d}{\rho\_c} \right) + \dotsb \right) \mathbf{v}\_{d\bar{\jmath}} \right) \cdot \mathbf{n}. \tag{A42}$$

As *α* 1 and *ρ<sup>d</sup>* ≈ 2*ρ<sup>c</sup>* =⇒ <sup>1</sup> <sup>−</sup> <sup>1</sup> 1−*α ρd ρc* <sup>&</sup>lt; 2 and ignoring terms containing squared and higher powers of *α*,

$$J\_{\Gamma} \approx \mathfrak{a} \rho\_d (\mathbf{v}\_m + \mathbf{v}\_{d\bar{j}}) \cdot \mathbf{n}.\tag{A43}$$

Applying the product rule,

$$\frac{\partial f\_{\Gamma}}{\partial \alpha} = \rho\_d (\mathbf{v}\_m + \mathbf{v}\_{dj}) \cdot \mathbf{n} + \alpha \rho\_d \frac{\partial \mathbf{v}\_{dj}}{\partial \alpha} \cdot \mathbf{n}$$

$$= \rho\_d \left(\mathbf{v}\_m + \mathbf{v}\_{dj} + \alpha \frac{\partial \mathbf{v}\_{dj}}{\partial \alpha}\right) \cdot \mathbf{n}.\tag{A44}$$

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
