Optimal control of diffuser shapes for non-uniform flow

A simplified model is used to identify the diffuser shape that maximises pressure recovery for several classes of non-uniform inflow. We find that optimal diffuser shapes strike a balance between not widening too soon, as this accentuates the non-uniform flow, and not staying narrow for too long, which is detrimental for wall drag. Three classes of non-uniform inflow are considered, with the axial velocity varying across the width of the diffuser entrance. The first case has inner and outer streams of different speeds, with a velocity jump between them that evolves into a shear layer downstream. The second case is a limiting case when these streams are of similar speed. The third case is a pure shear profile with linear velocity variation between the centre and outer edge of the diffuser. We describe the evolution of the time-averaged flow profile using a reduced mathematical model that has been previously tested against experiments and computational fluid dynamics models. The model consists of integrated mass and momentum equations, where wall drag is treated with a friction factor parameterisation. The governing equations of this model form the dynamics of an optimal control problem where the control is the diffuser channel shape. A numerical optimisation approach is used to solve the optimal control problem and Pontryagin’s maximum principle is used to find analytical solutions in the second and third cases. We show that some of the optimal diffuser shapes can be well approximated by piecewise linear sections. This suggests a low-dimensional parameterisation of the shapes, providing a structure in which more detailed and computationally expensive turbulence models can be used to find optimal shapes for more realistic flow behaviour.

pressure flow. Diffusers have numerous applications, from turbines in aerospace to hydropower [1][2][3] to automotive design [4]. There is a large body of literature on diffusers in the case where the inflow is uniform (see [5]), but only a limited literature available for non-uniform inlet flows [6].
In the case where the inflow is uniform, diffusers are usually designed to be straight sided, and the expansion angle is critical to performance [5]. The optimal angle strikes a balance between not being too shallow, since thin channels have larger wall drag, and not being too wide, since wide expansion angles result in boundary layer separation and poor consequent pressure recovery [7]. The optimum angle varies slightly, depending on the inflow boundary layer thickness, and whether the diffuser is two-dimensional or axisymmetric.
In the case where the inflow is non-uniform, there have been some experimental studies which have investigated diffuser performance for specific inflow profiles. Horlock and Lewis [8] studied an asymmetric shear flow and a symmetric wake flow in a linearly expanding two-dimensional diffuser using a simple inviscid model and experiments. They found that both of these non-uniform flow profiles became accentuated by the diffuser. The experiments of Wolf and Johnston [9] studied several non-uniform inlet flows for a two-dimensional diffuser. They considered an asymmetric flow profile with uniform shear, an asymmetric step-shear profile, a symmetric jet profile, and a symmetric wake profile. For each case, they studied straight sided diffusers with linearly expanding walls, using several different expansion angles and non-dimensional lengths. They found that these non-uniform inflows resulted in poorer diffuser performance than for a uniform inflow. They attributed this to the early onset of diffuser stall observed for these non-uniform inflows due to the accentuation of the flow profile. Other experimental studies are given by [10][11][12][13][14][15][16].
As is shown by Wolf and Johnston [9], the development of the non-uniform flow profile is strongly affected by the shape of the diffuser. An important feature in understanding non-uniform flows is the interplay between changes in the pressure and the kinetic energy flux factor, which is a normalised measure of how non-uniform a flow is [5]. A decrease in kinetic energy flux factor corresponds to a more uniform flow, and a rise in pressure. As shown by [8,9], diffusers with wide angles have the tendency to accentuate non-uniform flows and, in some extreme cases, create a jet-like outflow. In such cases, the outflow has a high-kinetic energy flux and, hence, a low-pressure recovery. On the other hand, diffusers with shallow angles have longer, narrower profiles, which create a lot of wall drag and consequently a larger drop in pressure. Optimal diffuser shapes, therefore, must strike a balance between mixing the flow in a narrow section and then widening the flow to decrease wall drag.
In this paper, we identify the optimal diffuser shape which satisfies these criteria. In contrast to diffusers with uniform inflow, where the channel shape is only restricted due to boundary layer separation, diffusers with nonuniform inflow have a shape which is also restricted due to the effect of accentuating the non-uniform flow. We find that in some cases, the optimal diffuser angle for non-uniform flow is smaller than that typically used for diffusers with uniform inflow. Furthermore, we show that, unlike for uniform flow, optimal diffuser shapes for non-uniform flow may contain an initial straight section that helps mix the flow before diffusing. Therefore, from a design perspective, the effect of the inflow profile cannot be ignored. In our analysis, we show how to optimise diffuser design based on the nature of the non-uniform inflow.
We investigate three different classes of non-uniform inflow, with the axial velocity varying across the width of the diffuser entrance. The first case has inner and outer streams of different speeds, with a velocity jump between them that evolves into a shear layer downstream, and the shear layer eventually interacts with the channel walls. The second case is the limit where the speeds of the streams are similar, creating a thin, slowly growing shear layer. In the third case, the inflow is a pure shear profile, with linear velocity variation between the centre and outer edge of the diffuser. These flow profiles are motivated by a low-head hydropower application, where the inner and outer streams are formed by a Venturi pipe which accelerates part of the flow in order to amplify the pressure drop across a turbine.
For these non-uniform flows which we consider, the development of the flow profile, which is fundamental to pressure recovery, can be described using a simple model for turbulent shear layers in confining channels [6]. The model assumes that the flow is composed of uniform streams separated by a linear shear layer. Wall drag is incorporated into the model with a friction factor, and the growth of the shear layers is modelled with a spreading parameter. As was shown by Benham et al. [6], the model predictions have good agreement with both CFD and experimental work for a range of channel shapes and Reynolds numbers. We restrict our attention to slender diffuser shapes, since our model, which is based on integrated equations of mass and momentum, applies to long and thin domains which are slowly varying. In this way, we avoid situations where there is boundary layer separation, which is not included in our model.
The effects of boundary layer separation have been described analytically in several recent studies [17][18][19][20][21][22]. Smith et al. [20] used a streamfunction formulation and asymptotic analysis to study the complex boundary layer structures that arise for a breakaway separation, and they calculated the corresponding separation point. Furthermore, for such separated flows, the analytical results of Kluwick and Scheichl [19] described how the velocity deficit between the outer flow and the flow in the boundary layer is related to the Reynolds shear stress. These studies are useful in determining the near-wall behaviour of the flow when there are strong adverse pressure gradients that cause separation. However, for the flows we consider in this study, we restrict our attention to diffuser shapes that are sufficiently slender that there is no risk of separation. Instead, we focus on the relationship between the diffuser shape and the development of the non-uniform bulk flow away from the walls.
The simple model which we use to describe the flow, was applied by Benham et al. [6] to investigate pressure recovery of a simple class of diffuser shapes by exhaustively searching a restricted design space. When we widen the parameter space and treat the diffuser shape as a continuous control, it is necessary to seek more complex tools to solve the problem. In this paper, we use the model as the basis for numerical optimisation of the diffuser shape, where the governing equations form the optimisation constraints. Such problems, as well as PDE-constrained optimisation problems, often arise in the field of flow control. With the advancement of computational power, these problems have become more feasible to solve. There are many different approaches to solving such problems which are discussed by Gunzburger [23].
In our approach, we exploit the fact that the model is one-dimensional and, upon discretisation, there are relatively few decision variables. This, in combination with the use of automatic differentiation to calculate gradients, allows us to use an interior point Newton method with relatively low computational effort [24]. In certain limiting cases, we solve the optimal control problem analytically using Pontragin's maximum principle [25] and these analytical results aid interpretation of the results from the numerical optimisation. We show that some of the optimal diffuser shapes look approximately like they are composed of piecewise linear sections. This motivates a low-dimensional parameterisation of the diffuser shapes, for which we use more detailed and computationally intensive CFD models to search for optima under more realistic flow behaviours. The two CFD models we use are a k-[26] and a k-ω Shear Stress Transport (SST) [27] turbulence model. We find that the optimal diffuser shapes for both these CFD models are very similar to those found using our reduced model. Section 2 outlines the model for the non-uniform flow profiles we consider and sets up the optimal control problem, discussing the choice of objective, the constraints, the number of parameters, and the optimisation approaches. In Sect. 3, we use a numerical optimisation approach to solve the optimal control problem and find optimal diffuser shapes in three different cases. In Sect. 4, we find analytical solutions to the optimal control problem in the last two of these cases. In Sect. 5, we present some CFD calculations and compare them to the results of the optimisation. Sect. 6 summarises the results of the paper and discusses the dependence of the optimal shapes on parameter choices.

Modelling turbulent shear layers in confining channels
In this section, we describe the flow scenarios which we consider and outline the simple model, previously presented by Benham et al. [6], which we use to describe these flows. This model is based on integrated conservation of mass and momentum equations in a long and thin geometry, as well as Bernoulli's equation, which govern an idealised time-averaged flow profile. A friction factor is used to parameterise the effect of wall drag, whilst a spreading parameter models the growth of shear layers. We restrict our attention to slender diffuser shapes with a small expansion angle, such that boundary layer separation is avoided.

(a) Developing shear layer
Centreline y = 0 y x (b) Small shear limit (c) Pure shear limit Fig. 1 Schematic diagram of the different flow cases. a Developing shear layer case, where the inflow has inner and outer streams of different speeds, with a velocity jump between them that develops into a shear layer. The channel length is sufficiently large such that the growing shear layer reaches across the channel. b Small shear limit case, which is the limiting case where the speeds of the streams are similar, such that the thin, slowly growing shear layer never reaches across the channel. c Pure shear limit, where the velocity varies linearly between the centre and the outer wall of the diffuser There are three different types of non-uniform channel flow we consider, all of which are symmetric about the channel centreline. In Fig. 1, we display each case, illustrating the axial velocity varying across the width of the diffuser. The first, which we call the Developing shear layer case, is an inflow composed of inner and outer streams of different speeds, with a velocity jump in between. In this case, a shear layer forms between the streams and grows downstream, eventually interacting with the channel walls. We restrict our attention to situations where the inner stream is slower than the outer stream. In other situations where the outer stream is slower than the inner stream, there is a greater risk of boundary layer separation, since the slowest region of flow is next to the wall [5]. Furthermore, it is well known that asymmetric flow instabilities, such as the Coanda effect [28], can occur in these situations, which we do not try to model here. The second case, called the Small shear limit, is similar to the first case, except the inner and outer streams have near-identical velocities, such that the shear between the flows is small and the thin shear layer grows slowly. In the third case, called the Pure shear limit, we consider a pure shear profile with linear velocity variation between the centre and outer edge of the diffuser. This corresponds to the first case in the downstream limit, where the shear layer has reached across the entire channel.
The simple model, presented by Benham et al. [6], is used to describe the idealised flow profiles for each of the three cases. The first two cases share the same formulation, whilst the third case is slightly different. Thus, we start by describing the governing equations for the first and second cases. Initially, we consider two-dimensional flow in a half channel 0 < y < h(x) and, later, we extend the model to axisymmetric channels (see Fig. 1a, where we indicate our coordinate axes). The inflow for the first two cases is composed of a slower moving central stream with speed U 2 and a faster outer stream with speed U 1 . A turbulent shear layer forms at the place where the parallel streams meet. We approximate the flow profile by decomposing it into two plug regions separated by a shear layer in which the velocity varies linearly between U 1 and U 2 (see Fig. 2). The approximate velocity profile is where h 1 and h 2 are the widths of the two plug regions, δ = h − h 1 − h 2 is the width of the shear layer, and ε y = (U 1 − U 2 )/δ is the shear rate. In the small shear limit, the plug flow speeds U 1 and U 2 are similar, such that the shear layer grows slowly. Whilst in the developing shear layer case, the shear layer may grow and interact with the channel walls, in the small shear limit, the channel is chosen to be sufficiently short that the slowly growing shear layer remains thin. However, in both cases, the shear rate decays with x as the shear layer grows [29]. We assume that the shear rate decays according to where S is a non-dimensional spreading parameter which must be determined from experiments or by comparison with CFD. Equation (2) can be derived from an entrainment argument (see Appendix in [6]), or by analogy with the growth of free shear layers. Assuming that the channel is long and thin, boundary layer theory [29] indicates that, to good approximation, the pressure does not vary across the channel width p = p(x). Averaged across the channel, conservation of mass and momentum equations are where ρ is the density, Q is the constant flow rate (per unit depth), and τ w is the wall shear stress. We parameterise the wall stress term with a friction factor f , such that τ w = −1/8 fρU 2 1 . We model the friction factor with the Blasius relationship, f = 0.316Re −1/4 , for flow in smooth pipes [30,31]. For all the examples we consider in this study, we use a finite value for the Reynolds number. For example, if Re = 10 6 , then f = 0.01. Although this friction factor is small, over the long diffuser length scales we consider, the effect of wall drag on pressure variations is significant.
Finally, we ignore viscous dissipation in the plug flow regions, since it is small compared to that at the walls and in the shear layer. Hence, in the plug regions, we assume Bernoulli's equation holds [32]. In certain cases, especially when the diffuser angle is wide, the speed of the slower plug region U 2 may decrease and reach zero. This has been observed in CFD simulations, which we display in Appendix A. In such cases, there is a portion of recirculating flow in the central part of the diffuser. We do not resolve the recirculation in these regions but since velocities are small, as observed in CFD, we treat the regions as stagnant zones with zero velocity (see Fig. 3).
Bernoulli's equation in each plug region holds along streamlines, ignoring transverse velocity components since they are small, and is implemented in a complementarity format for convenience The complementarity format of Eqs. (5) and (6) ensures that when either of the plug regions disappears, or if the slower plug region stagnates, Bernoulli's equation ceases to hold in that region. We find good comparison between our model predictions of the stagnant region and CFD calculations, which we discuss in Appendix A. To summarise the first and second cases, the simple model describes the evolution of the non-uniform velocity profile u(x, y), given by (1), and pressure p(x) in a symmetric confining channel. Equations (2)-(6) govern the variables U 1 , U 2 , h 1 , h 2 , δ, ε y and p, which are all functions of x. These equations can be solved for all x given inflow conditions at x = 0. Since the shear layer forms at x = 0, the inflow conditions for δ and ε y are δ(0) = 0 and ε y (0) = ∞. Pressure is measured with reference to the value at the inlet so we can take p(0) = 0 without loss of generality. All other inflow conditions form part of the set of parameters which we discuss in Sect. 2.2.
In the pure shear limit, the plug regions are non-existent, such that h 1 = h 2 = 0 and δ = h. Then the velocity profile takes the form In this case, the governing equations of the model reduce to (2)-(4) and (7), which govern the variables U 1 , U 2 , ε y and p.
We can extend the model to account for axisymmetric flows simply. For axisymmetric flow in a cylindrical channel 0 ≤ r ≤ h, we assume that that the velocity profile is identical to Eq. (1) in the first two cases, and (7) in the third case, except with y replaced by r . In the axisymmetric version of the model, Eqs.
The results of the axisymmetric and two-dimensional cases are compared in Sect. 3.

Formulation of the optimal control problem
In this section, we describe the optimal control problem by choosing an optimisation objective and formulating the control variables and contstraints. Starting with the objective, we note that diffuser performance can be measured in a number of different ways, for example using a pressure recovery coefficient or a loss coefficient [5]. The pressure recovery coefficient C p is a measure of the pressure gain in the diffuser from inlet to outlet, relative to the kinetic energy flux at the inlet. The loss coefficient K l is a measure of the total energy lost from inlet to outlet, relative to the kinetic energy flux at the inlet. For our optimal control problem, we could choose either of these coefficients as the objective. Maximising C p , for a given inflow, would produce the diffuser that converts the greatest amount of inflow kinetic energy into static pressure at the outflow. Minimising K l , for a given inflow, would produce the diffuser with the maximum amount of energy at the outflow. For this paper, we choose the pressure recovery coefficient as the objective. There are several ways to define the coefficient, but we shall use the so-called "mass-averaged" pressure recovery [11], which is defined as for the two-dimensional case and (11) for the axisymmetric case. The pressure recovery coefficient can take values C p ∈ (−∞, 1], where C p = 1 when all the kinetic energy of the inlet flow is converted into static pressure. For a given area ratio h(L)/ h(0) and inflow, there is a maximum possible pressure recovery C p I ≤ 1 [5]. For uniform inviscid flow this ideal limit is 4 in the axisymmetric case, but for non-uniform flow it is not known what the limit is. Now that we have chosen a suitable objective for the optimisation, we need to define a control. The diffuser shape is ultimately the control of the problem, but there are several different ways to formulate it. For example, we could use the shape function h(x) as the control, or we could use its derivative, or even the second derivative. To aid our choice of control, we consider the regularity requirements of the final shape. If the minimum requirement is that the shape be continuous, it will be convenient to choose the derivative of h as the control. If we also require smoothness (i.e. existence of the first derivative of h), then it will be convenient to choose the second derivative of h as the control. However, if no such requirements exist, then it is satisfactory to use h itself as the control. For this paper, we restrict ourselves to continuous but non-smooth shapes, and so we choose the shape derivative, or diffuser angle, as the control for optimising the diffuser shape. In reality, sudden expansions and sharp corners, if severe enough, can cause flow separation which is detrimental to pressure recovery [5]. Therefore, any such sharp corners must be rounded off with a suitable radius of curvature, upon construction. However, we neglect this concern from our mathematical analysis. An additional possible control of the problem is the channel length L. For now, we consider this fixed, but later we discuss the possibility of including L as a free parameter. After defining both the objective and the control of the optimisation, we now discuss the constraints. The most obvious constraints on the variables are the governing equations and inflow conditions. In addition, we may also want some constraints on the outflow. As mentioned earlier, constraining h(L)/ h(0) gives us a fixed maximum value for the pressure recovery. If h(L)/ h(0) is unconstrained, then the pressure recovery will be maximised with h(L)/ h(0) = ∞ [5]. However, this is impractical for construction and, due to Bernoulli's equation, we see that pressure recovery decays rapidly with h (like ∼ 1/ h 2 for two-dimensional flows and like ∼ 1/ h 4 for axisymmetric flows) so that a large majority of pressure is recovered for relatively small values of h(L)/ h(0). For example, if h(L)/ h(0) = 3 in uniform inviscid axisymmetric flow, the pressure recovery is C p I ≈ 0.99. Therefore, for practical considerations, we constrain the channel width at the outflow  List of the parameters of the optimal control problem. We treat the first 5 of these parameters as problem-specific, whereas the last 2 parameters are considered fixed. We also make use of the shorthand Another important constraint we need to consider is the boundedness of the control α. In particular, we note that for large values of the diffuser angle, boundary layers at the channel walls have the tendency to separate [7]. This phenomenon, which is often called 'diffuser stall', is not something that we attempt to capture with our model. However, it is known that diffuser stall has a detrimental effect on pressure recovery (because the flow does not slow down). Considering this, we give the control α an upper bound corresponding to the smallest diffuser angle which causes stall. The first appreciable stall of a straight walled diffuser is at α ≈ tan 7 • for the two-dimensional case and α ≈ tan 3.5 • for axisymmetric diffusers [5]. Furthermore, due to engineering constraints, it might not always be possible to construct channel shapes which contract more than a certain angle. Therefore, a lower bound on the control may also be necessary. If we denote the upper and lower bounds α max and α min , respectively, then α satisfies the box constraints It should be noted that, whilst Eq. (14) applies, the optimal control might not necessarily attain these bounding values. In such cases, Eq. (14) may be considered irrelevant.
To summarise the formulation of the optimal control problem, we seek to maximise the pressure recovery by manipulating the control α(x) within its bounds: with the constraints that Eqs. (2)-(6) hold, together with inlet conditions for all variables at x = 0, and the end constraint (13). Before moving on to the optimisation approaches, we note that there are several parameters which affect the solution. We list these parameters in Table 1 and discuss them in more detail in Sect. 6.

Optimisation approaches
Here, we describe the two main approaches we use to solve the optimisation problem. The first is a numerical approach, and the second is an analytical approach which uses Pontryagin's maximum principle. In each case, we give a brief description of the background theory, without going into depth. For a more thorough discussion, the reader is directed to the references in the text.

Numerical optimisation approach
For the numerical optimisation approach, we solve the optimisation problem (15) by discretising space, introducing values of the variables U 1 , U 2 , h 1 , h 2 , δ, ε y and p at each spatial point, and treating each discretised value as a degree of freedom. We use an interior point Newton method [24] (with the IpOpt library [33]) for nonlinear-constrained optimisation problems. Gradients are calculated using automatic differentiation in the JuMP package [34] of the Julia programming language [35]. The equality constraints we need to impose are Eqs. (2)-(6), inlet conditions and the terminal condition (13). There are also the inequality constraints (14) and those listed in the complementarity condition (5)- (6). As is often done, we impose the equality constraints using the quadratic penalty method [24], where we subtract their residual squared from the objective. For example, if our objective is to maximise the function g(x) subject to the equality constraint c(x) = 0, then we replace the objective with where μ is a penalty parameter. The interpretation of (16) is that by maximising g − μc 2 , we try to find the value of x that makes g as large as possible (to increase the objective value) and c 2 as small as possible (to impose the constraint equation). The value of μ must be chosen to be sufficiently large that the constraint equation is imposed accurately, but not too large that the problem becomes numerically ill-conditioned. Our inequality constraints are simple box constraints. These are dealt with by the interior point method using logarithmic barrier functions. More details, including how to choose the penalty parameter μ and barrier functions, are discussed by Nocedal and Wright [24]. It is not known whether this optimisation problem is convex so there may exist multiple solutions. In order to have confidence about the optimal solutions that are found, we use many different initial guesses to initialise the interior point method (although we have not yet found any multiple solutions).
The governing equations of the model consist of the algebraic equations, which are (3), (5)-(6) and inflow conditions, and the differential equations which are (2) and (4). We discretise space into n points and impose the algebraic equations exactly at every point. The differential equations are imposed using a second order backward finite difference scheme. It should be noted that whilst the complementarity conditions enforce a switch in the governing equations and may produce non-smooth behaviour in the solution, the equations themselves are smooth and can therefore be differentiated. Computation time is fast, owing to the use of automatic differentiation to calculate gradients (as opposed to finite differencing, for example). With 8 variables U 1 , U 2 , h 1 , h 2 , δ, ε y , p, h and one control α, the total number of degrees of freedom is 9n. For a discretisation of n = 100 grid points, and therefore 900 degrees of freedom, computation time is of the order of less than 10 s on a laptop computer.

Analytical optimisation approach
In the case of the small shear limit and the pure shear limit, we use an analytical optimisation approach. In this approach, we simplify the governing equations of our model (which are (2)-(6) for the small shear limit and (2)-(4), (7) for the pure shear limit), and we solve the simplified optimal control problem using Pontryagin's maximum principle [25,36].
Here, we briefly outline the mathematical details of Pontryagin's maximum principle which we will use later in Sects. 4.1 and 4.2. For the purpose of this outline, we consider a simplified optimal control problem in which a variable y(x) is controlled by a control α(x) and the domain is defined by 0 ≤ x ≤ L. We consider an objective functional for some given functions F and Φ, subject to the ordinary differential equation and initial/terminal conditions. We assume that the control α(x) is piecewise continuous, y(x) is absolutely continuous and F and G are of class C 2 in y and α, and piecewise continuous in x [36]. The Hamiltonian is defined as where λ(x) is the adjoint variable (or Lagrange multiplier). The set of equations which govern the state variable y, the adjoint variable λ and the control α, are known as Pontryagin's maximum principle, and these equations are The governing equations (20)- (22) are accompanied by the boundary conditions where y 0 and y L are initial and terminal conditions on the state variable y, where appropriate. If we treat the length of the domain L as a control variable, as well as α(x), then in addition to (20)- (24), there is the condition For both the small shear limit and the pure shear limit in Sects. 4.1 and 4.2, we find that F = 0 and both Φ and G are linear functions of the control α. As described by Pitcher [36] and McDanell and Powers [37], if in addition to the problem being linear in α, the control is also bounded above and below, as in (15), then (20) and (21) are exactly as before but (22) is replaced by where we have introduced the shorthand notation H α = ∂ H/∂α. When α(x) only takes its extreme values, then the control is said to be "bang-bang". If α(x) = α * (x) satisfies H α = 0 for a finite interval, then this is called a "singular arc". Robbins [38] showed that during the singular arc, α * is determined by solving where q is the smallest even number for which d q H α /dx q is not independent of α. Note that, if instead of a single-state variable y there are multiple state variables y i , as is the case in the examples in Sects. 4.1 and 4.2, then the analysis of this section still applies but with a Hamiltonian defined as Furthermore, (20), (21) and boundary conditions (23), (24) apply for each state variable and corresponding adjoint variable.

Numerical optimisation
This section, in which we solve the optimal control problem numerically using the approach outlined in Sect. 2.3.1 for several different cases, is divided into subsections for clarity. Firstly, we study the solution in the developing shear layer case (see Fig. 1a). In this case, we find that optimal diffuser shapes look approximately like they are composed of piecewise linear sections. This observation motivates us to introduce a low-dimensional parameterisation of the shapes that can be explored with contour plots, which is useful when comparing with CFD calculations later in Sect. 5. The small shear limit and the pure shear limit (see Fig. 1b, c) are discussed in Sects. 3.2 and 3.3, respectively.

The developing shear layer case
Having discussed the optimisation routine in Sect. 2.3.1, we use it to optimise channel shapes in several different cases, starting with the developing shear layer case. For plotting purposes, we maintain all variables in nondimensional form with reference to typical length scales and velocity scales. We use the initial channel half-width as a typical length scale h 0 = h(0) and the speed of the faster plug region at the inlet as a typical velocity scale In the developing shear layer case, we look at two-dimensional flow and choose parameter values The other parameters are taken as S = 0.11 (which we determine from comparison with CFD in Sect. 5) and f = 0.01 (which corresponds to a Reynolds number of Re = 10 6 and hydraulically smooth walls, using the Blasius relationship [30,31]). The number of grid points for discretising the simple model is n = 100. Simultaneously, we also investigate axisymmetric flow with the same parameter values, except with α max = tan 3.5 • . We plot the optimal diffuser angle for the twodimensional case in Fig. 4a, the corresponding optimal diffuser shape and velocity colour map in Fig. 4c, and pressure plot in Fig. 4g. The axisymmetric case is plotted in Fig. 4b, d, h.
In both cases, we observe that the optimal shape looks approximately like a piecewise linear function, which is divided into a straight part, followed by a widening part, followed by another straight part. In the two-dimensional case, the length of the first straight section aligns with the length it takes for the shear layer to spread completely across the channel. This suggests that mixing the flow to a more uniform profile is advantageous for the widening part to perform well. This is as expected because, as mentioned earlier, diffusers tend to accentuate non-uniform flow, producing an outflow with large kinetic energy flux (and therefore a low-pressure recovery). However, as discussed earlier, long thin channels cause large loss in pressure due to wall drag. Therefore, the optimal shape must have a straight section which is sufficiently long that the shear layer reaches across the channel, making the flow more well mixed, but no longer than that because of wall drag.
Interestingly, the widening part of the channel widens at a shallower angle than the maximum value (around tan 2.3 • compared to tan 7 • ). So the upper bound on α is not needed in this case. This behaviour is unexpected since diffusers are usually designed with a widening angle as close to tan 7 • as possible, regardless of the inflow. These results suggest that there is an optimal widening angle which is determined by the non-uniform inflow, rather than the risk of boundary layer separation.
Since we observe that the optimal shape in Fig. 4a-d looks approximately piecewise linear with three sections, we also try restricting the control α in this way to see if we can attain a near optimal solution with a piecewise linear shape. We parameterise α by splitting it into three parts: a straight part with α = 0 for 0 ≤ x < x 1 ; a widening part with constant α > 0 for x 1 ≤ x < x 2 , and a final straight part with α = 0 for x 2 ≤ x ≤ L. We treat x 1 and x 2 as control parameters and the value of α in the middle section is determined by the condition We optimise pressure recovery, using the same algorithm as before, but with α having only 2 degrees of freedom (DOF), the two parameters x 1 and x 2 , instead of 100 DOF. We plot the optimal diffuser angle in Fig. 4a, which is nearly identical to that obtained with 100 DOF. Moreover, the velocity colour map and pressure plot displayed in Fig. 4e, g both show a very close match. The pressure recovery coefficient for 2 DOF is C p = 0.5205, which is the same as for 100 DOF (up to 4 decimal places), suggesting that piecewise linear diffuser shapes are a very good approximation in this case. In Fig. 5a, we display a contour plot of C p for all possible values of x 1 and x 2 (we cut out part of the contour plot corresponding to α > tan 7 • ). This indicates that there is a clear unique optimum at x 1 / h 0 = 12.3 and x 2 / h 0 = 24.3. Note that an unoptimised shape, say with x 1 = 0 and x 2 = 5, gives a value of C p = 0.4254, which is 22% worse than the optimal shape. For the axisymmetric case, we find that the optimal shape has a similar structure and can also be well approximated by parameterisation with x 1 and x 2 . We find that the optimal value of x 1 / h 0 = 6.3 is a little bit shorter than the two-dimensional case. In fact, we can see in channel by x = x 1 . Instead, x 1 corresponds to the point where the shear layer reaches the centre of the channel. It reaches the outer wall of the channel slightly further downstream, during the widening section. This could be explained by the fact that pressure gradients due to wall drag are stronger (per unit flux) in the axisymmetric case (e.g. Poiseuille flow [32]), and therefore the optimal shape cannot afford a longer section of straight, narrow channel. Figure 4b, d, f, h shows a close comparison between the diffuser angle, velocity and pressure in the 2 DOF case and the 100 DOF case. The contour plot in Fig. 5b shows the optimal parameter values x 1 / h 0 = 6.3 and x 2 / h 0 = 19.8. The pressure recovery coefficient for 2 DOF is C p = 0.6886 compared to C p = 0.7018 for 100 DOF, suggesting that axisymmetric piecewise linear diffuser shapes are also a good approximation, but slightly less so than in the two-dimensional case. The contour plot in Fig. 5b has a steep gradient for small x 1 / h 0 , indicating that diffuser performance is poor in this case. This corresponds to situations in which the inner stream stagnates, and we model this with a zone which has zero velocity, as described earlier. We indicate all diffuser shapes where we do this by plotting a black line on top of the contours (there is no stagnation for the two dimensional case in Fig. 5a). The phenomenon of stagnation is particularly an issue when the inner stream is slow, and this analysis shows that the way to avoid such poor performance is to have a longer straight section in which the inner stream is accelerated before diffusing.

Small shear limit
For the small shear limit, we consider two-dimensional flow and choose parameter values U 2 (0)/U 0 = 0.8, h 2 (0)/ h 0 = 0.5, h(L)/ h 0 = 2.3, L/ h 0 = 20, α min = 0 and α max = tan 7 • . The other parameters S and f are taken at the same values as the previous cases. We display the optimal shape, velocity colour map and pressure plot in Fig. 6a, c, e. In this case, the optimal shape widens at the maximum angle α max until it reaches the exit width h(L) and then stays straight. Since the shear is small and the flow is almost uniform, there is no risk of accentuating the flow profile drastically. Therefore, wide angles initially are not penalised much, allowing the control to take its maximum value. This design is typically what is built in the diffuser industry for uniform inflow [5], where the maximum diffuser angle is set by the limit where boundary layer separation occurs. The optimal control only takes its extremal values, which is sometimes referred to as bang-bang control (see Sect. 2.3.2). In Sect. 4.1, we discuss analytical results for this limiting case and prove that the control must be bang-bang. The optimal channel, velocity colour map, pressure and control are displayed in Fig. 6b, d, f. The optimal shape is similar to those in Fig. 4, with a natural decomposition into two straight sections separated by a widening section. The widening section has an angle that increases from α ≈ tan 2.5 • to α ≈ tan 3 • and is nowhere greater than tan 7 • , showing that the upper bound α was not needed in this case. The optimal shape, like in Fig. 4, exhibits the balance between the necessity of a straight section that is long enough to allow some mixing, but not too long that wall drag dominates. In this case, which does not involve any of the switching behaviour that occurs when plug regions reach the wall, we derive some analytical results (discussed in Sect. 4.2) which support and help interpret the numerical optimisation. In particular, we investigate the nature of the optimal widening angle which lies in the interval (α min , α max ). This is of great interest because it indicates that the optimal design is unaffected by the conventional widening angle limit that exists due to boundary layer separation.

Analytical results
The numerical optimisation routine outlined in Sect. 2.3.1 can be applied to find optimal shapes for any choice of the parameters listed in Table 1. We have seen several examples of these in Figs. 4 and 6. In this section, we show that in the two limiting cases displayed in Fig. 6, the small shear limit and the pure shear limit, it is possible to make some analytical progress which aids our understanding and interpretation of the optimal control. Furthermore, the results discussed in this section include simple relationships that may be of instructive use for the purpose of diffuser design in industry. In both cases, we derive a reduced set of equations describing the dynamics, that is amenable to optimal control analysis using Pontryagin's maximum principle [25].

Small shear limit
To start with, we consider a two-dimensional diffuser where the inflow, given by (1), is almost uniform. We study this situation because it illustrates the limiting case of the earlier developing shear layer examples in Fig. 4, except where the inner and outer plug regions are of similar speeds. In such situations, the shear rate ε y is small, such that the shear layer develops slowly between the plug flow regions (see Eq. (2)). Furthermore, we restrict our attention to situations where the channel is sufficiently short that the shear layer never reaches across the channel (see Fig. 1b). Therefore, we consider an inflow given by (1) with small U 1 −U 2 , which is a perturbation to a uniform velocity profile. Then, we asymptotically expand our governing equations (1)-(6) about the leading order uniform flow solution, resulting in a reduced set of equations which are amenable to analysis using Pontryagin's maximum principle.
Hence, let us introduce the small parameter 1, which is defined by the difference in speed between the two plug regions at the inflow where V is a dimensional velocity scale. If the plug regions always exist, with positive width h 1 > 0, h 2 > 0, and we assume that the slower flow never stagnates U 2 > 0, then we need not consider the complementarity format for Bernoulli's equation. Therefore, we replace Eqs. (5) and (6) with We consider the distinguished limit where the friction factor f is small such that f = S F, where F = O (1).
Note that, whilst f is small, over the long diffuser length scales we consider, friction has a significant effect on the pressure variations. For a Reynolds number of Re = 10 6 and hydraulically smooth walls, the friction factor is calculated as f = 0.01, using the Blasius relationship [30,31]. Therefore, if = 0.1 and S = 0.11, then F = 0.91. Choosing these parameter values and setting V /U 0 = 2, we achieve the small shear limit example in Fig. 6a, c, e. By considering perturbations to uniform flow, we expand variables in powers of the small parameter , In the limit → 0, Eqs. (2)-(4) and (31) are satisfied by (43) The function h 1 0 , which represents the location of the centre of the shear layer to leading order, is determined at order O( ). Bernoulli's equation (31) for each plug region, at order O( ), iŝ From the relationship h 1 + h 2 + δ = h at order O( ), we find that Thus, the conservation of mass equation (3) iŝ and the momentum equation (4) is Thus, using Eqs. (44)-(46), we can simplify Eq. (48) to an equation purely involving h 1 0 and h

Equation (49) has solution
Combining (47) and (48), we obtain a differential equation for the pressure correctionp in terms of the channel shape and its derivative α = h (x), We now solve the optimal control problem outlined in Sect. 2.2, using the approach outlined in Sect. 2.3.2. Since the inflow conditions are fixed and we take p(0) = 0, maximising C p is equivalent to maximising pressure at the outlet p(L). Furthermore, the constraint equations have been reduced to (51). Therefore the optimal control problem, including terms up to and including order O( ), and written as a system of first order differential equations, is as follows: such that We now solve this reduced problem using Pontryagin's maximum principle [25], as outlined earlier. The Hamiltonian for this system is where λ p and λ h are the adjoint variables which satisfy the adjoint equations From (24), since the objective function only depends on pressure at the outlet Φ = p(L), we have the natural boundary condition Considering Eq. (61) and the fact that there is no dependance of the Hamiltonian (58) on the pressure p, Eq. (59) tells us that λ p = 1 for all values of x. There is no natural boundary condition for λ h since we are enforcing a condition on h at the outlet x = L (see (24)). The last condition from Pontryagin's maximum principle is the optimality condition, and since the Hamiltonian is linear in the control α, this takes the form of Eq. (26). Next we investigate whether the optimal control is bang-bang, or whether any singular arcs exist (see Sect. 2.3.2). Using (58) with (53), (54) and (60), and noting that λ p = 1, we see that which is negative for all values of x. Hence, it is impossible for singular arcs to exist in this case. Therefore the control is bang-bang with where the switching point γ is given by In Fig. 6c, e, we plot the solution to the optimal control problem found using the Hamiltonian approach on top of the solution found using the numerical optimisation routine outlined in Sect. 2.3.1. It is clear that the numerical optimisation routine has correctly found the bang-bang control which we have derived here, with γ / h 0 = 10.59 (the small discrepancy is probably due to the finite value of ).
It should be noted that the adjoint variable λ h is only solved for up to a constant of integration C (from integrating Eq. (60)) since it has no boundary condition. Instead, C is determined by the condition that In the case where we also allow the channel length L to be a control as well as α, as described in Sect. 2.3.2, we have the additional constraint on the Hamiltonian at the final point (25). By calculating H (58), it is straightforward to show that (65) and (25) are inconsistent unless γ = 0 or γ = L. For the case in Fig. 6a, c, e, it is clear that γ = 0 is impossible, so we conclude that the optimal diffuser length is L = γ . Therefore, including L as a control and taking α min = 0 • , the optimal diffuser shape for the small shear limit is one which expands at the maximum angle until h reaches h L , at which point the channel terminates (i.e. conventional diffuser design for uniform flow).

Pure shear limit
The next limiting case we investigate is the pure shear limit, in which the shear layer has already reached across the channel at the inflow, such that there are no plug regions (see Fig. 1c). The velocity profile is given by Eq. (7). For this velocity profile, conservation of mass and momentum equations (3), (4) reduce to We now solve the optimal control problem outlined in Sect. 2.2. As in Sect. 4.1, we maximise pressure at the outlet p(L). Furthermore, it is convenient to introduce a new variable, the scaled velocity difference between maximum and minimum velocitiesŨ = (U 1 − U 2 )/(U 1 + U 2 ). Using (66) and this new variable we can simplify Eqs. (2) and (67), which are the reduced system of constraint equations. Therefore, the optimal control problem is as follows: such that . Similarly to Sect. 4.1, the Hamiltonian for the system is constructed as which is linear in the control α. The adjoint equations are which have the natural boundary conditions There is no natural boundary condition for λ h since h is already prescribed at x = L (see (24)). Finally, as in Sect. 4.1, the optimality condition is (26). As described in Sect. 2.3.2, in order for there to be a singular arc, we must have H α = 0 and d H α /dx = 0 for a finite interval, where during this interval the value of the control is given by (27). In this case, we find that q = 2, such that (27) becomes Setting H α = dH α /dx = 0, and using Eq. (82), we find that there will only be a singular arc when for some x 2 > x 1 . We can solve the coupled system (26), (69)-(83) numerically for the optimal control and corresponding solution. There are, however, very special cases where the singular arc value α * is constant, when we can find analytical solutions, which we discuss at the end of this section. In practice, these solutions appear to be close to the behaviour observed in normal conditions, as can be seen in Fig. 6f where α is approximately constant during the widening section of the diffuser. In Fig. 6b, d, f, we plot the solution to the optimal control problem found using the Hamiltonian approach over the solution found using the numerical optimisation routine outlined in Sect. 2.3.1. It is clear that both approaches have found the same solution, with the singular arc lying between x 1 / h 0 = 2.9 and x 2 / h 0 = 32.9. The singular arc represents the balance between mixing and widening effects in the diffuser. Mixing reduces the non-uniformity of the flow, whilst widening tends to accentuate the non-uniform profile, and such non-uniformity can produce a high-kinetic energy, low-pressure outlet. Therefore, the control initially takes its minimum value α = 0 • to create good mixing. However, a straight section which is too long is detrimental to pressure recovery because of wall drag. Hence, a widening section is required after a critical length.
The optimal value of α in the widening section represents a balance between mixing and widening the flow. If α is too large, then the flow profile will become too non-uniform at the outlet. Conversely, if α is too shallow, wall drag losses are enhanced. The singular arc is interesting from both a mathematical point of view, but also from an engineering point of view. It clearly shows that diffuser designs for non-uniform inflow should take into account the nature of the non-uniform inflow profile. The optimal widening angle for manufacture is given by Eq. (83). Calculating α * (x) from (83) is difficult in general (Ũ is a variable), but we find that for certain parameter values, it takes a simpler and more useful form.
In general, during the singular arc α * (x) is not constant, yet in certain cases, such as Fig. 5b, it does not vary much over the singular arc interval. This raises the question of whether it is possible to find constant α * solutions. Noticing how the only variable in Eq. (83) isŨ , we seek solutions with constantŨ . Furthermore, we restrict our attention to solutions which begin on the singular arc. From Eq. (69), we see that constant α * solutions only exist if Therefore, reconciling Eqs. (83) and (84), it can be shown that constant singular arc solutions exist for parameters which satisfy ExcludingŨ = 0, and assuming that f and S are fixed, we are left with solving Eq. (85) forŨ . Substituting U 2 (0) and U 0 back into (85), we rewrite the equation in terms of the inlet velocity ratio U (0) = U 2 (0)/U 0 , so that we require the parameters to satisfy Similarly, (84) then indicates that For U (0) = 0.35, as in Fig. 6b, d, f, we find α * = tan 2.76 • , which is very close to the solution from the numerical optimisation, indicating that (87) is also a good approximation for non-constant singular arcs.
Considering that in typical situations f /S ≈ 0.1, and this is small, we expand (86) and (87) about f /S and ignore imaginary solutions, giving the approximate solution For the constant singular arc solution (87), we can solve the system (26), (69)-(83) analytically, though we do not include the details here. Note, we also need to ensure that both H α = 0 and d H α /dx = 0 for all values of x during the singular arc. This will enforce a further constraint on the other parameters of the problem. We do not include the details of this here, but instead, we simply state the constraint, written in terms of the rescaled non-dimensional channel length L/(Sh 0 ), and both h L / h 0 and f /S. Written in an asymptotic expansion in powers of f /S (similarly to before) this constraint takes the approximate form As an example of such a constant singular arc solution, we choose parameter values S = 0.11, f = 0.01 and α min = 0 • , from which, solving Eqs. (86) and (87), we have an inflow velocity ratio U (0) = 0.471 and a singular arc value of α * = tan 2.26 • . Using (90), we find that H α = 0 and dH α /dx = 0 along the singular arc if we choose the remaining parameter values h L / h 0 = 2 and L/ h 0 = 40. In Fig. 7, we display the velocity colour map for this solution, together with a plot of H α . We see that the solution starts on the singular arc until the expansion ratio h L / h 0 = 2 is reached, at which point H α becomes negative, such that the remaining length of the diffuser has angle α min = 0 • .
It is of further interest to investigate how the singular arc depends on the model parameters f and S. We plot the relationship between the constant singular arc value α * (87) and these parameters in Fig. 8. It is clear that increasing the friction factor f results in a higher α * . This is to be expected, since larger wall drag will penalise smaller angles more. Increasing the spreading parameter S also increases α * . This is because higher spreading rates results in better mixing, and hence, wider angles are more affordable.

Comparison with results from a k-model and a k-ω SST model
In this section, we discuss comparisons between the optimal shapes found using the simple model in the previous sections to calculations from CFD. Benham et al. [6] make comparisons between this model and a k-turbulence model [26], as well as experimental data generated with Particle Image Velocimetry (PIV). Here, we use both a k-and a k-ω Shear Stress Transport (SST) model [27] to compare with some of the simple model optimisation results. The k-model is one of the most popular computational turbulence models, whilst the k-ω SST model is Comparison between mathematical model and two computational turbulence models (k-and k-ω SST) for the optimal shape in Fig. 4c. a Velocity colour map calculated using the reduced model, with black dashed lines indicating the shear layer. b Corresponding velocity colour map calculated using the k-model. c Pressure, averaged across the channel width. d Velocity profiles at evenly spaced locations in the channel particularly robust in situations with strong adverse pressure gradients [27] (though for the small diffuser angles we consider in this study the adverse pressure gradients are not severe). Moreover, since a thorough comparison between the mathematical model and CFD has already been discussed in [6], we do not perform comparisons for all of the optimisation results of Sect. 3. Instead, we look at the geometry in Fig. 4e as a single example. Consider the example in Fig. 4e, as discussed in Sect. 3. In order to compare with the mathematical model, we use precisely the same inlet velocity profile in the CFD. Inlet conditions for the turbulence variables k, and ω are given by the free-stream boundary conditions [29] k = I 2 × 3/2 u 2 + v 2 , = 0.09k 3/2 / and ω = √ k/ , with turbulence intensity I = 10% and mixing length = 0.1h 0 (10% of the channel half-width). In both the CFD models, no slip boundary conditions are applied to the channel walls. Furthermore, we use all the standard turbulence parameter values, which are given by Launder and Spalding [26] for the k-model and Menter [27] for the k-ω SST model.
In Fig. 9a, b, we display colour plots of the time-averaged streamwise velocity u generated with both the simple model and the k-model. the simple model and both CFD models. There is good agreement between the models, with the simple model capturing the dominant features of the flow, such as maximum and minimum velocities, and the width of the shear layer. There is a slight discrepancy near the diffuser wall since our model does not resolve boundary layers, but instead parameterises their effect with a friction factor. However, we can see that our model accurately captures the effect of the boundary layers on the pressure by the close comparison between the models in Fig. 9c.
In Sect. 3, we investigated reducing the dimension of the control α by splitting it into three piecewise constant sections divided by x 1 and x 2 , which we treated as free parameters. Motivated by these piecewise linear shapes, here we make the same simplification, reducing the degrees of freedom of the control to 2. We explore the parameter space generated by x 1 and x 2 , using both CFD models to calculate pressure recovery. Each calculation made by the CFD models is much more computationally expensive than that of the simple model, but because of the low dimension of the degrees of freedom, we can still feasibly explore the different possible combinations of x 1 and x 2 . This would not be tractable, however, if we were to use 100 degrees of freedom, as we did with the simple model in Sect. 3. Hence, the numerical optimisation using the simple model is very useful for finding the general shape of the optimal channels, around which we can further search for optima using more realistic, yet more computationally intensive CFD models.
In Fig. 10, we plot contours of pressure recovery C p , given by (10), (11), as a function of the two parameters x 1 and x 2 , where C p is calculated using the k-model instead of the simplified model, as in Fig. 5. Similarly to the contour plots in Fig. 5, we exclude values of x 1 and x 2 which result in a diffuser angle larger than tan 7 • for the two-dimensional case and tan 3.5 • for the axisymmetric case.
Comparing Figs. 5 and 10, we see that the optimal diffuser shape using both the simple model and the kmodel, is characterised by similar values of x 1 and x 2 . In the two-dimensional case, according to the k-model, the optimal diffuser shape has x 1 / h 0 = 13 and x 2 / h 0 = 25, with a pressure recovery of C p = 0.5370. According to the simplified model, the optimal diffuser shape has x 1 / h 0 = 12.3 and x 2 / h 0 = 24.3, with a pressure recovery of C p = 0.5205, which is very close to that obtained with the k-model. Similarly, for the axisymmetric case, the k-model suggests an optimal diffuser shape with x 1 / h 0 = 7 and x 2 / h 0 = 21, giving a pressure recovery of C p = 0.7067, whereas the simplified model suggests x 1 / h 0 = 6.3 and x 2 / h 0 = 19.8, with a pressure recovery of C p = 0.6886. Considering that the diffuser angle is given by Eq. (29), it is clear that if the optimal values of x 1 and x 2 are similar, according to the simple model and the CFD, then the optimal diffuser angle is also similar. In fact, we can compare the value of C p as a function of diffuser angle by looking at intersections of the contour plots (Figs. 5, 10) with the lines x 1 − x 2 = const. We have also generated these pressure recovery data using the k − ω SST model, and we find the results very similar. The average discrepancy between the C p values calculated using the k-ω SST model and the k-model is 0.004 for the two-dimensional case, and 0.003 for the axisymmetric case.
These results indicate that the optimal shapes found using the numerical optimisation routine and the simplified model are similar to the optimal shapes that would be found if we were to use either of these computational turbulence models as a forward model. Hence, this gives us confidence that the optimal shapes generated using the simplified model are close to true optimal shapes in reality.

Discussion and conclusion
6.1 The effect of parameter values on the optimal shape Although we have investigated optimal diffuser shapes in a number of specific cases, we have not yet explored the various parameters of the model thoroughly, which are listed in Table 1. We now briefly discuss the effect that each of these parameters has on the optimal shapes. However, since there are many parameters, we do not provide plots for the analysis of every single parameter.
One of the most important parameters is the velocity ratio U 2 (0)/U 0 of the inflow. To explore this parameter, we investigate optimal diffuser shapes for a fixed inflow with h 2 (0)/ h 0 = 0.5, an expansion ratio h L / h 0 = 2.3 and a length ratio L/ h 0 = 40, and we vary the velocity ratio. The results of the optimisation are displayed in Fig. 11, for U 2 (0)/U 0 = 0.3 and U 2 (0)/U 0 = 0.7. We see that the effect of a velocity ratio which is closer to 1 is that an initial widening section becomes favourable. This is because when the inflow is more uniform, wider angles penalise pressure recovery less. Therefore, the balance is tipped in favour of reducing wall drag by expanding the channel a little. In the extreme case where U 2 (0)/U 0 becomes close to unity, we have seen in Sect. 4.1 that this initial widening section dominates throughout, such that the optimal control is purely bang-bang, with no singular arc. Notice how in the case of U 2 (0)/U 0 = 0.7 the shape can no longer be approximated with the parameterisation of x 1 and x 2 .
The effect of increasing or decreasing the ratio of the size of the plug regions from h 2 (0)/ h 0 = 0.5 is that the distance it takes for one of the plug regions to disappear becomes smaller. If the velocity ratio is small then, as described earlier, this distance is critical because it marks the point where the flow is sufficiently mixed that it is acceptable to expand thereafter at a wider angle. Hence, the effect of increasing or decreasing h 2 (0)/ h 0 is that this critical distance becomes smaller.
The effects of varying the diffuser expansion ratio h L / h 0 and length ratio L/ h 0 are more obvious and less interesting. Neither of them affect the optimal widening angle, but instead simply make the diffuser continue to expand wider and longer respectively. This is because they don't affect the crucial balance between wall drag and mixing effects.
Varying the upper and lower bounds on the diffuser angle, α max and α min , only affects the optimal solution if the diffuser angle touches the bounds over an interval. For example, in Fig. 11a, c, we see that α never touches α max . Therefore, in this case, raising α max would have no effect on the solution. However, α clearly lies on the lower-bound α min over an interval at the beginning and near the end of the domain. Therefore, varying α min here moves the optimal control along with it.
Throughout this manuscript, we have used a constant value of the spreading parameter S = 0.11. In Sect. 5, we showed that this parameter value is consistent with both a k-and a k-ω SST computational turbulence model using all standard turbulence parameter values [26,27]. In earlier work [6], we compared our simple model to PIV experiments in a three-dimensional geometry and found S = 0.18. The equivalent spreading parameter for free shear layers [39] has been reported to take a range of values for different experiments and CFD calculations. Since the spreading parameter S is associated with shear layer growth rate, for larger S, the shear layer will entrain the plug regions over a shorter distance. Similar to varying h 2 (0)/ h 0 , this decreases the critical distance after which expansion occurs. We have already discussed the effect of S on the singular arc in Sect. 4.2.
For rougher channels with a larger friction factor f , thinner channels and smaller angles will be penalised more. In such cases, the optimal widening angle is larger. Furthermore, if f is sufficiently large, it becomes more advantageous to have an initial widening section, similar to situations where the velocity ratio U 2 (0)/U 0 is close to 1. In the extreme case where wall drag dominates, the control becomes bang-bang because the penalty of worsening the non-uniform flow is eclipsed by the effect of wall drag. For very small wall drag, the optimal diffuser shape appears to prioritise mixing over widening the flow. In these cases, thin channels and small diffuser angles are not penalised very much, such that the critical distance after which expansion occurs is precisely the point where the shear layer has reached across the entire channel. At this point, the flow is sufficiently mixed and can afford expansion.

Conclusions
We have developed a numerical optimisation routine to find the diffuser shape which maximises pressure recovery for given non-uniform inflow, in both two-dimensional and axisymmetric cases. The optimisation uses a simplified mathematical model for the development of turbulent shear layers in confining channels. We find that some of the optimal diffuser shapes are well approximated by shapes which are composed of two straight sections separated by a widening section with a constant widening angle. This is in contrast to diffuser design for uniform flow, where diffusers do not typically have an initial straight section. Furthermore, we show that the optimal widening angle is less than the angle at which boundary layer separation typically occurs, which is usually the diffuser angle chosen for uniform flow. Therefore, we have shown that the effects of non-uniform inflow are critical to diffuser performance, and should not be ignored when it comes to diffuser design.
In two limiting cases, we use analytical techniques to interpret the optimal diffuser shapes found with the numerical optimisation. The first of these cases is the small shear limit, where the inflow is almost uniform, in which case the optimal control is bang-bang, such that the diffuser widens at the maximum possible angle until it reaches the desired cross-sectional area, and then remains at that area. The second case is the pure shear limit, where the inflow is a purely sheared flow with no plug regions. In this case, the optimal control may have a singular arc where, on an interval, the diffuser angle takes values between its upper and lower bounds. We show that in certain cases, the singular arc corresponds to a constant angle, and this angle depends on the friction factor f and the spreading parameter S. We compare some of the numerical optimisation results with CFD simulations using both a k-and a k-ω SST turbulence model (using all standard turbulence parameter values), and the comparion showed good agreement. In the case where we approximate the diffuser shape with piecewise linear sections, we show that both the simplified model and the CFD share almost identical optimal shapes. This suggests that the optimal shapes found using the numerical optimisation and the simplified model are indeed close to the true optimal shapes in reality.

Fig. 12
Comparison between our simple model and a k-CFD model for a two-dimensional diffuser with a stagnation region in the centre. a Velocity colour map calculated using the reduced model, with black dashed lines indicating the shear layer, and a solid black contour indicating the stagnated zone. b Corresponding velocity colour map calculated using the k-model. c Streamlines calculated from the k-model. d Pressure profile averaged across the width of the channel. An alternative colour scheme is used in the colour maps for the purposes of illustrating regions of zero velocity clearly To make the comparison, we choose an invelocity with U 2 (0)/U 0 = 0.75 and h 2 (0)/ h 0 = 0.5. We choose a value of U 2 (0)/U 0 fairly close to 1 to show that stagnation can occur for even moderately non-uniform inflow conditions. The diffuser that we select has constant widening angle α = tan 11 • and has non-dimensional length L/ h 0 = 30. For the CFD, we use the same k-turbulence model as in Sect. 5. In Fig. 12, we compare the results of the model and the CFD. Time-averaged velocity colour maps are compared in Fig. 12a, b, where we indicate the stagnated region in our model (a) with a black contour. The comparison is good, with the model capturing the dominant flow features, such as the width of the shear layer. In Fig. 12c, we display streamlines calculated using CFD. These indicate that there is indeed a stagnated region with recirculation in the centre of the diffuser. This region is located in approximately the same position, and has approximately the same size as the prediction from the mathematical model. The time-averaged pressure profile, averaged across the channel width is plotted in Fig. 12d and, again, it shows good agreement between the CFD and our model, suggesting that our model accurately captures the general behaviour of the diffuser when it has a stagnated region.