Derivation of the Adjoint Drift Flux Equations for Multiphase Flow

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.

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,4,5], aerospace [6,7] and turbomachinery [8,9,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{array}{cc}\hfill \frac{\partial {\rho }_m{\mathbf{v}}_m}{\partial t}+({\mathbf{v}}_m·\nabla )\left({\rho }_m{\mathbf{v}}_m\right)& =-\nabla \left({\rho }_m{p}_m\right)+\nabla ·\left(2{\mu }_m\mathrm{D}\left({\mathbf{v}}_m\right)\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & -\nabla ·\left(\frac{\alpha }{1-\alpha }\frac{{\rho }_c{\rho }_d}{{\rho }_m}{\mathbf{v}}_{dj}{\mathbf{v}}_{dj}\right)+{\rho }_m\mathbf{g}+F,\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_m}{\partial t}+\nabla ·\left({\rho }_m{\mathbf{v}}_m\right)& =0,\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill \frac{\partial \alpha }{\partial t}+\nabla ·\left(\alpha {\mathbf{v}}_m\right)& =-\nabla ·\left(\frac{\alpha {\rho }_c}{{\rho }_m}{\mathbf{v}}_{dj}\right)+\nabla ·K\nabla \alpha ,\hfill \end{array}$$

(1c)

where:

• $\alpha$ is the dispersed-phase volume fraction,
• ${\rho }_c$ is the continuum density,
• ${\rho }_d$ is the dispersed-phase density,
• ${\rho }_m$ is the mixture density, defined as $\alpha {\rho }_d+(1-\alpha ){\rho }_c$,
• ${\mathbf{v}}_m$ is the mixture velocity,
• $p_m$ is the mixture kinematic pressure,
• ${\mu }_m$ is the mixture viscosity, defined as the sum of the continuum, dispersed-phase and mixture turbulent viscosities, ${\mu }_c+{\mu }_d+{\mu }_m^t$,
• $\mathrm{D}\left({\mathbf{v}}_m\right)=\frac{1}{2}\left(\nabla {\mathbf{v}}_m+{(\nabla {\mathbf{v}}_m)}^T\right)$ is the mixture strain rate tensor,
• ${\mathbf{v}}_{dj}$ is the dispersed-phase settling velocity,
• $\mathbf{g}$ is the acceleration due to gravity,
• F is the capillary force and
• K is the turbulent diffusion coefficient, defined as the mixture eddy diffusivity, ${\nu }_m^t=\frac{{\mu }_m^t}{{\rho }_m}$.

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.1. 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 𝓡, the optimisation problem can be stated as,

$$\mathrm{optimise}J(\mathrm{x},\mathrm{y})\mathrm{subject}\mathrm{to}𝓡(\mathrm{x},\mathrm{y})=0,$$

(2)

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

$$ℒ=J+{\int }_{\Omega }\lambda 𝓡d\Omega ,$$

(3)

where $ℒ$ is the Lagrange function, $\lambda$ are the Lagrange multipliers (also referred to as the adjoint variables) and $\Omega$ 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_1,R_2,R_3,R_4,R_5)}^T$, as follows:

$$\begin{array}{cc}\hfill {(R_1,R_2,R_3)}^T& =({\mathbf{v}}_m·\nabla )\left({\rho }_m{\mathbf{v}}_m\right)+\nabla \left({\rho }_m{p}_m\right)-\nabla ·\left(2{\mu }_m\mathrm{D}\left({\mathbf{v}}_m\right)\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +\nabla ·\left(\frac{\alpha }{1-\alpha }\frac{{\rho }_c{\rho }_d}{{\rho }_m}{\mathbf{v}}_{dj}{\mathbf{v}}_{dj}\right)-{\rho }_m\mathbf{g}+\aleph {\rho }_m{\mathbf{v}}_m,\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill R_4& =-\nabla ·\left({\rho }_m{\mathbf{v}}_m\right),\hfill \end{array}$$

(4b)

$$\begin{array}{cc}\hfill R_5& =\nabla ·\left(\alpha {\mathbf{v}}_m\right)+\nabla ·\left(\frac{\alpha {\rho }_c}{{\rho }_m}{\mathbf{v}}_{dj}\right)-\nabla ·K\nabla \alpha ,\hfill \end{array}$$

(4c)

where ℵ is the porosity, associated with the Darcy term. The variation of the Lagrange function with respect to the primal variables, $({\mathbf{v}}_m,p_m,\alpha )$, and the design parameter, ℵ, is,

$$\delta ℒ={\delta }_{{\mathbf{v}}_m}ℒ+{\delta }_{p_m}ℒ+{\delta }_{\alpha }ℒ+{\delta }_{\aleph }ℒ,$$

(5)

where, for example, ${\delta }_{\alpha }ℒ=ℒ(\alpha +\delta \alpha )-ℒ\left(\alpha \right)$. We choose the adjoint variables, $(\mathbf{u},q,\beta )=(u_1,u_2,u_3,q,\beta )$, so that the variation with respect to the primal variables vanishes, i.e.,

$${\delta }_{{\mathbf{v}}_m}ℒ+{\delta }_{p_m}ℒ+{\delta }_{\alpha }ℒ=0,$$

(6)

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

$$\delta ℒ={\delta }_{\aleph }ℒ={\delta }_{\aleph }J+{\int }_{\Omega }(\mathbf{u},q,\beta ){\delta }_{\aleph }𝓡d\Omega .$$

(7)

Derivation of the Adjoint Drift Flux Equations

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

$$\begin{array}{c}{\delta }_{{\mathbf{v}}_m}J+{\delta }_{p_m}J+{\delta }_{\alpha }J\hfill \\ \hfill +{\int }_{\Omega }(\mathbf{u},q,\beta ){\delta }_{{\mathbf{v}}_m}𝓡d\Omega +{\int }_{\Omega }(\mathbf{u},q,\beta ){\delta }_{p_m}𝓡d\Omega +{\int }_{\Omega }(\mathbf{u},q,\beta ){\delta }_{\alpha }𝓡d\Omega =0,\end{array}$$

(8)

which can be expanded to,

$$\begin{array}{c}{\delta }_{{\mathbf{v}}_m}J+{\delta }_{p_m}J+{\delta }_{\alpha }J+{\int }_{\Omega }d\Omega \mathbf{u}·{\delta }_{{\mathbf{v}}_m}{(R_1,R_2,R_3)}^T+{\int }_{\Omega }d\Omega q{\delta }_{{\mathbf{v}}_m}{R}_4+{\int }_{\Omega }d\Omega \beta {\delta }_{{\mathbf{v}}_m}{R}_5\hfill \\ +{\int }_{\Omega }d\Omega \mathbf{u}·{\delta }_{p_m}{(R_1,R_2,R_3)}^T+{\int }_{\Omega }d\Omega q{\delta }_{p_m}{R}_4+{\int }_{\Omega }d\Omega \beta {\delta }_{p_m}{R}_5\\ \hfill +{\int }_{\Omega }d\Omega \mathbf{u}·{\delta }_{\alpha }{(R_1,R_2,R_3)}^T+{\int }_{\Omega }d\Omega q{\delta }_{\alpha }{R}_4+{\int }_{\Omega }d\Omega \beta {\delta }_{\alpha }{R}_5=0.\end{array}$$

(9)

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

$$\begin{array}{cc}\hfill {\delta }_{{\mathbf{v}}_m}{(R_1,R_2,R_3)}^T& =(\delta {\mathbf{v}}_m·\nabla )\left({\rho }_m{\mathbf{v}}_m\right)+({\mathbf{v}}_m·\nabla )\left({\rho }_m\delta {\mathbf{v}}_m\right)-\nabla ·\left(2{\mu }_m\mathrm{D}\left(\delta {\mathbf{v}}_m\right)\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & -\nabla ·\left(2{\delta }_{{\mathbf{v}}_m}{\mu }_d\mathrm{D}\left({\mathbf{v}}_m\right)\right)+\aleph {\rho }_m\delta {\mathbf{v}}_m,\hfill \end{array}$$

(10a)

$$\begin{array}{cc}\hfill {\delta }_{{\mathbf{v}}_m}{R}_4& =-\nabla ·\left({\rho }_m\delta {\mathbf{v}}_m\right),\hfill \end{array}$$

(10b)

$$\begin{array}{cc}\hfill {\delta }_{{\mathbf{v}}_m}{R}_5& =\nabla ·\left(\alpha \delta {\mathbf{v}}_m\right),\hfill \end{array}$$

(10c)

$$\begin{array}{cc}\hfill {\delta }_{p_m}{(R_1,R_2,R_3)}^T& =\nabla \left({\rho }_m\delta p_m\right),\hfill \end{array}$$

(10d)

$$\begin{array}{cc}\hfill {\delta }_{p_m}{R}_4& =0,\hfill \end{array}$$

(10e)

$$\begin{array}{cc}\hfill {\delta }_{p_m}{R}_5& =0,\hfill \\ \hfill {\delta }_{\alpha }{(R_1,R_2,R_3)}^T& =({\rho }_d-{\rho }_c)\left(({\mathbf{v}}_m·\nabla )\left(\delta \alpha {\mathbf{v}}_m\right)+\nabla \left(\delta \alpha p_m\right)+\delta \alpha (\aleph {\mathbf{v}}_m-\mathbf{g})\right)\hfill \end{array}$$

(10f)

$$\begin{array}{cc}\hfill & -\nabla ·\left(2{\delta }_{\alpha }{\mu }_d\mathrm{D}\left({\mathbf{v}}_m\right)\right)+\nabla ·{\delta }_{\alpha }\left(\alpha {\rho }_d{\mathbf{v}}_{dj}{\mathbf{v}}_{dj}\right),\hfill \end{array}$$

(10g)

$$\begin{array}{cc}\hfill {\delta }_{\alpha }{R}_4& =-({\rho }_d-{\rho }_c)\nabla ·\left(\delta \alpha {\mathbf{v}}_m\right),\hfill \\ \hfill {\delta }_{\alpha }{R}_5& =\nabla ·\left(\delta \alpha {\mathbf{v}}_m\right)+\nabla ·{\delta }_{\alpha }\left(\alpha {\mathbf{v}}_{dj}\right)-\nabla ·\left(\frac{{\mu }_m^t}{{\rho }_c}\nabla \delta \alpha \right)\hfill \end{array}$$

(10h)

$$\begin{array}{cc}\hfill & +\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)\nabla ·(\delta \alpha \nabla \alpha )+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)\nabla ·(\alpha \nabla \delta \alpha ).\hfill \end{array}$$

(10i)

Derivation of Equations (10a), (10g) and (10i) can be found in Appendix A, Appendix B, Appendix C, respectively, where the variation of ${\mu }_m^t$ 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,

$$\begin{array}{c}{\delta }_{{\mathbf{v}}_m}J+{\delta }_{p_m}J+{\delta }_{\alpha }J\hfill \\ +{\int }_{\Omega }d\Omega \mathbf{u}·((\delta {\mathbf{v}}_m·\nabla )\left({\rho }_m{\mathbf{v}}_m\right)+({\mathbf{v}}_m·\nabla )\left({\rho }_m\delta {\mathbf{v}}_m\right)-\nabla ·\left(2{\mu }_m\mathrm{D}\left(\delta {\mathbf{v}}_m\right)\right)\\ -\nabla ·\left(2{\delta }_{{\mathbf{v}}_m}{\mu }_d\mathrm{D}\left({\mathbf{v}}_m\right)\right)+\aleph {\rho }_m\delta {\mathbf{v}}_m)-{\int }_{\Omega }d\Omega q\nabla ·\left({\rho }_m\delta {\mathbf{v}}_m\right)+{\int }_{\Omega }d\Omega \beta \nabla ·\left(\alpha \delta {\mathbf{v}}_m\right)\\ +{\int }_{\Omega }d\Omega \mathbf{u}·\nabla \left({\rho }_m\delta p_m\right)+({\rho }_d-{\rho }_c){\int }_{\Omega }d\Omega \mathbf{u}·\left(({\mathbf{v}}_m·\nabla )\left(\delta \alpha {\mathbf{v}}_m\right)+\nabla \left(\delta \alpha p_m\right)+\delta \alpha (\aleph {\mathbf{v}}_m-\mathbf{g})\right)\\ -{\int }_{\Omega }d\Omega \mathbf{u}·\left(\nabla ·\left(2{\delta }_{\alpha }{\mu }_d\mathrm{D}\left({\mathbf{v}}_m\right)\right)\right)+{\int }_{\Omega }d\Omega \mathbf{u}·\left(\nabla ·{\delta }_{\alpha }\left(\alpha {\rho }_d{\mathbf{v}}_{dj}{\mathbf{v}}_{dj}\right)\right)\\ -({\rho }_d-{\rho }_c){\int }_{\Omega }d\Omega q\nabla ·\left(\delta \alpha {\mathbf{v}}_m\right)+{\int }_{\Omega }d\Omega \beta (\nabla ·\left(\delta \alpha {\mathbf{v}}_m\right)+\nabla ·{\delta }_{\alpha }\left(\alpha {\mathbf{v}}_{dj}\right)-\nabla ·\left(\frac{{\mu }_m^t}{{\rho }_c}\nabla \delta \alpha \right)\\ \hfill +\left(\frac{{\rho }_d}{{\rho }_c}-1\right)\nabla ·\left(\frac{{\mu }_m^t}{{\rho }_c}\delta \alpha \nabla \alpha \right)+\left(\frac{{\rho }_d}{{\rho }_c}-1\right)\nabla ·\left(\frac{{\mu }_m^t}{{\rho }_c}\alpha \nabla \delta \alpha \right))=0.\end{array}$$

(11)

Decomposing the objective function into contributions from the boundary, $\Gamma$, and interior, $\Omega$, of the domain,

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

(12)

Equation (11) can be reformulated as,

$$\begin{array}{c}{\int }_{\Gamma }d\Gamma \left(\mathbf{n}(\mathbf{u}·{\rho }_m{\mathbf{v}}_m)+\mathbf{u}({\rho }_m{\mathbf{v}}_m·\mathbf{n})+2{\mu }_m\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)-q{\rho }_m\mathbf{n}+\alpha \beta \mathbf{n}+\frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_m}\right)·\delta {\mathbf{v}}_m\hfill \\ -{\int }_{\Gamma }d\Gamma \left(2{\mu }_m\mathbf{n}·\mathrm{D}\left(\delta {\mathbf{v}}_m\right)·\mathbf{u}+2{\delta }_{{\mathbf{v}}_m}{\mu }_d\mathbf{n}·\mathrm{D}\left({\mathbf{v}}_m\right)·\mathbf{u}-2{\delta }_{{\mathbf{v}}_m}{\mu }_d\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_m\right)\\ +{\int }_{\Omega }d\Omega (-\nabla \mathbf{u}·\left({\rho }_m{\mathbf{v}}_m\right)-({\rho }_m{\mathbf{v}}_m·\nabla )\mathbf{u}-\nabla ·\left(2{\mu }_m\mathrm{D}\left(\mathbf{u}\right)\right)\\ +\aleph {\rho }_m\mathbf{u}+{\rho }_m\nabla q-\alpha \nabla \beta +\frac{\partial J_{\Omega }}{\partial {\mathbf{v}}_m})·\delta {\mathbf{v}}_m-{\int }_{\Omega }d\Omega \nabla ·\left(2{\delta }_{{\mathbf{v}}_m}{\mu }_d\mathrm{D}\left(\mathbf{u}\right)\right)·{\mathbf{v}}_m\\ +{\int }_{\Gamma }d\Gamma \left({\rho }_m\mathbf{u}·\mathbf{n}+\frac{\partial J_{\Gamma }}{\partial p_m}\right)\delta p_m+{\int }_{\Omega }d\Omega \left(-\nabla ·{\rho }_m\mathbf{u}+\frac{\partial J_{\Omega }}{\partial p_m}\right)\delta p_m\\ +{\int }_{\Gamma }d\Gamma (({\rho }_d-{\rho }_c)\left(\mathbf{u}({\mathbf{v}}_m·\mathbf{n})·{\mathbf{v}}_m+p_m\mathbf{u}·\mathbf{n}-q{\mathbf{v}}_m·\mathbf{n}\right)+\beta {\mathbf{v}}_m·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}(\mathbf{n}·\nabla )\beta \\ +\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)\left(\beta (\mathbf{n}·\nabla )\alpha -\alpha (\mathbf{n}·\nabla )\beta \right)+\frac{\partial J_{\Gamma }}{\partial \alpha })\delta \alpha \\ +{\int }_{\Gamma }d\Gamma \left(\mathbf{u}·\mathbf{n}{\delta }_{\alpha }\left(\alpha {\rho }_d{\mathbf{v}}_{dj}{\mathbf{v}}_{dj}\right)+\beta {\delta }_{\alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}-2{\delta }_{\alpha }{\mu }_d\mathbf{n}·\mathrm{D}\left({\mathbf{v}}_m\right)·\mathbf{u}+2{\delta }_{\alpha }{\mu }_d\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_m\right)\\ -{\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}·\nabla )\delta \alpha \\ +{\int }_{\Omega }d\Omega (({\rho }_d-{\rho }_c)\left(-({\mathbf{v}}_m·\nabla )\mathbf{u}·{\mathbf{v}}_m-p_m\nabla ·\mathbf{u}+\mathbf{u}·(\aleph {\mathbf{v}}_m-\mathbf{g})+({\mathbf{v}}_m·\nabla )q\right)\\ -({\mathbf{v}}_m·\nabla )\beta +\frac{{\mu }_m^t}{{\rho }_c}\left(\left(\frac{{\rho }_d}{{\rho }_c}-1\right)\left(\nabla ·(\alpha \nabla \beta )-\nabla \alpha ·\nabla \beta \right)-\nabla ·\nabla \beta \right)+\frac{\partial J_{\Omega }}{\partial \alpha })\delta \alpha \\ \hfill -{\int }_{\Omega }d\Omega \left(\nabla ·\mathbf{u}{\delta }_{\alpha }\left(\alpha {\rho }_d{\mathbf{v}}_{dj}{\mathbf{v}}_{dj}\right)+{\delta }_{\alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\nabla \beta +\nabla ·\left(2{\delta }_{\alpha }{\mu }_d\mathrm{D}\left(\mathbf{u}\right)\right)·{\mathbf{v}}_m\right)=0.\end{array}$$

(13)

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{array}{cc}\hfill -\nabla \mathbf{u}·\left({\rho }_m{\mathbf{v}}_m\right)-({\rho }_m{\mathbf{v}}_m·\nabla )\mathbf{u}-\nabla ·\left(2{\mu }_m\mathrm{D}\left(\mathbf{u}\right)\right)& =-{\rho }_m\nabla q+\alpha \nabla \beta -\aleph {\rho }_m\mathbf{u}-\frac{\partial J_{\Omega }}{\partial {\mathbf{v}}_m}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +\nabla ·\left(2\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathrm{D}\left(\mathbf{u}\right)\right)·{\mathbf{v}}_m,\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill \nabla ·\left({\rho }_m\mathbf{u}\right)& =\frac{\partial J_{\Omega }}{\partial p_m},\hfill \\ \hfill -\left({\mathbf{v}}_m+\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)\nabla \alpha \right)·\nabla \beta & =\nabla ·\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)\nabla \beta \hfill \end{array}$$

(14b)

$$\begin{array}{cc}\hfill & +S_1+S_2-\frac{\partial J_{\Omega }}{\partial \alpha },\hfill \end{array}$$

(14c)

where:

$$\begin{array}{cc}\hfill S_1& =({\rho }_d-{\rho }_c)\left(({\mathbf{v}}_m·\nabla )(\mathbf{u}·{\mathbf{v}}_m-q)-\mathbf{u}·(\aleph {\mathbf{v}}_m-\mathbf{g})\right)+\nabla ·\left(2\frac{\partial {\mu }_d}{\partial \alpha }\mathrm{D}\left(\mathbf{u}\right)\right)·{\mathbf{v}}_m,\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill S_2& =\left(({\rho }_d-{\rho }_c)p_m+\frac{\partial }{\partial \alpha }\left(\alpha {\rho }_d{\mathbf{v}}_{dj}{\mathbf{v}}_{dj}\right)\right)\nabla ·\mathbf{u},\hfill \end{array}$$

(15b)

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

$$\begin{array}{c}\hfill {\int }_{\Gamma }d\Gamma (\mathbf{n}(\mathbf{u}·{\rho }_m{\mathbf{v}}_m)+{\rho }_m{\mathbf{v}}_m·\mathbf{n}\mathbf{u}+2{\mu }_m\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)-q{\rho }_m\mathbf{n}+\alpha \beta \mathbf{n}+\frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_m}\end{array}$$

$$\begin{array}{cc}\hfill +2\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathbf{n}·\left(\mathrm{D}\left({\mathbf{v}}_m\right)·\mathbf{u}-\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_m\right))·\delta {\mathbf{v}}_m-{\int }_{\Gamma }d\Gamma 2{\mu }_m\mathbf{n}·\mathrm{D}\left(\delta {\mathbf{v}}_m\right)·\mathbf{u}& =0,\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma \left({\rho }_m{u}_{\mathrm{n}}+\frac{\partial J_{\Gamma }}{\partial p_m}\right)\delta p_m& =0,\hfill \\ \hfill {\int }_{\Gamma }d\Gamma (\left({\mathbf{v}}_m·\mathbf{n}+\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)(\mathbf{n}·\nabla )\alpha \right)\beta \end{array}$$

(16b)

$$\begin{array}{c}\hfill +\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)(\mathbf{n}·\nabla )\beta +C_1+C_2+\frac{\partial J_{\Gamma }}{\partial \alpha })\delta \alpha \\ \hfill -{\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}·\nabla )\delta \alpha & =0,\hfill \end{array}$$

(16c)

where:

$$\begin{array}{cc}\hfill C_1& =\left(({\rho }_d-{\rho }_c)(\mathbf{u}·{\mathbf{v}}_m-q){\mathbf{v}}_m+2\frac{\partial {\mu }_d}{\partial \alpha }\left(\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_m\right)\right)·\mathbf{n},\hfill \end{array}$$

(17a)

$$\begin{array}{cc}\hfill C_2& =\left(\left(({\rho }_d-{\rho }_c)p_m+\frac{\partial }{\partial \alpha }\left(\alpha {\rho }_d{\mathbf{v}}_{dj}{\mathbf{v}}_{dj}\right)\right)\mathbf{u}-2\frac{\partial {\mu }_d}{\partial \alpha }\left(\mathrm{D}\left({\mathbf{v}}_m\right)·\mathbf{u}\right)\right)·\mathbf{n}\hfill \end{array}$$

(17b)

and $u_{\mathrm{n}}=\mathbf{u}·\mathbf{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{array}{cc}\hfill -\nabla \mathbf{u}·\left({\rho }_m{\mathbf{v}}_m\right)-({\rho }_m{\mathbf{v}}_m·\nabla )\mathbf{u}-\nabla ·\left(2{\mu }_m\mathrm{D}\left(\mathbf{u}\right)\right)& =-{\rho }_m\nabla q+\alpha \nabla \beta -\aleph {\rho }_m\mathbf{u}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +\nabla ·\left(2\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathrm{D}\left(\mathbf{u}\right)\right)·{\mathbf{v}}_m,\hfill \end{array}$$

(18a)

$$\begin{array}{cc}\hfill \nabla ·\left({\rho }_m\mathbf{u}\right)& =0,\hfill \\ \hfill -\left({\mathbf{v}}_m+\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)\nabla \alpha \right)·\nabla \beta & =\nabla ·\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)\nabla \beta \hfill \end{array}$$

(18b)

$$\begin{array}{cc}\hfill & +S_1.\hfill \end{array}$$

(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), $\nabla ·\mathbf{u}=0$ [22] and, therefore, $S_2=0$.

For the adjoint boundary conditions, the terms in Equation (16a) involving ${\mu }_m$ can be rewritten as,

$$\begin{array}{c}\hfill {\int }_{\Gamma }d\Gamma 2{\mu }_m\mathbf{n}·\left(\mathrm{D}\left(\mathbf{u}\right)·\delta {\mathbf{v}}_m-\mathrm{D}\left(\delta {\mathbf{v}}_m\right)·\mathbf{u}\right)={\int }_{\Gamma }d\Gamma {\mu }_m\left((\mathbf{n}·\nabla )\mathbf{u}·\delta {\mathbf{v}}_m-(\mathbf{n}·\nabla )\delta {\mathbf{v}}_m·\mathbf{u}\right)\end{array}$$

(19)

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

$$\begin{array}{c}\hfill {\int }_{\Gamma }d\Gamma (\mathbf{n}(\mathbf{u}·{\rho }_m{\mathbf{v}}_m)+{\rho }_m{\mathbf{v}}_m·\mathbf{n}\mathbf{u}+{\mu }_m(\mathbf{n}·\nabla )\mathbf{u}-q{\rho }_m\mathbf{n}+\alpha \beta \mathbf{n}+\frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_m}\end{array}$$

$$\begin{array}{cc}\hfill +2\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathbf{n}·\left(\mathrm{D}\left({\mathbf{v}}_m\right)·\mathbf{u}-\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_m\right))·\delta {\mathbf{v}}_m-{\int }_{\Gamma }d\Gamma {\mu }_m(\mathbf{n}·\nabla )\delta {\mathbf{v}}_m·\mathbf{u}& =0,\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma \left({\rho }_m{u}_{\mathrm{n}}+\frac{\partial J_{\Gamma }}{\partial p_m}\right)\delta p_m& =0,\hfill \\ \hfill {\int }_{\Gamma }d\Gamma (\left({\mathbf{v}}_m·\mathbf{n}+\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)(\mathbf{n}·\nabla )\alpha \right)\beta \\ \hfill +\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)(\mathbf{n}·\nabla )\beta +C_1+C_2+\frac{\partial J_{\Gamma }}{\partial \alpha })\delta \alpha \end{array}$$

(20b)

$$\begin{array}{cc}\hfill -{\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}·\nabla )\delta \alpha & =0.\hfill \end{array}$$

(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. Primal boundary conditions.

	${\mathbf{u}}_{\mathbf{m}}$	$\mathit{\alpha }$	${\mathit{p}}_{\mathit{m}}$
Inlet	fixed value	fixed value	zero gradient
Wall	zero	zero gradient	zero gradient
Outlet	zero gradient	zero gradient	zero

. We will derive expressions for the three main 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\mathrm{and}\delta \alpha =0.$$

(21)

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

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma {\mu }_m(\mathbf{n}·\nabla )\delta {\mathbf{v}}_m·\mathbf{u}& =0,\hfill \end{array}$$

(22a)

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma \left({\rho }_m{u}_{\mathrm{n}}+\frac{\partial J_{\Gamma }}{\partial p_m}\right)\delta p_m& =0,\hfill \end{array}$$

(22b)

$$\begin{array}{cc}\hfill {\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}·\nabla )\delta \alpha & =0.\hfill \end{array}$$

(22c)

When both fluids are incompressible, $\nabla ·{\mathbf{v}}_m=0$ [22], and as $\delta {\mathbf{v}}_{m\mathrm{t}}=0$ along the inlet, $(\mathbf{n}·\nabla )\delta {\mathbf{v}}_m=(\mathbf{n}·\nabla )\delta {\mathbf{v}}_{m\mathrm{t}}$ [19], where ${\mathbf{v}}_{m\mathrm{t}}$ is the tangential component of the mixture velocity. Hence, Equation (22) reduces to:

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma {\mu }_m(\mathbf{n}·\nabla )\delta {\mathbf{v}}_{m\mathrm{t}}·{\mathbf{u}}_{\mathrm{t}}& =0,\hfill \end{array}$$

(23a)

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma \left({\rho }_m{u}_{\mathrm{n}}+\frac{\partial J_{\Gamma }}{\partial p_m}\right)\delta p_m& =0,\hfill \end{array}$$

(23b)

$$\begin{array}{cc}\hfill {\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}·\nabla )\delta \alpha & =0,\hfill \end{array}$$

(23c)

where ${\mathbf{u}}_{\mathrm{t}}$ is the tangential component of the adjoint velocity, from which we deduce the boundary conditions for the adjoint variables at the inlet to be:

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{t}}& =0,\hfill \end{array}$$

(24a)

$$\begin{array}{cc}\hfill u_{\mathrm{n}}& =-\frac{1}{{\rho }_m}\frac{\partial J_{\Gamma }}{\partial p_m},\hfill \end{array}$$

(24b)

$$\begin{array}{cc}\hfill \beta & =0\iff {\mu }_m^t\ne 0.\hfill \end{array}$$

(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 $p_m$ enters the primal drift flux equations, the zero gradient boundary condition of $p_m$ at the inlet is applied to q as well,

$$(\mathbf{n}·\nabla )q=0.$$

(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,

$${\mathbf{v}}_m=0,\delta {\mathbf{v}}_m=0\mathrm{and}(\mathbf{n}·\nabla )\delta \alpha =0.$$

(26)

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 ${\mathbf{v}}_m$, go to zero. Equation (20) therefore reduces to:

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma {\mu }_m(\mathbf{n}·\nabla )\delta {\mathbf{v}}_m·\mathbf{u}& =0,\hfill \end{array}$$

(27a)

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma \left({\rho }_m{u}_{\mathrm{n}}+\frac{\partial J_{\Gamma }}{\partial p_m}\right)\delta p_m& =0,\hfill \\ \hfill {\int }_{\Gamma }d\Gamma (\left(\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)(\mathbf{n}·\nabla )\alpha \right)\beta \end{array}$$

(27b)

$$\begin{array}{cc}\hfill +\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)(\mathbf{n}·\nabla )\beta +C_2+\frac{\partial J_{\Gamma }}{\partial \alpha })\delta \alpha & =0.\hfill \end{array}$$

(27c)

As at the inlet, the primal velocity does not diverge and $\delta {\mathbf{v}}_{m\mathrm{t}}=0$ along the wall, so Equation (27) reduces to:

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma {\mu }_m(\mathbf{n}·\nabla )\delta {\mathbf{v}}_{m\mathrm{t}}·{\mathbf{u}}_{\mathrm{t}}& =0,\hfill \end{array}$$

(28a)

$$\begin{array}{cc}\hfill {\int }_{\Gamma }d\Gamma \left({\rho }_m{u}_{\mathrm{n}}+\frac{\partial J_{\Gamma }}{\partial p_m}\right)\delta p_m& =0,\hfill \\ \hfill {\int }_{\Gamma }d\Gamma (\left(\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)(\mathbf{n}·\nabla )\alpha \right)\beta \end{array}$$

(28b)

$$\begin{array}{cc}\hfill +\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)(\mathbf{n}·\nabla )\beta +C_2+\frac{\partial J_{\Gamma }}{\partial \alpha })\delta \alpha & =0,\hfill \end{array}$$

(28c)

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

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{t}}& =0,\hfill \end{array}$$

(29a)

$$\begin{array}{cc}\hfill u_{\mathrm{n}}& =-\frac{1}{{\rho }_m}\frac{\partial J_{\Gamma }}{\partial p_m},\hfill \\ \hfill \left(\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)(\mathbf{n}·\nabla )\alpha \right)\beta \end{array}$$

(29b)

$$\begin{array}{cc}\hfill +\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)(\mathbf{n}·\nabla )\beta & =-C_2-\frac{\partial J_{\Gamma }}{\partial \alpha }.\hfill \end{array}$$

(29c)

Equation (29c) is used to determine $\beta$ 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 p_m=0,(\mathbf{n}·\nabla )\delta {\mathbf{v}}_m=0\mathrm{and}(\mathbf{n}·\nabla )\delta \alpha =0.$$

(30)

The second integral in Equations (20a) and (20c) therefore goes to zero and, with $\delta p_m=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{array}{c}\hfill \mathbf{n}(\mathbf{u}·{\rho }_m{\mathbf{v}}_m)+{\rho }_m{\mathbf{v}}_m·\mathbf{n}\mathbf{u}+{\mu }_m(\mathbf{n}·\nabla )\mathbf{u}-q{\rho }_m\mathbf{n}+\alpha \beta \mathbf{n}+\frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_m}\end{array}$$

$$\begin{array}{cc}\hfill -2\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_m& =0,\hfill \\ \hfill \left({\mathbf{v}}_m·\mathbf{n}+\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)(\mathbf{n}·\nabla )\alpha \right)\beta \end{array}$$

(31a)

$$\begin{array}{cc}\hfill +\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)(\mathbf{n}·\nabla )\beta +C_1+C_2+\frac{\partial J_{\Gamma }}{\partial \alpha }& =0.\hfill \end{array}$$

(31b)

Note that the term containing $\mathrm{D}\left({\mathbf{v}}_m\right)=0$, because $(\mathbf{n}·\nabla ){\mathbf{v}}_m=0$. Decomposing Equation (31a) into its normal and tangential components yields:

$$\begin{array}{c}\hfill {\rho }_m\mathbf{u}·{\mathbf{v}}_m+{\rho }_m{u}_{\mathrm{n}}{\mathbf{v}}_m·\mathbf{n}+{\mu }_m(\mathbf{n}·\nabla )u_{\mathrm{n}}-{\rho }_{m}q+\alpha \beta +\frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_m·\mathbf{n}}\end{array}$$

$$\begin{array}{cc}\hfill -2\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_m·\mathbf{n}& =0,\hfill \end{array}$$

(32a)

$$\begin{array}{cc}\hfill {\rho }_m{\mathbf{v}}_m·\mathbf{n}{\mathbf{u}}_{\mathrm{t}}+{\mu }_m(\mathbf{n}·\nabla ){\mathbf{u}}_{\mathrm{t}}+\frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_{m\mathrm{t}}}-2\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_{m\mathrm{t}}& =0.\hfill \end{array}$$

(32b)

Equations (31b), (32a) and (32b) are used to determine $\beta$, q and ${\mathbf{u}}_{\mathrm{t}}$, respectively. Since $u_{\mathrm{n}}$ is prescribed at the inlet, the adjoint continuity equation, Equation (18b), is used to calculate $u_{\mathrm{n}}$ at the outlet, ${\Phi }_{\Gamma }$. The boundary conditions for the adjoint variables at the outlet are summarised as:

$$\begin{array}{cc}\hfill {\rho }_m{\mathbf{v}}_m·\mathbf{n}{\mathbf{u}}_{\mathrm{t}}+{\mu }_m(\mathbf{n}·\nabla ){\mathbf{u}}_{\mathrm{t}}-2\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_{m\mathrm{t}}& =-\frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_{m\mathrm{t}}},\hfill \end{array}$$

(33a)

$$\begin{array}{cc}\hfill u_{\mathrm{n}}& ={\Phi }_{\Gamma },\hfill \\ \hfill \left({\mathbf{v}}_m·\mathbf{n}+\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)(\mathbf{n}·\nabla )\alpha \right)\beta \end{array}$$

(33b)

$$\begin{array}{cc}\hfill +\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)(\mathbf{n}·\nabla )\beta & =-C_1-C_2-\frac{\partial J_{\Gamma }}{\partial \alpha },\hfill \\ \hfill \mathbf{u}·{\mathbf{v}}_m+u_{\mathrm{n}}{\mathbf{v}}_m·\mathbf{n}+{\nu }_m(\mathbf{n}·\nabla )u_{\mathrm{n}}+\frac{\alpha \beta }{{\rho }_m}+\frac{1}{{\rho }_m}\frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_m·\mathbf{n}}\end{array}$$

(33c)

$$\begin{array}{cc}\hfill -\frac{2}{{\rho }_m}\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_m·\mathbf{n}& =q,\hfill \end{array}$$

(33d)

where ${\nu }_m={\displaystyle \frac{{\mu }_m}{{\rho }_m}}$ is the mixture kinematic viscosity. A summary of the boundary conditions for the adjoint variables is presented in

Table 2. Adjoint boundary conditions for ducted flows.

	${\mathbf{u}}_{\mathbf{t}}$	u${}_{\mathbf{n}}$	$\mathit{\beta }$	q
Inlet	zero	Equation (24b)	zero	zero gradient
Wall	zero	Equation (29b)	Equation (29c)	zero gradient
Outlet	Equation (33a)	Equation (33b)	Equation (33c)	Equation (33d)

.

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·\mathbf{n},$$

(34)

where $\alpha {\rho }_d$ is the dispersed-phase mass fraction and ${\mathbf{v}}_d·\mathbf{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,$$

(35)

where ${}_o$ refers to the outlet. The derivatives of the objective function, Equation (34), with respect to the primal variables are:

$$\begin{array}{cc}\hfill \frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_m·\mathbf{n}}& =\alpha {\rho }_d,\hfill \end{array}$$

(36a)

$$\begin{array}{cc}\hfill \frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_{m\mathrm{t}}}& =0,\hfill \end{array}$$

(36b)

$$\begin{array}{cc}\hfill \frac{\partial J_{\Gamma }}{\partial \alpha }& ={\rho }_d\left({\mathbf{v}}_m+{\mathbf{v}}_{dj}+\alpha \frac{\partial {\mathbf{v}}_{dj}}{\partial \alpha }\right)·\mathbf{n},\hfill \end{array}$$

(36c)

$$\begin{array}{cc}\hfill \frac{\partial J_{\Gamma }}{\partial p_m}& =0.\hfill \end{array}$$

(36d)

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

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{t}}& =0,\hfill \end{array}$$

(37a)

$$\begin{array}{cc}\hfill u_{\mathrm{n}}& =0,\hfill \end{array}$$

(37b)

$$\begin{array}{cc}\hfill \beta & =0\iff {\mu }_m^t\ne 0,\hfill \end{array}$$

(37c)

$$\begin{array}{cc}\hfill (\mathbf{n}·\nabla )q& =0.\hfill \end{array}$$

(37d)

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

$$\begin{array}{cc}\hfill {\mathbf{u}}_{\mathrm{t}}& =0,\hfill \end{array}$$

(38a)

$$\begin{array}{cc}\hfill u_{\mathrm{n}}& =0,\hfill \\ \hfill \left(\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)(\mathbf{n}·\nabla )\alpha \right)\beta \end{array}$$

(38b)

$$\begin{array}{cc}\hfill +\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)(\mathbf{n}·\nabla )\beta & =0,\hfill \end{array}$$

(38c)

$$\begin{array}{cc}\hfill (\mathbf{n}·\nabla )q& =0.\hfill \end{array}$$

(38d)

Note that $C_2=0$ because $\mathbf{u}=0$. At an outlet, to satisfy the adjoint continuity equation, Equation (18b), $u_{\mathrm{n}}=0$, so:

$$\begin{array}{cc}\hfill {\mathbf{v}}_m·\mathbf{n}{\mathbf{u}}_{\mathrm{t}}+{\nu }_m(\mathbf{n}·\nabla ){\mathbf{u}}_{\mathrm{t}}-\frac{2}{{\rho }_m}\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_{m\mathrm{t}}& =0,\hfill \end{array}$$

(39a)

$$\begin{array}{cc}\hfill u_{\mathrm{n}}& =0,\hfill \\ \hfill \left({\mathbf{v}}_m·\mathbf{n}+\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)·\mathbf{n}+\frac{{\mu }_m^t}{{\rho }_c}\left(\frac{{\rho }_d}{{\rho }_c}-1\right)(\mathbf{n}·\nabla )\alpha \right)\beta \end{array}$$

(39b)

$$\begin{array}{cc}\hfill +\frac{{\mu }_m^t}{{\rho }_c}\left(1-\alpha \left(\frac{{\rho }_d}{{\rho }_c}-1\right)\right)(\mathbf{n}·\nabla )\beta & =-C_1-\frac{\partial J_{\Gamma }}{\partial \alpha },\hfill \end{array}$$

(39c)

$$\begin{array}{cc}\hfill \mathbf{u}·{\mathbf{v}}_m+{\nu }_m(\mathbf{n}·\nabla )u_{\mathrm{n}}+\frac{\alpha \beta }{{\rho }_m}+\frac{1}{{\rho }_m}\frac{\partial J_{\Gamma }}{\partial {\mathbf{v}}_m·\mathbf{n}}-\frac{2}{{\rho }_m}\frac{\partial {\mu }_d}{\partial {\mathbf{v}}_m}\mathbf{n}·\mathrm{D}\left(\mathbf{u}\right)·{\mathbf{v}}_m·\mathbf{n}& =q.\hfill \end{array}$$

(39d)

Note that $C_2=0$ because $u_{\mathrm{n}}=0$ and $\mathrm{D}\left({\mathbf{v}}_m\right)=0$ when $(\mathbf{n}·\nabla ){\mathbf{v}}_m=0$. A summary of the adjoint boundary conditions, using the objective function defined in Equation (34), is presented in

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

	${\mathbf{u}}_{\mathbf{t}}$	${\mathit{u}}_{\mathbf{n}}$	$\mathit{\beta }$	q
Inlet	zero	zero	zero	zero gradient
Wall	zero	zero	Equation (38c)	zero gradient
Outlet	Equation (39a)	zero	Equation (39c)	Equation (39c)

.

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 ${\mathbf{v}}_{dj}$ is modelled using,

$${\mathbf{v}}_{dj}={\mathbf{v}}_0{10}^{-k\alpha },$$

(40)

where ${\mathbf{v}}_0$ is the maximum theoretical settling velocity and k is a settling parameter, and its partial derivative with respect to $\alpha$ is given by,

$$\frac{\partial {\mathbf{v}}_{dj}}{\partial \alpha }=-kln10{\mathbf{v}}_{dj}.$$

(41)

5.2. Takacs Model

In this formulation ${\mathbf{v}}_{dj}$ is modelled using,

$$\begin{array}{c}{\mathbf{v}}_{dj}={\mathbf{v}}_0\left(e^{-a(\alpha -{\alpha }_r)}-e^{-a_1(\alpha -{\alpha }_r)}\right),\end{array}$$

(42a)

$$\begin{array}{c}\mathbf{0}⩽{\mathbf{v}}_{dj}⩽{\mathbf{v}}_{00},\end{array}$$

(42b)

where:

• a is the hindered settling parameter,
• $a_1$ is the flocculent settling parameter,
• ${\alpha }_r$ is the volume fraction of non-settleable solids at the inlet and
• ${\mathbf{v}}_{00}$ is the maximum practical settling velocity,

and its partial derivative with respect to $\alpha$ is given by,

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

(43)

For both models the partial derivative of $\alpha {\mathbf{v}}_{dj}$ with respect to $\alpha$ is given by,

$$\frac{\partial }{\partial \alpha }\left(\alpha {\mathbf{v}}_{dj}\right)={\mathbf{v}}_{dj}+\alpha \frac{\partial {\mathbf{v}}_{dj}}{\partial \alpha }.$$

(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.