2.1.1. Inflow, outflow and wall boundary conditions

Along the channel walls no slip and no penetration conditions give us

(2.6)\begin{equation} u(x,\pm{1})=v(x,\pm{1})=0,\quad \text{for all}\ x\leq 0, \end{equation}

and along the line of symmetry a vanishing shear stress and no penetration mean that

(2.7)\begin{equation} \frac{\partial u}{\partial y}(x,0) = v(x,0) =0. \end{equation}

By choosing the constant pressure gradient far upstream to be such that the volume flux of fluid in the half-channel is equal to, say, $1$, we get, assuming fully developed flow ($u=u(y)$ and $v=0$ as $x\rightarrow -\infty$) the usual parabolic velocity profile

(2.8)\begin{equation} u(x,y)=\frac{3}{2}(1-y^2), \ v(x,y)=0,\quad \text{as}\ x\rightarrow -\infty. \end{equation}

This immediately leads to the conclusion that

(2.9a,b)\begin{equation} \psi=\frac{3}{2}y\left(1-\frac{y^2}{3}\right) \quad\text{and}\quad \omega ={-}3y, \end{equation}

as $x\rightarrow -\infty$, where we have chosen to fix $\psi =0$ on $y=0$. The line of symmetry and the solid walls $y=\pm 1$ are streamlines of the flow (see § 2.1.2). As $x\rightarrow \infty$, uniform flow conditions are assumed to hold and therefore $\lim _{x\rightarrow \infty }u(x,y)=1/\eta _\infty$ and $\lim _{x\rightarrow \infty }v(x,y)=0$, where $\eta _\infty :=\lim _{x\rightarrow \infty }\eta (x)$. Thus, the fluid vorticity vanishes infinitely far downstream and

(2.10)\begin{equation} \lim_{x\rightarrow\infty}\psi(x,y)=\frac{y}{\eta_\infty}, \end{equation}

the free surface being a streamline, with $\psi (x,\eta (x))=1$.

2.1.2. Free surface conditions

We require that the free surface $\varGamma$ does not move, that the free surface shear stress vanishes and that there is a zero net normal force at each point on $\varGamma$.

Let $\hat {\boldsymbol {n}}$ denote the outward pointing unit normal vector and $\hat {\boldsymbol {t}}$ a unit tangent vector to the boundary of $\varOmega$, oriented in the clockwise direction. Then the kinematic condition imposed on the free surface is

(2.11)\begin{equation} \hat{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{v}=0, \end{equation}

which is equivalent to the free surface being a streamline ($\psi (x,\eta (x))=1$ in the present case).

Denoting the Cauchy stress tensor by

(2.12)\begin{equation} \boldsymbol{\sigma}:={-}p\boldsymbol{\delta} + (\boldsymbol{\nabla}\boldsymbol{v}+(\boldsymbol{\nabla}\boldsymbol{v})^{\rm T}), \end{equation}

where $\boldsymbol {\delta }$ is the identity tensor, the vanishing of the shear stress on the free surface may be written as

(2.13)\begin{equation} \sigma_{nt}:=\hat{\boldsymbol{t}}\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot}\hat{\boldsymbol{n}} = 0. \end{equation}

From (2.3a,b) and the definition of $\boldsymbol {\sigma }$ in (2.12) we get

(2.14)\begin{equation} \frac{\partial^2\psi}{\partial n^2} - \frac{\partial^2\psi}{\partial t^2} =0, \quad\text{on}\ \varGamma, \end{equation}

and thence, on the free surface,

(2.15)\begin{equation} \omega:= \Delta\psi = \frac{\partial^2\psi}{\partial n^2} + \frac{\partial^2\psi}{\partial t^2} = 2\frac{\partial^2\psi}{\partial t^2} ={-}2\kappa\frac{\partial\psi}{\partial n}, \end{equation}

(see, for example, Batchelor (Reference Batchelor1967) or Longuet-Higgins (Reference Longuet-Higgins1953)), where $\kappa$ denotes the signed curvature of the free surface

(2.16)\begin{equation} \kappa(x) ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\boldsymbol{n}}. \end{equation}

The relationship (2.15) will be seen to be useful in § 3.

Finally, the normal stress balance equation that must be satisfied on $\varGamma$ is

(2.17)\begin{equation} \left[\sigma_{nn}\right]:=\hat{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot}\hat{\boldsymbol{n}} + p_a = \gamma\kappa, \end{equation}

where $\gamma$ denotes the dimensionless surface tension (the reciprocal of the capillary number) and $[\sigma _{nn}]$ denotes the jump $\sigma _{nn} + p_a$ where $p_a$ is atmospheric pressure, assumed constant. Again, using (2.3a,b) and the definition of $\boldsymbol {\sigma }$ in (2.12), we are led to the result

(2.18)\begin{equation} [\sigma_{nn}] ={-}[p]-2\frac{\partial^2\psi}{\partial n\partial t} ={-}[p]-2\frac{\partial^2\psi}{\partial n\partial s} = \gamma\kappa, \quad \text{on}\ \varGamma. \end{equation}

Here, $\partial /\partial s$ denotes differentiation with respect to the arclength $s$ and $[p]$ is the excess pressure (fluid pressure $p$ minus $p_a$).

2.2.1. Shape of the free surface near $(0,1)$

Michael (Reference Michael1958) assumed, at least sufficiently close to the separation point $(0,1)$, that the equation of the free surface $\varGamma$ was $\theta =\alpha$ in plane polar coordinates $(r,\theta )$ centred at $(0,1)$ with $\theta$ measured in the anticlockwise direction from the channel wall $y=1$; see figure 2. We now seek to generalize this by assuming only the existence of a right-hand tangent line to the free surface at the separation point $(0,1)$, and making an angle $\alpha$ with the channel wall $y=1$. Thus, $\alpha$ is the limit, as $r\rightarrow 0$, of $\theta$, when $(-r\cos \theta,1-r\sin \theta )$ is a point on the free surface.

Figure 2. Polar angle $\theta$ measured in the anticlockwise direction from the channel wall $y=1$.

Drawing inspiration from section A.4 of the appendix of Ramalingam (Reference Ramalingam1994) but with some differences of notation from his, we write the equation of the free surface $\varGamma$ near $(0,1)$ as the sum of the equation of the right-hand tangent line at $(0,1)$ and terms accounting for the possible nonlinearity of the free surface

(2.19)\begin{align} y & = 1+h(r) \nonumber\\ & = 1+ r\sin(\alpha-{\rm \pi}) + h_1r^{p_1+1} + h_2r^{p_2+1}-\frac{c_pr^2}{2\gamma}\cos\alpha +O(r^{p_3+1}) \nonumber\\ & = 1 - r\sin\alpha + h_1r^{p_1+1} + h_2r^{p_2+1}-\frac{c_pr^2}{2\gamma}\cos\alpha +O(r^{p_3+1}), \end{align}

with $p_3>1,\ p_2>p_1>0$ and $0< \gamma <\infty$, as shown in figure 1, where $r$ denotes the distance from $(0,1)$. The coefficients $h_1$ and $h_2$ and exponents $p_1$ and $p_2$ are to be determined from the kinematic condition (2.11), the vanishing of the shear-stress (2.13) and the normal stress balance (2.18) on the free surface. Here $c_p$ is a constant (actually, a function $c_p=c_p(\gamma )$) appearing in the singular expression for the pressure (see § 2.2.6) and chosen so as to make the singular normal stress and the far-field normal stress continuous on the free surface. Note that in his thesis, Ramalingam (Reference Ramalingam1994) did not attempt to calculate $h_2$ and that the term involving $c_p$ is absent.

2.2.2. Singular solution in a sector centred at $(0,1)$

In the sector

(2.20)\begin{equation} \mathcal{D} =\{(r,\theta): 0\leq{r}\leq{r_c},\; 0\leq \theta \leq {\rm \pi}+ \arcsin(h(r)/r)\}, \end{equation}

at the separation point $(0,1)$ (for a certain $r_c>0$) we write the stream function in the form

(2.21)\begin{equation} \psi = 1 + \sum_{n=1}^\infty r^{1+\lambda_n}f_n(\theta), \end{equation}

where the parameters $\{\lambda _n\}_{n=1}^\infty$ are eigenvalues to be determined and the non-zero function $f_n(\theta )$ is a linear combination of trigonometric functions

(2.22)\begin{equation} f_n(\theta) = A_n\cos((\lambda_n+1)\theta) + B_n\sin((\lambda_n+1)\theta)+C_n\cos((\lambda_n-1)\theta)\!+\!D_n\sin((\lambda_n-1)\theta).\end{equation}

Lugt & Schwiderski (Reference Lugt and Schwiderski1965) showed that the functions $r^{1+\lambda _n}f_n(\theta )$ ($n=1, 2, 3, \ldots$) in (2.21) form a complete set of solutions to the biharmonic equation in the plane, provided that $\lambda _n\not =1$ and $Re(\lambda _n)>0$. In all that follows we have ordered the $\{\lambda _n\}$ so that

(2.23)\begin{equation} 0< Re(\lambda_1)< Re(\lambda_2)< Re(\lambda_3)< \ldots . \end{equation}

The radial and transverse components of the velocity $\boldsymbol {v}$ are given by

(2.24a,b)\begin{equation} v_r = \frac{1}{r}\frac{\partial\psi}{\partial\theta} = \sum_{n=1}^\infty r^{\lambda_n}f_n'(\theta) \quad \text{and}\quad v_\theta ={-}\frac{\partial\psi}{\partial r} ={-}\sum_{n=1}^\infty (\lambda_n+1)r^{\lambda_n}f_n(\theta), \end{equation}

so that, as remarked for example by Sturges (Reference Sturges1979), the weak regularity condition $Re(\lambda _n)>0$ may be seen to mean that the velocity vanishes as $r\rightarrow 0$. Following Michael (Reference Michael1958), we also observe that $Re(\lambda _n)>0$ is a physical requirement in order that the integrated stresses remain finite as $r\rightarrow 0$.

Note that additional terms of the form

(2.25)\begin{equation} r(\mathcal{A}_0\cos\theta + \mathcal{B}_0\sin\theta + \mathcal{C}_0\theta\cos\theta + \mathcal{D}_0\theta\sin\theta) \end{equation}

and

(2.26)\begin{equation} r^2(\mathcal{A}_1\cos(2\theta) + \mathcal{B}_1\sin(2\theta) + \mathcal{C}_1\theta + \mathcal{D}_1), \end{equation}

could have been added to the series (2.21) (being solutions to the biharmonic equation) but these are excluded. The terms in (2.25) are omitted for reasons that are the same as those given above for the insistence that the weak regularity condition $Re(\lambda _n)>0$ hold. For the terms in (2.26), it may be shown that if the boundary conditions applied are the same as those for $f_1(\theta )$ in §§ 2.2.3–2.2.5, non-trivial solutions for the coefficients $\mathcal {A}_1$, $\mathcal {B}_1$, $\mathcal {C}_1$ and $\mathcal {D}_1$ are possible only if $\alpha$ satisfies

(2.27)\begin{equation} \sin(2\alpha)-2\alpha\cos(2\alpha)=0, \end{equation}

(see too, (2.31) of Anderson & Davis (Reference Anderson and Davis1993) or (32) of Lugt & Schwiderski (Reference Lugt and Schwiderski1965)) and therefore only for very specific values of the separation angle: $\alpha \approx 0.715{\rm \pi}$, $1.230{\rm \pi}$ or $1.735{\rm \pi}$.

If the equation $y=1+h(r)$ of $\varGamma$ is linear in $r$, as assumed by Michael (Reference Michael1958), then obviously conditions on this surface are imposed at $\theta =\alpha$. More generally, however, in order to impose (kinematic and stress) conditions at a point $(r,\theta )$ on the free surface (2.19) we need to be able to develop all pertinent functions of $\theta$ about $\alpha$. We note from (2.19) and $h=r\sin (\theta -{\rm \pi} ) = -r\sin \theta$ that

(2.28)\begin{equation} \sin\theta=\sin\alpha-h_1r^{p_1} -h_2r^{p_2} + \frac{c_pr}{2\gamma}\cos\alpha + {\rm h.o.t.}, \end{equation}

where h.o.t. refers to higher-order terms. Choosing $(r,\theta )$ to be the polar coordinates of a point on the free surface $\varGamma$, and recalling that $\alpha$ is the limit in this case as $r\rightarrow 0$ of $\theta$, we have $\sin (\theta -\alpha ) \sim \theta -\alpha$ which means, using (2.28) and a Taylor series centred at $\alpha$, that

(2.29)\begin{align} \cos\theta & =\cos\alpha-(\theta-\alpha)\sin\alpha + {\rm h.o.t.} \nonumber\\ & = \cos\alpha -\sin\alpha\left[\left(\sin\alpha-h_1r^{p_1} -h_2r^{p_2} + \frac{c_pr}{2\gamma}\cos\alpha\right)\cos\alpha - \sin\alpha\cos\theta\right]\nonumber\\ &\quad + {\rm h.o.t.}\nonumber\\ \Rightarrow \cos\theta & = \cos\alpha + \tan\alpha\left(h_1r^{p_1} +h_2r^{p_2} -\frac{c_pr}{2\gamma}\cos\alpha\right) + {\rm h.o.t.} \end{align}

More generally, any function $F(\theta )$, differentiable at $\theta = \alpha$ (for example, $f_n$ or $f_n'$), may be written as a Taylor series centred at $\alpha$ as follows:

(2.30)\begin{equation} F(\theta) = F(\alpha) -\frac{1}{\cos\alpha}\left(h_1r^{p_1} + h_2r^{p_2} - \frac{c_pr}{2\gamma}\cos\alpha\right)F'(\alpha) + {\rm h.o.t.} \end{equation}

In what follows, we describe the boundary conditions on the solid wall $\{(x,1): x\leq 0\}$ as well as the kinematic, shear-free and normal stress balance conditions that must hold on the free surface $\varGamma$.

2.2.3. Wall boundary conditions (2.6)

From (2.6) and (2.24a,b), no slip and no penetration on the solid wall $\theta =0$ lead to the conditions

(2.31)\begin{equation} f_n(0) = f_n'(0) =0, \quad n=1,2,3, \ldots \end{equation}

2.2.4. Kinematic condition (2.11)

Using (2.19) and (2.24a,b), the vanishing of the normal component $\hat {\boldsymbol {n}}\boldsymbol {\cdot }\boldsymbol {v}$ of velocity on the free surface $y=1+h(r)$ yields,

(2.32)\begin{align} &\sum_{n=1}^\infty \left(h'\cos\theta - \frac{r-hh'}{\sqrt{r^2-h^2}}\sin\theta\right)r^{\lambda_n}f_n'(\theta) \nonumber\\ & \quad + \sum_{n=1}^\infty\left(h'\sin\theta + \frac{r-hh'}{\sqrt{r^2-h^2}}\cos\theta\right)(\lambda_n+1)r^{\lambda_n}f_n(\theta) = 0. \end{align}

Developing all functions of $\theta$ in (2.32) about $\alpha$, as explained in (2.28)–(2.30), the leading-order term for the first line of (2.32) may be shown to be

(2.33)\begin{equation} \left(\frac{h_1p_1r^{p_1}}{\cos\alpha}\right)r^{\lambda_1}f_1'(\alpha), \end{equation}

whereas that for the second line of (2.32) is $-(\lambda _1+1)r^{\lambda _1}f_1(\alpha )$. To $O(r^{\lambda _1})$, the kinematic condition therefore gives

(2.34)\begin{equation} f_1(\alpha)=0, \end{equation}

in agreement with (A.27) of Ramalingam (Reference Ramalingam1994).

Suppose now that $p_1= \lambda _2-\lambda _1$, as assumed by Ramalingam (Reference Ramalingam1994). Then, equating the $O(r^{\lambda _2})$ terms of (2.32) gives

(2.35) \begin{align} \frac{h_1p_1}{\cos\alpha}f_1'(\alpha) & ={-}\frac{h_1}{\cos\alpha}(\lambda_1+1)f_1'(\alpha) + (\lambda_2+1)f_2(\alpha),\notag\\ \Rightarrow f_2(\alpha) & = \frac{h_1f_1'(\alpha)}{\cos\alpha}. \end{align}

and the expression (2.35) corrects the result in (A.28) of Ramalingam (Reference Ramalingam1994).

2.2.5. Vanishing free surface shear stress (2.13)

The evaluation of the leading-order ($O(r^{\lambda _1-1})$) terms in (2.13) leads (again, using (2.28)–(2.30)) to

(2.36)\begin{equation} \frac{1}{2}((1-\lambda_1^2)f_1(\alpha) +f_1''(\alpha)) = 0, \end{equation}

which, combined with (2.34), gives

(2.37)\begin{equation} f_1''(\alpha) =0, \end{equation}

as also obtained by Ramalingam (Reference Ramalingam1994) in his (A.30). From (2.22), (2.31), (2.34) and (2.37) we have (cf. Michael (Reference Michael1958), equations (1)–(2))

(2.38)$$\begin{gather} A_n+C_n=0,\quad n=1,2,3,\ldots, \end{gather}$$

(2.39)$$\begin{gather}(\lambda_n+1)B_n + (\lambda_n-1)D_n = 0, \quad n=1,2,3,\ldots, \end{gather}$$

and, when $n=1$, (cf. Michael (Reference Michael1958), equations (3)–(4))

(2.40)\begin{align} &A_n\cos((\lambda_n+1)\alpha) + B_n\sin((\lambda_n+1)\alpha)+C_n\cos((\lambda_n-1)\alpha)+D_n\sin((\lambda_n-1)\alpha) = 0, \end{align}

(2.41)\begin{align} &A_n(\lambda_n+1)^2\cos((\lambda_n+1)\alpha) + B_n(\lambda_n+1)^2\sin((\lambda_n+1)\alpha)\nonumber\\ &\quad+C_n(\lambda_n-1)^2\cos((\lambda_n-1)\alpha)+D_n(\lambda_n-1)^2\sin((\lambda_n-1)\alpha) =0. \end{align}

We see, from calculation of the determinant of the system (2.38)–(2.41) when $n=1$ that for non-trivial solutions to exist in this case, $\lambda _n$ has to satisfy the eigenvalue problem

(2.42)\begin{equation} \lambda_n\sin{2\alpha} - \sin{2\alpha\lambda_n}=0, \end{equation}

and, for $\alpha \in [{\rm \pi}, 3{\rm \pi} /2]$, this equation always has a root $\lambda _1\in [1/3,1/2]$. Equation (2.42) has already been obtained by Dean & Montagnon (Reference Dean and Montagnon1949), Lugt & Schwiderski (Reference Lugt and Schwiderski1965) and Moffatt (Reference Moffatt1964) in the context of corner and wedge flows, and, for example, by Sturges (Reference Sturges1979) and Ramalingam (Reference Ramalingam1994) in the present context. Huilgol & Tanner (Reference Huilgol and Tanner1977) also derived (2.42) in their study of the separation at a sharp edge of a second-order fluid under creeping flow conditions.

From (2.38)–(2.42), three of the coefficients $A_1$, $B_1$, $C_1$ or $D_1$ in the development (2.22) of $f_1(\theta )$ may be expressed in terms of the fourth. Expressing $A_1$, $B_1$ and $C_1$ in terms of $D_1$, for example, leads to an expression for $r^{\lambda _1+1}f_1(\theta )$ identical to that for the leading-order perturbation approximation to the stream function in (2.22) of Anderson & Davis (Reference Anderson and Davis1993).

Again, assuming that $p_1 = \lambda _2 - \lambda _1$ and using (2.28)–(2.30), the $O(r^{\lambda _2-1})$ terms in (2.13) give

(2.43)\begin{equation} 2\lambda_1f_1'(\alpha)\frac{h_1p_1}{\cos\alpha} - \frac{(1-\lambda_1^2)h_1}{2\cos\alpha}f_1'(\alpha) - \frac{h_1}{2\cos\alpha}f_1'''(\alpha) + \frac{1}{2}((1-\lambda_2^2)f_2(\alpha) +f_2''(\alpha))= 0. \end{equation}

Together with (2.35) we then get

(2.44)\begin{equation} f_2''(\alpha) = \frac{h_1}{\cos\alpha}(f_1'''(\alpha)-p_1(2\lambda_1-p_1) f_1'(\alpha)). \end{equation}

2.2.6. Normal stress balance (2.17)

Observing that $x(r) = (r^2-h(r)^2)^{1/2}$ and $y=1+h(r)$ and using (2.19), it may be shown that the signed curvature (2.16)

(2.45)\begin{align} \kappa(r) & ={-}\frac{h_1p_1(p_1+1)}{\cos\alpha} r^{p_1-1} -\frac{h_2p_2(p_2+1)}{\cos\alpha} r^{p_2-1}\nonumber\\ &\quad + h_1^2\tan\alpha \sec^2\alpha p_1(2p_1+1)r^{2p_1-1} + \frac{c_p}{\gamma} + {\rm h.o.t.} \end{align}

Let us now integrate with respect to $r$ the $r$-component of the conservation of linear momentum equation (2.2):

(2.46)\begin{equation} 0={-}\frac{\partial p}{\partial r} + \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial v_r}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 v_r}{\partial \theta^2} -\frac{v_r}{r^2} - \frac{2}{r^2}\frac{\partial v_\theta}{\partial\theta}. \end{equation}

Then, supposing that $\lambda _n\not =1$ and in agreement with Michael (Reference Michael1958) up to an additive constant, we get

(2.47)\begin{equation} p = \sum_{n=1}^\infty\frac{r^{\lambda_n-1}}{(\lambda_n-1)}\left((\lambda_n+1)^2f_n'(\theta) + f_n'''(\theta)\right) +p_a - c_p. \end{equation}

Using this, (2.12) and (2.24a,b) we end up with the following expression for the normal stress on $y=1+h(r)$:

(2.48)\begin{align} \|\boldsymbol{n}\|_2^2 \sigma_{nn} & ={-}\left(\sum_{n=1}^\infty\frac{r^{\lambda_n-1}}{(\lambda_n-1)}\left((\lambda_n+1)^2f_n'(\theta) + f_n'''(\theta)\right)+p_a-c_p\right)\overbrace{\left(h'^2+\frac{(r-hh')^2}{r^2-h^2}\right)}^{n_r^2+n_\theta^2} \nonumber\\ &\quad + 2\sum_{n=1}^\infty \lambda_nr^{\lambda_n-1}f_n'(\theta)\underbrace{\left(h'^2\cos{2\theta} -2\frac{h'(r-hh')}{\sqrt{r^2-h^2}}\sin{2\theta} - \frac{(r-hh')^2}{{r^2-h^2}}\cos{2\theta} \right)}_{n_r^2-n_\theta^2}\nonumber\\ &\quad + \sum_{n=1}^\infty\left((1-\lambda_n^2)r^{\lambda_n-1} f_n(\theta) + r^{\lambda_n-1}f_n''(\theta)\right)\nonumber\\ &\quad \times \underbrace{\left({-}h'^2\sin{2\theta} -2\frac{h'(r-hh')}{\sqrt{r^2-h^2}}\cos{2\theta} + \frac{(r-hh')^2}{{r^2-h^2}}\sin{2\theta} \right)}_{n_rn_\theta}, \end{align}

where

(2.49)\begin{equation} \|\boldsymbol{n}\|_2^2 =n_r^2+n_\theta^2 = h'^2+\frac{(r-hh')^2}{r^2-h^2} = 1+O(r^{2p_1}) \end{equation}

is the square of the (Euclidean) norm of a normal vector $\boldsymbol {n}=(n_r,n_\theta )$ to the free surface.

2.2.7. Determination of $h_1$, $h_2$ and $f_2(\theta )$

Determination of $h_1$. With the choice $p_1=\lambda _1$ (that is, assuming that curvature effects enter into a dominant balance with the normal stress on the free surface) and using the developments (2.28)–(2.30) we can equate the dominant $O(r^{\lambda _1-1})$ terms in (2.18) from (2.45) and (2.48), to obtain

(2.50) \begin{align} -\frac{\gamma h_1\lambda_1(\lambda_1+1)}{\cos\alpha} & ={-}\frac{1}{(\lambda_1-1)}\left((\lambda_1+1)^2f_1'(\alpha) + f_1'''(\alpha)\right) - 2\lambda_1f_1'(\alpha),\notag\\ \Rightarrow \frac{\gamma h_1}{\cos\alpha} & = \frac{(3\lambda_1^2+1)f_1'(\alpha) + f_1'''(\alpha)}{\lambda_1(\lambda_1^2-1)}. \end{align}

From (2.38)–(2.40), for example, the coefficients $A_1$, $C_1$ and $D_1$ appearing in the definition (2.22) of $f_1$ may be re-expressed in terms of $B_1$. The subsequent evaluation of $f_1'(\alpha )$ and $f_1'''(\alpha )$, the use of the eigenvalue equation (2.42) and some simplification then allows us to rewrite (2.50) as

(2.51)\begin{equation} \gamma h_1= \frac{4B_1\lambda_1\sin\alpha\cos\alpha}{(1-\lambda_1)\sin\lambda_1\alpha} = \frac{4B_1\cos\lambda_1\alpha}{(1-\lambda_1)}. \end{equation}

This bears some resemblance to (A.33) of Ramalingam (Reference Ramalingam1994), but note that his $B_1$ corresponds to our $A_1$. Using the numerical scheme to be presented in the next section we have chosen to calculate $C_1$ rather than $B_1$ and in terms of this coefficient, $h_1$ is written as

(2.52)\begin{equation} \gamma h_1 ={-}\frac{4C_1\cos\lambda_1\alpha\sin\alpha}{(\cos\alpha-\lambda_1\cot\lambda_1\alpha\sin\alpha)}. \end{equation}

Determination of $h_2$. With the choices $p_1=\lambda _1$, $p_1=\lambda _2-\lambda _1$ and $p_2=\lambda _2$ and using the developments (2.28)–(2.30) and the results (2.34) and (2.37) we can also equate the $O(r^{\lambda _2-1})$ terms in (2.18) from (2.45) and (2.48), to obtain

(2.53) \begin{align} & -\frac{1}{(\lambda_2-1)}((3\lambda_2^2+1)f_2'(\alpha) + f_2'''(\alpha)) + \frac{h_1}{(\lambda_1-1)\cos\alpha}f_1^{(iv)}(\alpha) \notag\\ & \quad ={-}\frac{\gamma h_2\lambda_2(\lambda_2+1)}{\cos\alpha} + \gamma h_1^2 \tan\alpha\sec^2\alpha \lambda_1(2\lambda_1+1), \notag\\ \Rightarrow & \frac{\gamma h_2}{\cos\alpha} = \frac{(3\lambda_2^2+1)f_2'(\alpha) + f_2'''(\alpha)}{\lambda_2(\lambda_2^2-1)} + \frac{\gamma h_1^2}{2}\tan\alpha\sec^2\alpha. \end{align}

In deriving (2.53) we have made use of the result that $f_1^{(iv)}(\alpha )=0$; something that follows immediately from the definition of $f_1$ from (2.22) with coefficients $A_1$ to $D_1$ and $\lambda _1$ and $\alpha$ satisfying (2.38)–(2.42).

Determination of $f_2(\theta )$. Referring to the linear combination (2.22) in the case $n=2$ we see that (2.35), (2.38), (2.39) and (2.44) represent a system of four equations in the four unknowns $A_2$, $B_2$, $C_2$ and $D_2$, which may therefore be found in terms of $h_1$, $\alpha$ and one of $A_1$, $B_1$, $C_1$ or $D_1$. The coefficient matrix for this system is non-singular since $p_1=\lambda _1=\lambda _2-\lambda _1 \Rightarrow \lambda _2=2\lambda _1$ (in agreement with the $r$ exponent in the correction $Ca\psi _1$ to the leading-order term in the perturbation expansion (2.13a) of the stream function in Anderson & Davis (Reference Anderson and Davis1993)) and therefore, unless $\alpha =3{\rm \pi} /2$, does not satisfy (2.42).

The assumptions that $p_1=\lambda _1$ and $p_1=\lambda _2-\lambda _1$ (and hence that $\lambda _2=2\lambda _1$) are necessary in order that the separation angle $\alpha$ be different from ${\rm \pi}$. This is proved in Lemma A.1 (see Appendix A).

In summary of the results of the analysis of § 2.2, the singular expansion (2.21) for the stream function may now be written as

(2.54)\begin{equation} \psi = 1 + r^{1+\lambda_1}f_1(\theta) + r^{1+2\lambda_1}f_2(\theta) + \ldots, \end{equation}

where the terms on the right-hand side depend upon $r$, $\lambda _1$, $\theta$ and one of $A_1$, $B_1$, $C_1$ or $D_1$. We will find $C_1$ as part of the solution of the boundary integral discretization of the system of (2.5a,b) and in § 3.1.1 explain precisely how. Ramalingam (Reference Ramalingam1994) correctly stated that $\lambda _1$ and $\alpha$ cannot be decided by a local analysis but must be determined by matching with the far-field solution. The same is true of the pressure constant $c_p$. However, he did not perform such a matching. The coefficients $h_1$ and $h_2$ may be determined from (2.50) and (2.53), respectively. Details of the system of nonlinear equations to be solved for $\lambda _1$, $\alpha$, $c_p$, $h_1$ and $h_2$ will be described in § 3.2.2.

3.1.1. Discretization

A boundary element method for the determination of $(\psi,\omega )$ in the domain $\varOmega$ with a given free surface $\varGamma$ will consist of a quadrature evaluation of (3.3)–(3.4) and incorporate the boundary conditions given in § 2.1.1 and free surface conditions (2.11) and (2.15). Referring to figure 3 we divide $BC'$ into $N_1$ line segments. Unlike Kelmanson (Reference Kelmanson1983b), the singular solution

(3.7)\begin{equation} \psi(r, \theta) \sim 1 + r^{1+\lambda_1}f_1(\theta) + r^{1+2\lambda_1}f_2(\theta), \end{equation}

(see (2.21)) and its derivatives are used along $C'C''$.

(N.B. An alternative to this approach, allowing us to take $r_c\rightarrow 0$, would be to use a reformulation of the boundary integral equation. We note that in Stokes flow the tangential derivative of the pressure $p$ is the same as the normal derivative of $\omega$ along a curve,

(3.8)\begin{equation} \frac{\partial p}{\partial s} = \frac{\partial\omega}{\partial n}, \end{equation}

so that, following Hansen & Kelmanson (Reference Hansen and Kelmanson1994), the first terms in the integrals on the right-hand sides of (3.3) and (3.4) may be replaced by $({\partial p}/{\partial s})G_2$ and $({\partial p}/{\partial s})G_1$, respectively. An integration by parts may now be performed along the domain boundary, the integrability of $-p({\partial G_i}/{\partial s})$ $(i=1,2)$ at the separation point allowing us to shrink the arc radius to zero. This merits further investigation.)

Again, unlike Kelmanson (Reference Kelmanson1983b) $C''D$ is divided into ($N_2$) arcs rather than chords. Finally, $EA$ is divided into $N_3$ segments. We denote the total number of arcs and line segments by $M:=N_1+N_2+N_3$ and, ordering these in a clockwise direction, beginning with the segment having a left-hand end point at $B$ we have

(3.9)\begin{equation} BC'\cup C''D\cup EA = \bigcup_{i=1}^M \partial\varOmega_i, \end{equation}

where the $i$th arc or line segment is denoted by $\partial \varOmega _i$. The abscissa of the $i$th point $(x_i',y_i')$ of (3.3)–(3.4) is chosen to be that of the midpoint of the segment $\partial \varOmega _i$ in $BC'$ or $EA$ or that of the midpoint of the chord corresponding to the arc $\partial \varOmega _i$ in $C''D$. The corresponding ordinate $y_i'$ is chosen so that $(x_i',y_i')\in \partial \varOmega _i$. Then, one last point $(x_{M+1}',y_{M+1}') \in \partial \varOmega$, different from the $M$ others, is chosen. The $2M+1$ unknowns to be determined from the discretization of (3.3)–(3.4) are then numerical approximations to:

1. (i) $\omega (x_i',y_i')$ and $({\partial \omega }/{\partial n})(x_i',y_i')$, for $i=1,2,\ldots, N_1$ (along $BC'$);

2. (ii) $({\partial \psi }/{\partial n})(x_i',y_i')$ and $({\partial \omega }/{\partial n})(x_i',y_i')$, for $i=N_1+1,\ldots, M$ (along $C''D$ and $EA$); and

3. (iii) one of the coefficients, ($C_1$, say), appearing in the function $f_1(\theta )$ of (3.7) (the other coefficients being found from (2.38)–(2.40)).

Note that along $C''D$, $\omega$, appearing in the integrals on the right-hand sides of (3.3)–(3.4), is replaced by $\displaystyle -2\kappa ({\partial \psi }/{\partial n})$, as explained in (2.15).

We denote the numerical approximations to the variable evaluations in (i) and (ii) in the list above by $\omega _i$, $\omega _{n,i}$ and $\psi _{n,i}$. The $j$th equations in a discretized form of (3.3) and (3.4), leading to a full rank linear system of $2M+1$ equations, may now be written as

(3.10)\begin{align} &\frac{\xi(x_j',y_j')}{2{\rm \pi}}\psi(x_j',y_j') \nonumber\\ &\quad = \sum_{i=1}^{N_1}\left(\omega_{n,i}\int_{\partial\varOmega_i}G_2(x,1,x_j',y_j')\,{{\rm d} x} -\omega_i\int_{\partial\varOmega_i}\frac{\partial G_2}{\partial y}(x,1,x_j',y_j')\,{{\rm d}x}-\frac{\partial G_1}{\partial y}(x,1,x_j',y_j')\,{{\rm d} x}\right) \nonumber\\ &\qquad + \sum_{i=N_1+1}^{M}\left(\omega_{n,i}\int_{\partial\varOmega_i}G_2(x,y,x_j',y_j')\,{\rm d}s +\psi_{n,i}\int_{\partial\varOmega_i}\left(2\kappa\frac{\partial G_2}{\partial n}(x,y,x_j',y_j')+G_1(x,y,x_j',y_j')\right)\,{\rm d}s\right.\nonumber\\ &\qquad \left.-\psi_i\int_{\partial\varOmega_i}\frac{\partial G_1}{\partial n}(x,y,x_j',y_j')\,{\rm d}s\right) \nonumber\\ &\qquad + \int_{C'C''\cup AB\cup DE} \left(\frac{\partial\omega}{\partial n}G_2 - \omega\frac{\partial G_2}{\partial n} +\frac{\partial\psi}{\partial n}G_1 -\psi\frac{\partial G_1}{\partial n}\right)\,{\rm d}s,\quad j=1,2,3,\ldots,M+1, \end{align}

(3.11)\begin{align} \frac{\xi(x_j',y_j')}{2{\rm \pi}}\omega_j & = \sum_{i=1}^{N_1}\left(\omega_{n,i}\int_{\partial\varOmega_i}G_1(x,1,x_j',y_j')\,{{\rm d} x} -\omega_i\int_{\partial\varOmega_i}\frac{\partial G_1}{\partial y}(x,1,x_j',y_j')\,{{\rm d} x}\right) \nonumber\\ &\quad + \sum_{i=N_1+1}^{M}\left(\omega_{n,i}\int_{\partial\varOmega_i}G_1(x,y,x_j',y_j')\,{\rm d}s +2\psi_{n,i}\int_{\partial\varOmega_i}\kappa\frac{\partial G_1}{\partial n}(x,y,x_j',y_j')\,{\rm d}s\right)\nonumber\\ &\quad + \int_{C'C''\cup AB\cup DE} \left(\frac{\partial\omega}{\partial n}G_1 - \omega\frac{\partial G_1}{\partial n}\right)\,{\rm d}s,\quad j=1,2,3,\ldots,N_1, \end{align}

(3.12)\begin{align} -\frac{\xi(x_j',y_j')}{\rm \pi}\kappa(x_j',y_j') & = \sum_{i=1}^{N_1}\left(\omega_{n,i}\int_{\partial\varOmega_i}G_1(x,1,x_j',y_j')\,{{\rm d} x} -\omega_i\int_{\partial\varOmega_i}\frac{\partial G_1}{\partial y}(x,1,x_j',y_j')\,{{\rm d} x}\right) \nonumber\\ &\quad + \sum_{i=N_1+1}^{M}\left(\omega_{n,i}\int_{\partial\varOmega_i}G_1(x,y,x_j',y_j')\,{\rm d}s +2\psi_{n,i}\int_{\partial\varOmega_i}\kappa\frac{\partial G_1}{\partial n}(x,y,x_j',y_j')\,{\rm d}s\right)\nonumber\\ &\quad + \int_{C'C''\cup AB\cup DE} \left(\frac{\partial\omega}{\partial n}G_1 - \omega\frac{\partial G_1}{\partial n}\right)\,{\rm d}s,\quad j=N_1+1,\ldots,M. \end{align}

3.1.2. Remarks

1. (i) For the sake of simplicity, in our presentation (3.10)–(3.12) have been written out in a form that is fuller than necessary. It should be clear that a number of considerable simplifications can be made: $\psi =0$ along $AE$, for example, and $\kappa =0$ along $BC'$ and $AE$.

2. (ii) All integrals over elements $\partial \varOmega _i$ have either been evaluated exactly or, as in the case of elements along $\varGamma$, using an adaptive quadrature method (Shampine Reference Shampine2008).

3. (iii) Integrals over $AB$ were calculated exactly using the known inflow conditions (see § 2.1.1).

4. (iv) Now consider the right-hand sides of (3.10)–(3.12). At the point $D$ there is a singularity in $\partial \psi /\partial n$ and this can lead to unphysical oscillations in the calculated values of the $x$-component of velocity in the elements on the free surface nearest to $D$. Therefore, rather than evaluate the integrals

   (3.13)\begin{equation} \int_{DE} \left(\frac{\partial\omega}{\partial n}G_2 - \omega\frac{\partial G_2}{\partial n} +\frac{\partial\psi}{\partial n}G_1 -\psi\frac{\partial G_1}{\partial n}\right)\,{\rm d}s \end{equation}
   and
   (3.14)\begin{equation} \int_{DE}\left(\frac{\partial\omega}{\partial n}G_1-\omega\frac{\partial G_1}{\partial n}\right)\,{\rm d}s, \end{equation}
   using the known outflow conditions (again, see § 2.1.1), the numerical solution behaviour near outflow was seen to improve considerably by replacing the integrals with ones over $DD'\cup D'E'\cup E'E$ where $D'$ is the point $(\hat {x}_\infty,\eta _\infty )$ and $E'$ the point $(\hat {x}_\infty,0)$ for some suitable $\hat {x}_\infty > x_\infty$. This is justified on the grounds that given any region $\mathcal {E}\subset \mathbb {R}^2$ (for example, the rectangle $(x_\infty, \hat {x}_\infty )\times (0,\eta _\infty )$) bounded by a positively oriented, piecewise smooth, simple closed curve $\mathcal {S}$ and point $(x',y')\notin \bar {\mathcal {E}}$
   (3.15)\begin{equation} \oint_\mathcal{S}\left(\frac{\partial\omega}{\partial n}G_2 - \omega\frac{\partial G_2}{\partial n} +\frac{\partial\psi}{\partial n}G_1 -\psi\frac{\partial G_1}{\partial n}\right)\,{\rm d}s = \oint_\mathcal{S}\left(\frac{\partial\omega}{\partial n}G_1-\omega\frac{\partial G_1}{\partial n}\right)\,{\rm d}s = 0. \end{equation}
   Over $DD'$, since
   (3.16a,b)\begin{equation} \frac{\partial^2\psi}{\partial n\partial s}(x,\eta(x))\rightarrow 0, \quad \frac{\partial^2\omega}{\partial n\partial s}(x,\eta(x)) \rightarrow 0, \end{equation}
   as $x\rightarrow \infty$, the values of ${\partial \psi }/{\partial n}$ and ${\partial \omega }/{\partial n}$ were assumed to be the same as in $\partial \varOmega _{N_1+N_2}$ and thus the integral
   (3.17)\begin{equation} -\int_{DD'}G_1 \,{\rm d}s, \end{equation}
   appeared in the coefficient matrix when multiplying either $({\partial \psi }/{\partial n})(x',y')$ or $({\partial \omega }/{\partial n})(x',y')$ and $(x',y')\in \partial \varOmega _{N_1+N_2}$. All other contributions to the integral over $DD'\cup D'E'\cup E'E$ were calculated exactly using the known fully developed solution. For all computational results presented in § 4 it was found to be adequate to choose $\hat {x}_\infty = 2{x}_\infty$. In doing this we effectively extended the downstream flow domain to twice its original length.

5. (v) The integrals over the arc $C'C''$ in (3.10)–(3.12) create the elements of the coefficient matrix that multiply $C_1$ in the function $f_1(\theta )$ of (3.7).

6. (vi) A well known disadvantage of boundary element methods is that the coefficient matrix in the resulting system of equations is typically ill-conditioned and non-symmetric. Although there has been recent progress in developing preconditioned GMRES (generalized minimal residual) iterative methods for systems arising in the context of boundary element methods (see, for example, Benedetti, Aliabadi & Davì (Reference Benedetti, Aliabadi and Davì2008)), the small scale of the present problem allowed us to use a direct method such as Gaussian elimination with pivoting.

3.2.1. Free surface $\{(r,\theta ): r>r_c\}$

To determine numerically the free surface shape we use a Levenberg–Marquardt method (Levenberg Reference Levenberg1944; Marquardt Reference Marquardt1963) to solve the least-squares problem,

(3.18)\begin{equation} \min_{\boldsymbol{c}} \sum_{i=1}^{N+1} R_i(\boldsymbol{c})^2, \end{equation}

where $R_i$ ($i=1,2,3,\ldots, N+1$) is the residual of the normal stress balance (see (2.18))

(3.19)\begin{equation} R_i:=[p(x_i,\eta(x_i))]+2\frac{\partial^2\psi}{\partial n\partial s}(x_i,\eta(x_i)) + \gamma\kappa(x_i,\eta(x_i)), \end{equation}

evaluated at a point $(x_i,\eta (x_i))$ on the free surface (for some $x_i\geq x_0$) and depending on some vector of parameters $\boldsymbol {c}$.

Calculation of $R_i$. In order to calculate $R_i$ we need both the pressure excess over atmospheric pressure, $[p]$, and the tangential derivative of the tangential free surface velocity $\hat {\boldsymbol {t}}\boldsymbol {\cdot }\boldsymbol {v}(={\partial \psi }/{\partial {n}})$. From (2.2) we get

(3.20)\begin{equation} \hat{\boldsymbol{t}}\boldsymbol{\cdot}\boldsymbol{\nabla} [p] = \hat{\boldsymbol{t}}\boldsymbol{\cdot}\Delta\boldsymbol{v} \Rightarrow \frac{\partial p}{\partial s} = \frac{\partial\omega}{\partial n} \end{equation}

along $\varGamma$. Let us define the arc $\partial \varOmega _i$ ($i=N_1+1,\ldots,N_1+N_2$) along $\varGamma$ by

(3.21)\begin{equation} \partial\varOmega_i =\left\{(x,\eta(x)): a_i\leq x\leq a_{i+1}\right\}, \end{equation}

so that $\partial \varOmega _i$ has end points $(a_i,\eta (a_i))$ and $(a_{i+1},\eta (a_{i+1}))$, for a certain choice of $\{a_i\}$ such that

(3.22)\begin{equation} x_0=a_{N_1+1}< a_{N_1+2} < \ldots < a_{N_1+N_2+1} = x_\infty. \end{equation}

Then, if the outflow excess pressure $[p(D)]=[p(a_{N_1+N_2+1},\eta (a_{N_1+N_2+1}))] = 0$ we calculate

(3.23)\begin{equation} [p(a_i)] = [p(a_{i+1})]-\omega_{n,i}\int_{\partial\varOmega_i}\,{\rm d}s,\quad i=N_1+N_2,\ldots,N_1+1. \end{equation}

At $D$ we set the tangential derivative of the tangential free surface velocity equal to zero and to calculate

(3.24)\begin{equation} \frac{\partial^2\psi}{\partial n\partial s}(a_i,\eta(a_i)),\quad i=N_1+1,N_1+2,N_1+3,\ldots,N_1+N_2, \end{equation}

we use simple finite difference formulae.

The vector of parameters $\boldsymbol {c}$. The vector $\boldsymbol {c}$ of (3.18) consists in the present paper of the parameters $\{\alpha _i, \beta _i, \gamma _i, \epsilon _{0,i}, \epsilon _{\infty,i}\}$ appearing in a representation of $\eta$ as

(3.25)\begin{equation} \eta(x) = \sum_{i=1}^n \beta_i \eta_i(x), \end{equation}

of $n$ linearly independent functions

(3.26)\begin{equation} \eta_i(x) = 1 + \alpha_i\tanh(x(\epsilon_{\infty,i}+(\epsilon_{0,i}-\epsilon_{\infty,i})\exp(-\gamma_ix))), \end{equation}

(see Kelmanson Reference Kelmanson1983b). Since we require that the linear combination (3.25) be equal to $1$ at $x=0$, we have the constraint that

(3.27)\begin{equation} \sum_{i=1}^n \beta_i = 1, \end{equation}

so that the number of free parameters in the linear combination (3.25) totals $5n-1$ and thus $\boldsymbol {c}\in \mathbb {R}^{5n-1}$. Having chosen $n$ we then set $N=5n-2$.

3.2.2. Free surface $\displaystyle 1 - r\sin \alpha + h_1r^{\lambda _1+1} + h_2r^{\lambda _2+1} - \frac {c_pr^2}{2\gamma }\cos \alpha$ with $r\leq r_c$

Coupled with the numerical determination of the free surface $\eta (x)$ away from the singularity, as described above, we have, from (2.19), (2.42), (2.45), (2.50), (2.52) and (2.53) the problem of finding the solution $(h_1, h_2, \alpha, \lambda _1, c_p)$ to the nonlinear system,

(3.28)\begin{gather} 1 - r_c\sin\alpha + h_1r_c^{\lambda_1+1} + h_2r_c^{2\lambda_1+1}-\frac{c_pr_c^2}{2\gamma}\cos\alpha = \eta(x_0), \end{gather}

(3.29)\begin{gather} -\frac{h_1\lambda_1(\lambda_1+1)}{\cos\alpha} r^{\lambda_1-1} -\frac{2h_2\lambda_1(2\lambda_1+1)}{\cos\alpha} r^{2\lambda_1-1}\nonumber\\ + h_1^2 \lambda_1(2\lambda_1+1)(\tan\alpha \sec^2\alpha) r^{2\lambda_1-1} + \frac{c_p}{\gamma} = \kappa(x_0),\end{gather}

(3.30)\begin{gather}\frac{\gamma h_1}{\cos\alpha} = \frac{(3\lambda_1^2+1)f_1'(\alpha) + f_1'''(\alpha)}{\lambda_1(\lambda_1^2-1)} ={-}\frac{4C_1\cos\lambda_1\alpha\tan\alpha}{\cos\alpha-\lambda_1\cot\lambda_1\alpha\sin\alpha}, \end{gather}

(3.31)\begin{gather}\frac{\gamma h_2}{\cos\alpha} = \frac{(12\lambda_1^2+1)f_2'(\alpha) + f_2'''(\alpha)}{2\lambda_1(4\lambda_1^2-1)} + \frac{\gamma h_1^2}{2}\tan\alpha\sec^2\alpha, \end{gather}

(3.32)\begin{gather}\lambda_1\sin{2\alpha} - \sin{2\alpha\lambda_1} =0, \end{gather}

where, on the right-hand side of (3.29),

(3.33)\begin{equation} \kappa(x_0):= \frac{\eta''(x_0)}{(1+\eta'^2(x_0))^{3/2}}, \end{equation}

is the signed curvature of the surface $y=\eta (x)$ at $x=x_0$. The solution of the system (3.28)–(3.32) is achieved using the trust-region dogleg algorithm fsolve in MATLAB, see Powell (Reference Powell1970).