Nonlinear two-dimensional free surface solutions of flow exiting a pipe and impacting a wedge

This paper concerns the flow of fluid exiting a two-dimensional pipe and impacting an infinite wedge. Where the flow leaves the pipe there is a free surface between the fluid and a passive gas. The model is a generalisation of both plane bubbles and flow impacting a flat plate. In the absence of gravity and surface tension, an exact free streamline solution is derived. We also construct two numerical schemes to compute solutions with the inclusion of surface tension and gravity. The first method involves mapping the flow to the lower half-plane, where an integral equation concerning only boundary values is derived. This integral equation is solved numerically. The second method involves conformally mapping the flow domain onto a unit disc in the s-plane. The unknowns are then expressed as a power series in s. The series is truncated, and the coefficients are solved numerically. The boundary integral method has the additional advantage that it allows for solutions with waves in the far-field, as discussed later. Good agreement between the two numerical methods and the exact free streamline solution provides a check on the numerical schemes.


Introduction
In this paper, we consider a classical two-dimensional potential free surface flow. Fluid comes in via a pipe, which is placed above a wedge with internal angle β. The incoming flow travels in the direction of gravity. The flow is forcibly separated from the pipe, and a free surface is formed. The flow configuration is shown in Fig. 1. This model is a generalisation of plane bubbles (when β = π ), the subject of a variety of previous studies [1][2][3][4][5][6][7][8]. While some of these works are more analytic in their approach (for example, see [2,3]), the majority of previous investigations make use of a numerical method referred to in literature as series truncation methods. The flow domain is mapped surface tension is of the same order as inertia (for a review, see [30]). For gravity-capillary solutions, it is found that any value of μ is permissible (see Chap. 3 of [30]). This is incorporated into the numerical method by taking μ as an unknown which must be found as part of the solution.
To summarise, in this paper, we compute solutions to the flow geometry shown in Fig. 1 with β ∈ [π/2, π] for three cases: when there is no gravity or surface tension, when there is gravity but no surface tension, and with the inclusion of both gravity and surface tension. An exact free streamline solution is found. Solutions with gravity and surface tension are found using either a series truncation method or a boundary integral method.
The paper is organised as follows. Section 2 is the formulation. Section 3 is the derivation of the exact free streamline solution. Section 4 is a description of the boundary integral method. Section 5 is a description of the series truncation method. Section 6 contains the results of the paper. Section 7 is a conclusion.

Formulation
Consider the problem of flow exiting a pipe and impacting a wedge. We choose Cartesian coordinates, with x pointing in the direction of gravity, and the origin placed at the tip of the wedge. Flow comes in via a two-dimensional pipe of width 2H at x → −∞, where the fluid travels with constant velocity U in the positive x direction. We take H as the reference length and U as the reference velocity. The flow separates where the pipe ends with separation angle μ. Below the pipe is a wedge, placed such that the tip lies beneath the middle of the pipe. There exists a central streamline which meets the wedge at a stagnation point. The angle between the central streamline and the wedge is given by β, and the wedge is at a nondimensional distance W in the x direction below the end of the pipe. We assume the flow is steady, and the flow configuration is shown in Fig. 1.
We denote the velocity vector as u = (u, v), where u and v are the components of the velocity in the x and y directions. We will assume that the fluid is incompressible and the flow is irrotational, and hence we can write the velocity field in terms of two harmonic functions, the velocity potential φ and the streamfunction ψ, defined such that Without loss of generality, we choose φ = 0 at the separation point, and ψ = 0 on the central streamline. We denote the value of φ at the tip of the wedge as φ B . In this scaling, the streamline along the wall at y = 1 and free surface are given by ψ = 1. We also define the complex potential f = φ + iψ. When solving this problem, we will conformally map the flow domain to some auxiliary space with preferable geometry. Via the theory of analytic functions, the governing equation can be handled by demanding the function ξ = u − iv is an analytic function of Fig. 1 Flow configuration in the z-plane. The distance in the x direction from the pipe to the origin is defined as W H so that the corresponding dimensionless quantity is W z = x + iy. The problem then reduces to finding an analytic function such that the relevant boundary conditions are satisfied.
Since the free surface is unknown a priori, as well as the kinematic boundary condition, we require a dynamic boundary condition. The nondimensionalised Bernoulli equation on the free surface (denoted x = η(y)) is given by Here q = |u|, κ is the mean curvature (counted positive when the centre of curvature lies inside the fluid), and F and α are the Froude number and Weber number, respectively, given by Due to the symmetry of the problem about the line y = 0, we can restrict our attention to the problem in 0 ≤ ψ ≤ 1 (y ≥ 0). Recalling that x points in the direction of gravity, the boundary conditions in the z-plane are For the case of β = π , this problem becomes that of flow exiting a pipe. An example of a solution is shown in Fig.  2a. The central streamline ψ = 0 no longer hits a wedge, but becomes a line of symmetry. Since we are dealing with inviscid potential flow, we can take this line of symmetry to be a wall. Furthermore, we can use the symmetry about the streamline ψ = −1. Viewed this way, we see that this problem is then equivalent to that of a plane bubble. The separation angle μ is now the angle between the central streamline and the bubble surface. This problem was solved numerically for finite F both when α −1 = 0 (i.e. no surface tension) and with finite α in a sequence of papers [6][7][8]. The free streamline solution (F → ∞ and α → ∞) solution is simply a uniform stream. When β = π/2, this problem becomes that of flow exiting a pipe onto a flat plate. This was solved numerically when α −1 = 0 by Christodoulides and Dias [11]. We will provide a modification to their numerical scheme to solve the more general problem with arbitrary β. Figure 2(b) shows a typical solution. We also present a new numerical scheme based upon representing the unknowns via an integral equation concerning solely values on the boundary. The necessity for this new method is that it allows for waves in the far-field, a feature solutions with finite surface tension exhibit for this model. One such solution is shown in Fig. 2c. However, first, we will present an exact solution to this problem, under the assumption that both gravity and surface tension are negligible (F → ∞, α → ∞).

Free streamline solution
We will present an exact solution for the flow configuration shown in Fig. 1 under the assumption that gravity and surface tension are both negligible, using a method devised by Love [21]. The method was modified by Hopkinson [22] to allow for internal singularities, and revisited recently by Eggers and Smith [31]. This problem was solved by [24] under the assumption that the wedge was of finite length (see Chap. 2.7). Bernoulli's equation (2) now gives that We restrict our attention to 0 ≤ ψ ≤ 1, noting that the solution for −1 ≤ ψ < 0 is found by reflecting the flow across the streamline AB D. Throughout this paper, points A, B, C, and D will refer to the points as shown in Fig. 3 (A is the flow upstream, B where the two walls meet, C the separation point, and D is downstream). We note that the configuration shown in Fig. 3 is the same as that shown in Fig. 1, but with μ = π . No other value of μ is possible, since if μ < π, then the value of q at the separation point C is zero, while if μ > π, then it is infinite. This can be shown by noting that the local behaviour at the separation point is a flow inside a corner with interior angle μ, which is given by where f and ξ are the complex potential and complex velocity, respectively. From this, one can see that as z → 0, |ξ | → ∞ if μ > π, while if μ < π, ξ → 0. In either case, Eq. (8) cannot be satisfied. We note here that the inclusion of gravity and capillarity allows for different values of μ. Observing that there exist infinite velocities at the separation point when μ > π, it must be that μ ≤ π when surface tension is ignored, or else (7) cannot be satisfied. In fact, it is known that when ignoring surface tension, only three values of μ can satisfy Eq. (7): μ = π/2, 2π/3, and π [25]. The effect of capillarity on the separation angle has been the subject of intense study ( [26][27][28][29][32][33][34], and for a review, [30]). It has been found that when surface tension is included, any value of μ is permissible (that is, local analysis does not restrict the possible values of μ to a discrete set). This is because the infinite velocity occurring when μ > π can be balanced with infinite values of the curvature. Hence, when numerically computing gravity-capillary flows, the value μ must be included as an unknown. We solve the free streamline problem by conformally mapping the problem on to an auxilliary t-plane, which we take to be the lower half-plane, such that all boundaries map onto the real axis. The t-plane is shown in Fig. 4. This mapping is found by guessing a complex potential f (t) which will give us the desired properties of the flow.  Full details can be found in [22], but for this problem, we take f (t) to be a point sink at the origin, given by The streamline AB D is mapped onto the positive real axis, while the streamline AC D is mapped onto the negative real axis in the t-plane. The point C maps to t = −1, and B maps to t = d, where Denoting ξ = e τ −iθ , we define Ω to be the function The flow maps to a semi-infinite strip in the Ω-plane,where the wall B D is given by θ = π − β, and the vertical walls AB and AC are mapped onto the horizontal line θ = 0, as shown in Fig. 5. One finds that the Schwarz-Christoffel mapping from the Ω-plane to the t-plane is given by When β = π , we find that Ω = 0, and the flow is simply a uniform stream. One can then obtain streamlines by integrating the identity along lines of constant ψ. We perform this integration numerically by fixing ψ =ψ and discretising φ. We note that we truncate φ sufficiently far upstream, where φ A is large enough such that the flow is approximately a uniform stream there. Hence, one can write z(−φ A + iψ) = iψ, and then we use the trapezoidal rule to integrate Eq. (14). Care must be taken when integrating at B (and later in the paper at C when μ = π ) due to the stagnation point singularity. For example, from Eq. (9) we know that at B We remove the singularity from the integral by writing The last integral can be evaluated analytically, while we evaluate the first one numerically using the trapezoidal rule. The constant C 1 in this case can be found exactly using Eq. (13). Later in the paper, the value is approximated by matching the solution at the meshpoint closest to the singularity. In the following section, we discuss a boundary integral method that can be used to solve the model with gravity and surface tension.

Boundary integral method
In this section, we will present the first numerical scheme used to solve the model formulated in Sect. 2. The method allows us to compute solutions including the effects of gravity and surface tension. We first express ξ in the form We seek τ ( f ) and θ( f ) as real valued functions of the complex variable f = φ + iψ. As done in Sect. 3, we again map the flow to the lower half t-plane by writing the complex potential in the form of Eq. (10). The method involves using Cauchy's integral formula to express the unknowns in terms of boundary values. Writing t = γ + iδ, we relate values of τ in the f -plane with values in the t-plane by writing with a similar equation forθ . The boundary conditions (4-6) can then be written in the form Applying Cauchy's integral formula to the functionτ − iθ in the t-plane, where we take a contour consisting of the real axis and a semi-circle of arbitrarily large radius in the lower half t-plane, one can show that where γ 0 is a point on the real axis in the t-plane, and the integrals are evaluated in the Cauchy principle value sense. The real component of the above equation giveŝ The boundary conditions (19)(20)(21) reduce Eq. (23) tô The right-hand side of the above integral equation concerns only values ofθ on the free surface (γ ∈ [−1, 0]). Hence, setting γ 0 ∈ [−1, 0], we can find values ofτ on the free surface if we know values ofθ . We note that the integral has a removable singularity at γ = γ 0 . We instead write We perform a change of variables on the integral, such that it is expressed in terms of φ. This results in where γ 0 = −e −πφ 0 . Bernoulli's equation (2) now reads It suffices to find θ(φ) and τ (φ) that satisfy (26) and (27), which is done numerically, as described below.
where φ D must be chosen to be sufficiently large such that the solution becomes invariant to further increase in φ D . We discretise φ into N equally spaced points, given by We We take as unknown N values of θ(φ I + i), which we shall denote θ I . One can then find values of θ(φ M I + i), denoted θ M I , via fourth-order interpolation formulae. We note that the separation angle μ is related to θ : One can find τ at the midpoints φ M I by evaluating Eq. (26) at φ M I . The integral is evaluated numerically using the trapezoidal rule. We satisfy Bernoulli's equation at the midpoints φ M I , resulting in N − 1 equations. When β = π , we also fix W , the height of the pipe (we describe shortly how the value of W is recovered). Thus far, satisfying the boundary integral equation (26) and Bernoulli's equation ensure τ and θ are analytic, and satisfy the boundary conditions. The final equations required to close the system depend on the type of solution sought. For example, if seeking a pure gravity solution with β = π/2, then we have N + 2 unknowns: θ 1 , . . ., θ N , φ B , and B. Given we already have N equations, the two equations which close the system are fixing μ, and setting θ N = θ N −1 (such that the far-field is flat). If we seek gravity-capillary solutions with β = π/2 (these solutions have waves in the far-field), it is found that more quantities have to be fixed to obtain a unique solution. In this case, the N + 1 unknowns are θ I and B. The additional equation is given by fixing μ. A more detailed discussion on what must be fixed and be allowed to vary for each case is given in Sect. 6.
In the above, we need values of x and ∂θ/∂φ on the interface to satisfy Bernoulli's equation. Furthermore, we need to evaluate W . Values of ∂θ/∂φ are found using four-point finite difference formulae on the discretised solution θ I . Values of x and W are found as described below. The values of x I = x(φ I + i) are evaluated by integrating the identity We find that This integral is evaluated numerically using the midpoint rule, ensuring to remove the singularity at C from the integral. Recall that locally, the flow behaves like This gives where the constant C 2 is approximated by equating the local behaviour (33) to the numerical solution (31) at the first midpoint φ M 1 . The height of the pipe W is found in the same manner as discussed in Sect. 3. We numerically integrate Eq. (31) along the streamlines ψ = 0 and ψ = 1. We again truncate φ → −∞ to φ = −φ A , ensuring φ A is sufficiently large that further increase in φ A results in negligible change in the value of W . We then discretise φ ∈ [−φ A , φ B ] for ψ = 0 into M + 1 equally spaced points φ I via and φ ∈ [−φ A , 0] for ψ = 1 into L + 1 equally spaced points φ I via where k = (φ A + φ B )/M and l = φ A /L. We know the values of θ on the boundaries ψ = 0 and ψ = 1. The values of τ ( f ) can be found on AB (ψ = 0, φ < φ B ) via the identity Meanwhile, the values of τ can be found on AC (ψ = 1, φ < 0) via the identity The above integral is written such that the singularity that occurs as φ 0 → 0 is removed. The value W is then given by The above integrals are computed numerically via the midpoint rule, where the singularities at B and C are again removed. When computing solutions, we must ensure that the number of unknowns matches the number of equations.
The number of parameters we require to fix to obtain a unique solution depends on the parameter values considered, as discussed in Sect. 6. This concludes the description of the boundary integral method. In the following section, we discuss the series truncation method.

Series truncation method
In this section, we will present a second numerical scheme based on series truncation to compute fully nonlinear solutions to the system of equations (4)- (7). As with the method from the previous section, it also allows for the inclusion of gravity and surface tension. However, one weakness is that it can only compute solutions which far downstream approach a uniform stream. Hence, this method cannot compute solutions with gravity and surface tension when β = π/2, since such solutions are found to approach a periodic wavetrain downstream. For a review of series truncation methods applied to steady potential flow, see §3 of [35]. This section will follow closely the work of [6,7], who solved this problem for β = π with surface tension and gravity, and [11], who likewise devised a method to find solutions with gravity for β = π/2. The method once again involves conformally mapping the flow domain to an auxiliary plane. This time we choose the s-plane to be a unit semi-circle (in the upper half-plane). Given that the flow domain in the f -plane is given by an infinite strip 0 ≤ ψ ≤ 1, the mapping from f to s is found to be The f -plane and s-plane are shown in Figs. 6 and 7, respectively. The point A maps to s = 0, C to s = −1, D to s = 1, and B to s = s B , where The reflected portion of the flow, that is −1 ≤ ψ ≤ 0, is mapped to the semi-circle in the lower half-plane obtained when reflecting the image in Fig. 7 across the real axis. To ensure that the series representation of ξ converges, we must remove all singularities of ξ in {s : |s| ≤ 1, Im(s) ≥ 0}. The first singularity we will consider comes from the separation point C. The flow here behaves like the flow in a corner of angle μ. Using Eqs. (9) and (40), one finds Next, we consider the singularity at the tip of the wedge B. Similar to Eq. (42), the leading order singular behaviour is that of flow in a corner of interior angle β. Again, using Eqs. (9) and (40), one finds The final singularity is associated with the flow in the far-field D. The singularity has two possible behaviours. First, when β ∈ (π/2, π), it was shown by Birkhoff and Carter [1] that Using Eq. (40), we find ξ ∼ (− log(1 − s)) 1/3 as s → 1.
If β = π/2, the far-field can approach either a uniform stream or a periodic train of waves. Assume for now that the flow approaches a uniform stream with velocity U f . The method in this section requires ξ to be defined at s = 1, and hence can only solve for solutions with flat far-fields. The downstream Froude number F f and Weber number α f relate to F and α via the equations It is found that where λ is a root of the equation Here, B f is the downstream Bond number, given by B f = F 2 f /α f . The above behaviour is found by linearising the flow about a uniform stream with speed U f . It is the dispersion relation for gravity-capillary waves, with a wavelength of iλ (and hence a wavenumber k = 2π i/λ). Note that the far-field has oscillations unless λ is real (see Eq. (47)). Equations (10) and (47) imply that Therefore, we propose two different representations for ξ . In the first instance, when β = π/2, we write where C is a fixed constant satisfying 0 < C < 0.5, such that ξ is purely real on the line s ∈ [−1, s B ]. It can be seen that (50) satisfies the singular behaviour (42), (43), and (45), and that ξ(0) = 1. This series, with β = π , was used in [7,8,36]. If β = π/2, then we instead take where A and λ are unknown, and have to be found as part of the solution. We note here that with β = π/2, no waveless solutions were found when F f = 0 and B f = 0. Hence, the series (51) can only be used to compute solutions with no surface tension. When B f = 0, then from Eq. (48) we find that λ is real for F f > 1 (i.e. solutions are waveless when the flow in the far-field is super-critical). The constant λ is taken to be the smallest positive root of (48). We note that [11] used a similar series representation for ξ when considering β = π/2, except they did not include the term exp(A(1 − s) 2λ ) in their expansion. They instead used We found that, truncating the series after a sufficient number of coefficients, their results agree within graphical accuracy with those computed using the series (51). However, the coefficients a n decay at a much slower rate without this new exponential term, as discussed in Sect. 6.3. We truncate both series after N − 1 terms. When β = π/2, this results in N + 1 unknowns (a 1 , . . ., a N −1 , μ, B). The free surface in the s-plane is given by s = e iσ with σ ∈ [0, π]. We discretise σ into N + 1 equally space points, given by The series representation of ξ has been constructed such that the kinematic boundary conditions and the governing equation are satisfied. It is left to satisfy the dynamic boundary condition (7) at each meshpoint σ I , which is given by This gives us N + 1 equations for N + 1 unknowns, which are solved via Newton's method. When β = π/2, we have two additional unknowns, A and λ. Therefore, we satisfy Eq. (7) at N + 2 equally spaced points in σ , given by We also satisfy Eq. (48), where F f is calculated using Eq. (46), and U f can be found by evaluating (51) at s = 1. Hence, we have obtained N + 3 equations for N + 3 unknowns, a system which can be solved numerically using Newton's method. The profile of the interface is obtained by numerically integrating the identity We do this using the trapezoidal rule on the discritised domain σ I . The height of the pipe W is found in the same manner as was discussed in Sect. 4. We numerically integrate equation along the streamlines ψ = 0 and ψ = 1 from φ = −φ A to φ = φ B and φ = 0, respectively. Values of ξ can be found for a given φ and ψ by finding the the corresponding value of s from Eq. (40), and then using the relevant series solution for ξ . The same discretisation and numerical integration as was done in Sect. 4 is used to compute W . This concludes the description of the series truncation method. In the following section, we discuss the results of the paper.

Results
We begin by briefly recapping the results when β = π , which are found in [6][7][8], since the results here help inform us on what to expect for the solution space for other values of β.

Results for β = π
First, consider the case when α −1 = 0. The flow configuration is as shown in Fig. 2a. The model corresponds to flow exiting a pipe with no obstruction, or plane bubbles rising at a constant velocity. For each value of F, there is a unique solution. It is found that there exists a critical Froude number F C ≈ 0.51, such that when F < F C the separation angle μ = π/2, while if F > F C , then μ = π , as illustrated in Table 1. The solutions with μ = π/2 are referred to as smooth bubbles, while μ = π solutions are cusped bubbles. There is also a unique pointed bubble with μ = 2π/3 when F = F C . When surface tension is nonzero (that is α −1 = 0), the angle μ becomes continuously dependent on F. This is shown in Fig. 8 for α = 5 and α = 20. It can be seen that as F decreases, the value μ oscillates about μ = π/2. Furthermore, as α is increased, the amplitude and wavelength of the oscillations decrease. Let us denote F 1 (α) as the largest value of F where the curves intersect for a given α, F 2 (α) as the second largest such value, and so on. These branches F i (α) are branches of smooth bubbles, and are monotonically increasing with α. As α → ∞, all of these branches converge onto a single value, F = F * . It is found that F * ≈ 0.318. The mechanism by which a unique solution is chosen from a possible set of solutions, by including surface tension (or other forces), and taking said forces to zero, is known as solution selection. A similar phenomenon is seen in Saffman-Taylor fingering [37,38]. The solutions on the primary branch F 1 (α) smooth profiles with no oscillations. Fig. 8 Relationship between μ and F for α = 5 (solid curve) and α = 20 (dashed curve) for β = π . The dotted lines are μ = π/2 and μ = π Table 1 A table for the value of μ when α −1 = 0 On the higher order branches of smooth bubbles, the profiles develop oscillations [39]. It is of interest to note that the solution space of axisymmetric Taylor bubbles exhibit much of the same behaviour [39,40].

Results for β ∈ (π/2, π)
The flow configuration for β ∈ (π/2, π) is shown in Fig. 2b. The results are qualitatively similar to the results found for β = π . First, consider the case when α −1 = 0. For a fixed β, we now have two free parameters, F and φ B . In the computations, instead of fixing φ B and recovering W a posteriori, we chose to fix W and allow φ B to vary. Note that the parameter φ B is related to s B in the series truncation method via Eq. (41), and d in the boundary integral method via Eq. (11). Unsurprisingly, it is found that φ B monotonically increases as W is increased. When surface tension is ignored, we again see a critical value F C , where the solutions have the properties given in Table 1. It is now the case that F C has dependence on β and W . We conjecture that as W → ∞, F C has a finite limiting value equivalent to the F C found by Vanden-Broeck [8] for plane bubbles. This comes from the fact that as W increases, the pipe is moved further away from the wedge, and the behaviour of the flow near the separation point becomes closer to that of flow exiting a pipe without an obstruction (β = π ). The dotted and dashed curves in Fig. 9 show the dependence of F c on W for β = 2π/3 and β = 5π/6, respectively. We experience difficulties with computing F = F C solutions for values of W greater than that shown in the figure due to our inability to get enough meshpoints distributed along the free surface for these more extreme profiles. Streamlines to typical solutions are shown in Fig. 10 for β = 2π/3 and W = 1. The bold lines in the figure correspond to solid boundaries, the dashed curves interior streamline, and the solid curve the free surface. Figure  - Fig. 10c is for F = 20, one would expect reasonable agreement with the free streamline solution (F → ∞). It can be seen that the agreement is very good near the separation point. However, as one moves further along the profile, the solutions start to deviate, due to the different singular behaviour downstream (the free streamline solution approaches a uniform stream, as opposed to the behaviour (45)). This is shown in Fig. 10d, where we have plotted the flow further downstream for the F = 20 solution. The solutions presented are computed using the series truncation method described in Sect. 5. However, we also computed the solutions using the boundary integral method from Sect. 4, as shown by the circles in figure (a) and (b). As can be seen, the profiles almost overlap. The agreement between the analytic solution and two numerical schemes provides a convincing check on the methods used. When surface tension is included, similar results to the plane bubble are found. One finds that the separation angle μ has a continuous dependence on F. For a given α, the value of μ oscillates about π/2 for values of F < F C . Meanwhile, when F → ∞, μ → π . This is the same behaviour as shown in Fig. 8 for β = π . One could again use the solution selection procedure to find a unique solution to the zero surface tension problem with a separation angle of μ = π/2. However, there is no clear physical significance to the selected solution. Another similarity with the β = π solution space is that the solutions on the higher mode branches also develop oscillations on the free surface, as shown in Fig. 11. Profiles of solutions with β = 2π/3, W = 0.5, and α = 5 are shown. The values of F are chosen such that μ = π/2.  Finally, we shall discuss the solution space when β = π/2. This problem was already solved numerically using a series truncation method when surface tension is ignored [11]. However, as mentioned previously in Sect. 5, the authors failed to remove an additional singularity caused by the far-field decay of the free surface (see Eq. (47)). The results presented in the paper agree with our calculations within graphical accuracy. However, the decay of the coefficients of the truncated series a n is weaker, as shown in Table 2. The table shows the order of the coefficients for the solutions with φ B = 1.5 for two values of F. It can be seen this additional singularity is important to remove for convergence of the series representation of ξ . It was found by [11] that, as with values of β in the range β ∈ (π/2, π], fixing F and W results in a unique solution. In this instance, the far-field is not an infinitesimal jet but instead a uniform stream. There again exists a critical value F C , dependent on W , such that the separation angle μ agrees with the behaviour seen in Table 1. The relation between F C and W is shown by the solid curve in Fig. 9. The solution space when surface tension is included is more complicated than when β ∈ (π/2, π]. We find that fixing F and α is not sufficient to produce a unique solution. We find we must also fix two additional parameters, which we choose to be μ and φ B . Solutions with waves in the far-field are witnessed. We denote the amplitude of the waves in the far-field as ω. We compute some solution branches, fixing F, α, φ B , and W , and then allowing μ to vary. Having fixed the parameters as described above, it is found that solutions do not exist for all values of μ. Instead, solutions exist within a range of values of μ, the range depending on the fixed parameters. Once one solution has been found, one can find the complete solution branch via continuation in μ. The solution branches start and end on solutions with larger ω. For all parameters tested, starting from the smallest value of μ for which there is a solution, ω monotonically decreases to a minimum, and then monotonically increases to the largest possible value of μ. This can be seen in Fig. 12, where we show some typical solution branches for various values parameter  Fig. 13 The profiles corresponding to the points a-i in Fig. 12 values. The branches are plotted in (μ, ω) space. Some typical profiles are shown in Fig. 13. It was found for the parameters tested that the branch never touched the ω axis. Hence, no waveless gravity-capillary profiles were found.
To check the convergence of the solutions with waves in the far-field, we vary both the mesh spacing h and domain truncation φ D . It can be seen in Fig. 14 that the agreement upon mesh refinement is good. We also notice that the last wavelength on the figure is spurious: this is expected, due to the truncation of the computational domain. This ensures that domain truncation does not cause errors to the flow close to the separation point, and that our numerical approximation of a finite length far-field is valid, as long as we compute enough wavelengths down from the pipe.

Conclusion
In conclusion, we have provided an extensive numerical study of the model shown in Fig. 1. We have considered three case. First, when surface tension and gravity are negligible, we derived an exact free streamline solution. Next, when gravity is included but surface tension is neglected, we recovered previously found results for β = π and β = π/2, as well as finding results for β ∈ (π/2, π). Finally, when surface tension is also included, we recomputed the results for β = π , while the results for β ∈ [π/2, π) are novel. The case when β = π/2 has a richer solution space, since the far-field condition allows for waves to occur on the free surface. Good agreement between the two numerical methods employed, as well as the exact free streamline solution, provides a check on the accuracy of the numerical results presented.