The separation angle of the free surface of a viscous fluid at a straight edge

Abstract

We rework some of the die-swell singularity analysis for Stokes flow, originally by Ramalingam (Ramalingam, 1994 Fiber spinning and rheology of liquid-crystalline polymers, PhD thesis, Massachusetts Institute of Technology) in appendix A of his PhD thesis, in an attempt to demonstrate that for capillary numbers in the range $(0,\infty )$ the curvature may enter into the normal stress balance on the free surface and lead to separation angles exceeding $180^{\circ }$ and infinite curvature at the separation point. The singular coefficients in the asymptotic solution and the free surface shape in a neighbourhood of the separation point cannot be determined by a local analysis of the Michael type (Michael, Mathematika, vol. 5, 1958, pp. 82–84) but must be found from matching with the solution valid away from the die edge. The numerical method that we use in the truncated die-swell domain is a boundary element method incorporating the singular solution near the separation point. Although there is some variation in the extrudate swell ratios at different capillary numbers reported in the numerical literature, our results for capillary numbers $Ca$ from $1$ to $1000$ are within the range of values published in earlier papers. The computed separation angles at different values of $Ca$ agree well with the range of separation angles to be found in experimental and numerical papers. The separation angle appears to converge to a value different from $180^{\circ }$ as $Ca$ increases, leading us to conclude that the case of zero surface tension ($Ca=\infty$, with corresponding separation angle of $180^{\circ }$), is a singular limit.

1. Introduction

The extrusion of a fluid, Newtonian or non-Newtonian, from a die in the form of a right-circular tube or a straight channel, say, may give rise to a jet whose free surface is different in shape from that of the die orifice. In the case of a viscoelastic fluid the extrudate may have a diameter several times that of the die and gives rise to a phenomenon known as die swell or extrudate swell, the correct prediction of which is of crucial importance in maintaining the quality of processing applications of polymer melts such as blow moulding or fibre spinning (Koopmans 1999). Even in the case of Newtonian fluids where polymer chains, giving rise to hoop and extra normal stresses, are absent the jet will typically have non-zero Gaussian curvature, the cross-sectional diameter of the jet at any distance from the orifice being dependent on factors including the Reynolds number ($Re$) and capillary number ($Ca$: the ratio between viscous drag forces and surface tension forces acting across the free surface). The expansion or contraction of jets of Newtonian fluids under different flow conditions and fluid properties seems to have first been reported experimentally by Middleman & Gavis (1961).

One of the simplest extrusion flows of practical interest is that of the creeping planar flow of a Newtonian fluid from a straight channel into air, in the absence of gravity. Although there is now broad agreement in the literature, both experimental and computational, on the final extrudate swell ratios of the jet at different capillary numbers, the study of the asymptotic solution of Stokes flow near the separation point remains one of great importance, for at least two reasons.

First, this simple flow is relevant to more complex cases.

1. (i) The neglect of inertial terms remains valid even near the separation point since, if the components of the velocity $\boldsymbol {v}$ vary radially like $r^\lambda$ for some $\lambda$ with $Re(\lambda )>0$ (see (2.24a,b)), then the Cauchy stress $\sim r^{\lambda -1}$ and dominates over $\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {v} \sim r^{2\lambda -1}$.

2. (ii) As pointed out by Ramalingam (1994), the analysis of the singular behaviour of the flow variables near the orifice of an axisymmetric flow may be expected to appear locally planar. Therefore, the asymptotic behaviour in the present case should not be expected to be significantly different from the axisymmetric case.

3. (iii) The neglect of gravity may be justified following the argument of, for example, Michael (1958): gravity gives variations in the free surface stresses of $O(r)$ and this, of course, corresponds to a contribution of $O(r^{\lambda -1})$ with $\lambda =2$. However, as will be shown later, the strongest modes have exponents $\lambda < 2$, so that gravity plays no part in forming the free surface very close to the point of separation.

4. (iv) The asymptotic form of the velocity field and even the stresses of this simple flow are of great relevance to some viscoelastic fluid models and this is of interest because of the polluting effects that stress singularities can have on the computation of extrusion flows of polymer solutions and melts. For, example, Evans, Palhares Junior & Oishi (2017) and Evans & Evans (2019) have proved that for extrudate flow of both the Phan-Thien–Tanner (PTT) and Giesekus models in the presence of a solvent viscosity, the velocity field is dominated by the Newtonian contribution near the join of the die wall and free surface. A modified upper-convected Maxwell (MUCM) model was derived by Apelian, Armstrong & Brown (1988) from network theory and the fluid relaxation time made to decrease at increasing values of the trace of the stress tensor. A consequence of this was that the viscoelastic stress fields reduced to the asymptotic expressions for a Newtonian fluid near singularities at non-deformable boundaries, including therefore those in ‘stick-slip’ flow. It should be made clear, however, that viscoelastic stresses can have a dramatic effect on the shape of the free surface in extrusion flows: see, as one example of many articles that could be cited, the recent paper of Varchanis et al. (2020). Since the angle of separation is determined by matching the local singular solution with the far-field solution this too is expected to be different from the Newtonian value, even if the Reynolds number and the surface tension are the same.

The second reason for the ongoing importance of studying the asymptotic solution of Stokes flow near the separation point is that an apparent discrepancy for this problem between theoretical results on the one hand and computational and experimental results on the other persists to the present day. Michael (1958), in his landmark paper of more than 60 years ago, demonstrated using a local analysis and under the assumption of either a planar free surface near the separation point or zero surface tension, that the angle of separation $\alpha$ between the free surface and the channel wall (see figure 1) is $180^{\circ }$. However, there has been a wealth of results in the literature in the decades since then which seem to contradict his theoretical result. Indeed, going back as far as some of the earliest finite element (Nickell, Tanner & Caswell 1974) and boundary element (Kelmanson 1983b; Tanner, Lam & Bush 1985; Tanner 1986) calculations and experimental results (Batchelor, Berry & Horsfall 1973; Nickell et al. 1974), agreement on the shape of the free surface was good and consistently showed that at the channel exit the free jet surface forms an angle in excess of $180^{\circ }$. The experimental results reported in Nickell et al. (1974), Tanner (1986) and Tanner et al. (1985) were all at Reynolds numbers below $10^{-3}$ and those of Batchelor et al. (1973) were for $Re\approx 10^{-8}$. The data from all these studies showed separation angles of around $192^{\circ }$ for the computations and between $189^{\circ }$ and $194^{\circ }$ for the experiments.

Figure 1. Problem geometry of planar extrusion.

Schultz & Gervasio (1990) modified the earlier matched eigenfunction method of Trogdon & Joseph (1980, 1981) to study the planar extrusion flow of a Newtonian fluid at zero Reynolds number. Results for the free surface slope at the separation point seemed to show this tended to $0$ when $Ca^{-1}=0$, as the number of eigenfunctions in their expansions was increased sufficiently, consistent with the analysis of Michael (1958). However, for non-zero surface tensions, no evidence of a zero free surface slope was found, and indeed for the larger surface-tension curve $Ca^{-1} = 1$ it was found to be quite clear that the slope was not converging to zero with an increase in the number of eigenfunctions.

More recently, mixed finite element methods with local irregular mesh refinement near the separation point were used in computations by Salamon et al. (1995) of planar extrudate flow at differing capillary numbers and Reynolds numbers, with and without slip along the die wall. Verificatory calculations at $Ca=0$ (the ‘stick-slip’ problem) and $Re=0$ were stated to show excellent agreement for the computed pressure along the wall with the analytical results of Richardson (1970), obtained with the Weiner–Hopf method, and with the finite element results of Nickell et al. (1974), referred to already above. The coefficients in the asymptotic expansion for the $x$-component of the velocity along the free surface calculated by Salamon et al. (1995) were also pronounced to be in excellent agreement with those in the literature (Richardson 1970; Ingham & Kelmanson 1984; Georgiou et al. 1989; Tanner & Huang 1993). It should be noted, however, that the leading coefficient ascribed to Richardson (1970) in table III of Salamon et al. (1995) is in fact that of Tanner & Huang (1993) and that Richardson's leading coefficient is some $16\,\%$ smaller than the correct value. For larger values of $Ca$ the calculations of Salamon et al. (1995) predicted that the slope of the free surface at the die edge was strictly positive and, following an earlier conjecture by Schultz & Gervasio (1990) that the free surface $y=\eta (x)$ near the separation point was of the form

(1.1)\begin{equation} \eta(x)-1 = a + bx + cx^n, \end{equation}

determined that the best-fit coefficients were $a=0$, $b=0.176$, $c=0.0263$ and $n=1.43$ for $Ca=1$ and $Re=0$. This gives an angle of separation of $(180+ \text {arctan}(0.176))\approx 189.98^{\circ }$ and, since $1< n<2$, predicts both infinite curvature and integrable stresses at the separation point, in keeping with the hypothesis of Schultz & Gervasio (1990). The authors stated themselves to be in agreement with the earlier arguments of Ramalingam (1994) that the contact angle $\alpha$ could be determined from the matching of the asymptotic solution (valid close to the separation point) with the bulk flow.

Anderson & Davis (1993) performed a regular perturbation expansion in the limit of small capillary number, similar to that of Richardson (1970) but, unlike his expansion, theirs was valid near the corner singularity. Consistent with the results of Schultz & Gervasio (1990) the authors showed that their local solution predicted a free surface having infinite curvature at the corner which balances the normal force along the free surface.

Despite the weight of evidence presented above for a contact angle different from $180^{\circ }$ for all but the cases $Ca=0,\infty$ it would be disingenuous to suggest that all authors have been of this persuasion. In an attempt to reconcile what were seen to be contradictory results from theory, computation and experiments for the separation angle between the free surface and the die wall, Tanner (1986) and Tanner et al. (1985) used boundary element methods to numerically investigate the effect on the jet shape of rounding the tube exit. By setting the normal stress to zero just upstream of the separation point on the tube exit the authors found that the initial departure of the free surface from the solid, rounded wall was tangential, as demanded by the Michael theory. The authors were not able to explain satisfactorily, however, why computations with sharp edged walls did not conform to the condition of Michael. Trogdon & Joseph (1981) attempted to solve the problem of determining the swell of a low speed planar jet by expressing the solution for the flow as biorthogonal eigenfunction series and expressing the boundary conditions $F_i(x,\eta (x))=0$ ($i=1,2,3$) on the free surface in terms of conditions evaluated on $y=1$ by using Taylor series:

(1.2)\begin{equation} F_i(x,\eta(x)) \approx F_i(x,1) + (\eta(x)-1)\frac{\partial F_i}{\partial y}(x,1) =0. \end{equation}

The interface shape that resulted was a regular perturbation to the flat jet surface and the authors predicted two possible values of the separation angle for any value of $Ca$, including $\alpha =180^{\circ }$. This they chose as the physically more reasonable choice and arrived at the conclusion that the dominant singular behaviour for the stresses was $O(r^{-1/2})$, irrespective of the value of $Ca$. However, the approach of Trogdon & Joseph has come under criticism by both Ramalingam (1994) and Salamon et al. (1995), who point out that the Taylor series (1.2) for the kinematic and stress conditions used by Trogdon & Joseph are not valid near the singularity $(0,1)$ because velocity gradients and the pressure become singular as $x\rightarrow 0$. Sturges (1979) has drawn attention to the fact that the local analysis of Michael (1958) gave the angle of separation between a planar free surface of a Newtonian fluid and a planar solid boundary and that therefore the source of apparent discrepancy between the prediction of the Michael theory ($\alpha =180^{\circ }$) and the results of computations and experiments of the die-swell problem (even with very small but non-zero surface tensions) resides in the fact that the jet surface for the extrusion of a fluid is not planar.

An outline of the present paper is as follows. After a brief reminder of the governing equations for Stokes flow and an introduction to the problem to be solved and the boundary conditions (§ 2.1) we rework some of the die-swell singularity analysis for Stokes flow, originally by Ramalingam (1994) in appendix A of his PhD thesis, in an attempt to demonstrate that for capillary numbers in the range $(0,\infty )$ the curvature enters into the normal stress balance on the free surface and leads to separation angles different from ${\rm \pi}$ and infinite curvature at the separation point. In this case, the singular coefficients and the free surface shape in a neighbourhood of the separation point cannot be determined by a local analysis of the Michael type (Michael 1958) but must be found from matching with the solution valid away from the die edge. These conclusions are entirely consistent with one scenario given by Ramalingam but we differ from him in some of the details of our calculations, and these are presented in § 2.2.

The numerical method that we use for the solution of Stokes flow in the truncated die-swell domain is a boundary element method incorporating the singular solution near the separation point. Costabel (1987), Ingham & Kelmanson (1984, 1986), Katsikadelis (2016), Kelmanson (1983a,b), Sauter & Schwab (2011) and a huge number of other authors have already described in some considerable detail the derivation, advantages and disadvantages of boundary element methods for the solution of a wide range of boundary value problems in physics and engineering (including ones with boundary singularities) and we accordingly limit our description of the basic method in § 3.1.1 to a brief treatment, in order to give proper emphasis to more novel aspects of our work. Kelmanson (1983a) used a boundary element method to solve the Stokesian ‘stick-slip’ problem. Because the separation angle in this problem is fixed equal to $180^{\circ }$ and the exponents $\{\lambda _n\}_{n=1}^\infty$ in the local asymptotic solution for the stream function $\psi$ are known (cf. (2.21)), Kelmanson was able to rewrite the governing biharmonic equation for $\psi$ as one for $\chi$ where

(1.3)\begin{equation} \chi = \psi - g, \end{equation}

and $g$ was a truncated asymptotic solution chosen so that $\chi$ had no singularities up to fourth derivatives. This allowed the boundary used in the computations to pass right through the separation point, since all terms appearing in the boundary integral formulation for the problem (even the normal derivative of the Laplacian of $\chi$) remained bounded in a neighbourhood of the separation point. The same approach was not possible, however, when the author again used boundary element methods to solve the Stokesian die-swell problem (Kelmanson 1983b) because neither the separation angle nor (therefore) the exponents in the asymptotic solution were known a priori. This led the author to admit that he was unable to take account of the true nature of the solution near the point of separation and that therefore in this region the numerical results should be expected to be in error.

In the context of finite element methods, Georgiou et al. (1989) used special singular finite elements around the separation point in the Newtonian planar ‘stick-slip’ problem. In these elements the field shape functions embodied the known form of the singularity in the radial direction and were compatible with the adjacent ordinary quadrilateral elements. The singular elements for the pressure used lower-order representations than those in the velocity singular elements and because the pressure is singular at the separation point itself no node was placed there. Compared with ordinary finite elements, the singular elements gave more accurate results for relatively coarse meshes and normal stress oscillations on the free surface were practically eliminated. However, as with Kelmanson (1983b), an analogous approach for the die-swell problem, using the correct asymptotic solution near the separation point, was out of reach since the angles of separation and the associated spectra of exponents were not available. Therefore, Georgiou, Schultz & Olson (1990) used singular elements around the separation point having field shape functions in the radial direction of the same form as in the ‘stick-slip’ case. Despite this approximation, the authors found that the use of singular elements speeded up the convergence of the free surface considerably compared with the cases where ordinary finite elements were used. More recently, Georgiou & Boudouvis (1999) used singular finite element methods to solve both the axisymmetric and planar Newtonian extrudate-swell problems. The authors found that their singular finite element method performed better than a standard finite element method for flows at low Reynolds numbers and with small to moderate surface tensions. We compare the results of our simulations with theirs in § 4 of the present paper.

In this paper we show that it is possible to calculate the leading exponent in the local asymptotic solution for $\psi$ as the solution to a transcendental equation involving the separation angle. We use the leading-order asymptotic solution for all variables on an arc of radius $0< r_c\ll 1$ centred at the separation point, the leading singularity expansion coefficient in the asymptotic solution being one of the unknowns to be determined. This is explained in greater detail in § 3.1.

By requiring that the normal component of velocity and the shear stress should both vanish on the free surface, and by imposing suitable conditions elsewhere on the domain boundary (see §§ 2.1.1 and 2.1.2) we are led to a well-posed problem. The normal stress balance (see (2.18)) is then used to determine the correct position of the free surface itself. A Levenberg–Marquardt method, described by Levenberg (1944) and Marquardt (1963) and implemented in lsqnonlin of MATLAB, is used to determine numerically an optimum free surface shape by minimizing the sum of squares of the normal surface stress residual evaluated at certain points along the jet. More details are given in § 3.2.1. Once the jet surface has been determined numerically away from the separation point, as described above, we then proceed to find the parameters in the equation of the free surface near the separation point by matching the two surfaces and their curvatures at a distance $r_c$ from the separation point using the trust-region dogleg algorithm (implemented in fsolve in MATLAB: see Powell (1970)). This is explained in greater detail in § 3.2.2.

2. Stokes flow

The steady inertialess flow of a Newtonian fluid in an $n$-dimensional domain $\varOmega$ may be described by the non-dimensionalized equations of continuity and linear momentum as follows:

(2.1)\begin{align} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v} & = 0, \quad \text{in}\ \varOmega, \end{align}

(2.2)\begin{align} -\boldsymbol{\nabla}{p} + \Delta\boldsymbol{v} & = \boldsymbol{0}, \quad \text{in } \varOmega, \end{align}

where $\boldsymbol {v}$ denotes the fluid velocity and $p$ is the pressure.

From this point onwards we confine our attention to domains in the Cartesian plane. With $n=2$, (2.1) may be satisfied by the introduction of a twice continuously differentiable stream function $\psi =\psi (x,y)$, where $x$ and $y$ are Cartesian coordinates, and the Cartesian components $u$ and $v$ of the velocity are written as

(2.3a,b)\begin{equation} u=\frac{\partial\psi}{\partial y},\quad v={-}\frac{\partial\psi}{\partial x}. \end{equation}

By calculating the curl of (2.2) and thus eliminating the pressure, we arrive at the biharmonic equation for the stream function

(2.4)\begin{equation} \Delta^2\psi = 0. \end{equation}

Introducing the (negative scalar) vorticity $\omega :=\Delta \psi$ we see from (2.4) that $\omega$ satisfies Laplace's equation and that we may therefore rewrite (2.1)–(2.2) in terms of $\psi$ and $\omega$ as

(2.5a,b)\begin{equation} \Delta\omega=0,\quad \Delta\psi=\omega. \end{equation}

2.1. Planar extrusion: geometry and boundary conditions

We now consider the creeping planar extrusion flow of a Newtonian fluid, as illustrated in figure 1. A viscous fluid is driven by an imposed pressure difference, in the absence of gravity, in the channel/die between two parallel solid walls, $x\in (-\infty,0],\; y=\pm {1}$, and is then extruded at $x=0$. Here $(0,1)$ is the separation point and, as will be seen later, where the fluid stress is singular. The free surface of the fluid jet $y=\eta (x)$ is labelled $\varGamma$ in the figure. Sufficiently far upstream the fluid flow is Poiseuille flow and in the jet downstream the flow becomes uniform. Figure 1, for symmetry reasons, shows only the upper portion of the flow geometry.

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 (1967) or Longuet-Higgins (1953)), 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. Planar extrusion: the solution and the shape of the free surface near the point of separation

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

Michael (1958) 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 (1994) 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 (1994) 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 (1965) 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 (1979), the weak regularity condition $Re(\lambda _n)>0$ may be seen to mean that the velocity vanishes as $r\rightarrow 0$. Following Michael (1958), 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 (1993) or (32) of Lugt & Schwiderski (1965)) 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 (1958), 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 (1994).

Suppose now that $p_1= \lambda _2-\lambda _1$, as assumed by Ramalingam (1994). 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 (1994).

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 (1994) in his (A.30). From (2.22), (2.31), (2.34) and (2.37) we have (cf. Michael (1958), 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 (1958), 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 (1949), Lugt & Schwiderski (1965) and Moffatt (1964) in the context of corner and wedge flows, and, for example, by Sturges (1979) and Ramalingam (1994) in the present context. Huilgol & Tanner (1977) 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 (1993).

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 (1958) 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 (1994), 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 (1993)) 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 (1994) 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. Numerical scheme

3.1. An integral formulation of (2.5a,b) using fundamental solutions

Let $\varOmega$ be the truncated domain whose boundary $\partial \varOmega$ is $AB\cup BC'\cup C'C''\cup C''D\cup DE\cup EA$, as shown in figure 3. The inflow boundary $AB$ is set at a distance $-x_{-\infty }$ upstream and the outflow boundary $DE$ a distance $x_\infty$ downstream of the die exit. Now $C'C''$ is defined to be the arc

(3.1)\begin{equation} \{(r_c,\theta): 0\leq \theta \leq {\rm \pi}+ \arcsin(h(r_c)/r_c)\}, \end{equation}

for some choice of $0< r_c\ll 1$. Thus, $C'$ is the point $(-r_c,1)$ and $C''$ the point $(x_0, \eta (x_0))$ where $x_0$ is the positive root of

(3.2)\begin{equation} x_0^2 + (\eta(x_0)-1)^2 = r_c^2. \end{equation}

Suppose that $(x',y')$ is an arbitrary point on $\partial \varOmega$. Then, using Green's identities it may be shown that $\psi (x',y')$ and $\omega (x',y')$, the solutions to (2.5a,b), may be expressed in integral form as

(3.3)\begin{gather} \frac{\xi(x',y')}{2{\rm \pi}}\psi(x',y') = \int_{\partial\varOmega}\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{gather}

(3.4)\begin{gather}\frac{\xi(x',y')}{2{\rm \pi}}\omega(x',y') = \int_{\partial\varOmega}\left(\frac{\partial\omega}{\partial n}G_1-\omega\frac{\partial G_1}{\partial n}\right)\,{\rm d}s, \end{gather}

where the normal derivative $\partial /\partial n := \hat {\boldsymbol {n}}\boldsymbol {\cdot }\boldsymbol {\nabla }$ and $\xi (x',y')$ is the angle between the left- and right-hand tangents at $(x',y')$. The functions $G_1=G_1(x,y,x',y')$ and $G_2=G_2(x,y,x',y')$ are fundamental solutions of, respectively, the two-dimensional Laplace's equation and biharmonic equation, chosen to be

(3.5)\begin{gather} G_1={-}\frac{1}{4{\rm \pi}}\log((x-x')^2+(y-y')^2), \end{gather}

(3.6)\begin{gather}G_2 ={-}\frac{1}{16{\rm \pi}}(((x-x')^2+(y-y')^2)(\log((x-x')^2+(y-y')^2)-2)). \end{gather}

Figure 3. The integrals on the right-hand sides of (3.3) and (3.4) are evaluated along $AB\cup BC'\cup C'C''\cup C''D\cup DE\cup EA$. Along the arc $C'C''$ the singular solution (3.7) and its derivatives are used.

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 (1983b), 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 (1994), 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 (1983b) $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 2008).

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ì (2008)), the small scale of the present problem allowed us to use a direct method such as Gaussian elimination with pivoting.

3.2. Determination of the free surface $\varGamma$

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

To determine numerically the free surface shape we use a Levenberg–Marquardt method (Levenberg 1944; Marquardt 1963) 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 1983b). 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 (1970).

4. Results

4.1. Parameter values and boundary element meshes

Calculations were performed for values of the dimensionless surface tension $\gamma = 0.001$, $0.01$, $0.1$ and $1$, on meshes varying from $M=624$ ($1249$ unknowns) to $M=1854$ ($3709$ unknowns). The radius $r_c$ of the arc $C'C''$ (see figure 3) was set equal to $10^{-5}$. With a choice of $n=10$ basis functions $\eta _i$ (see (3.26)) the number of points $N+1$ where the normal stress balance residual (3.19) was evaluated was equal to $5n-1 = 49$. For any given value of $\gamma$ and on any mesh, the Levenberg–Marquardt algorithm for the determination of the vector of parameters $\boldsymbol {c}\in \mathbb {R}^{5n-1}$ solving the least-squares problem (3.18) was considered to have converged when $\boldsymbol {c}^{(k+1)}$, the $(k+1)$th iterative value of $\boldsymbol {c}$, satisfied

(4.1)\begin{equation} \|\boldsymbol{c}^{(k+1)}-\boldsymbol{c}^{(k)}\|_2 < s\|1+\boldsymbol{c}^{(k)}\|_2, \end{equation}

the relative step tolerance $s$ being set equal to $10^{-6}$. Our choice of the relative step tolerance $s$ led to values of the objective function (3.18) that were always (sometimes significantly) less than $O(10^{-5})$. Numerical continuation was used in order to obtain solutions on increasingly refined meshes. That is, the converged values of $\boldsymbol {c}$ calculated on a given mesh were used as the starting values for computations on finer meshes.

4.2. Numerical results

4.2.1. Extrudate swell ratio

In table 1 and figure 4 we present the results of calculations of the extrudate swell ratio for values of the dimensionless surface tension $\gamma = 1$, $0.1$, $0.01$ and $0.001$ using meshes varying from $M=1434$ ($2869$ unknowns) to $M=1854$ ($3709$ unknowns).

Table 1. Values of the extrudate swell ratio $\eta _\infty$ for different values of the non-dimensional surface tension $\gamma$, computed using the boundary integral equation method and $M$ boundary elements.

Figure 4. Meshes with $N_1=N_2=N_3$ varying from $208$ to $618$. Extrudate swell ratios $\eta (\infty )$ computed at values $1$, $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

Comparison is made in table 1 with the results of the singular finite element calculations of Georgiou & Boudouvis (1999) and those of the finite element method of Mitsoulis, Georgiou & Kountouriotis (2012), the method of fundamental solutions of Poullikkas et al. (1998) and the high-resolution finite element method of Salamon et al. (1995). The results of Poullikkas et al. (1998), although in reasonable agreement over the range $\gamma \in [0.001, 1]$ with some earlier numerical results in the literature, such as those of Kelmanson (1983b), are rather different from those of the other authors shown in this table. Although there is a small amount of spread in the values obtained by Georgiou & Boudouvis (1999), Mitsoulis et al. (2012) and Salamon et al. (1995), the result of Salamon et al. (1995) for $\gamma =1$ was computed on a highly refined finite element mesh involving $171\,258$ unknowns and may be considered, we believe, to be a benchmark value. The value of the extrudate swell ratio at this value of $\gamma$ computed by Georgiou & Boudouvis (1999) is within $0.026\,\%$ of the benchmark value and was obtained on a finite element mesh leading to $30\,866$ unknowns. Our own best calculation of $\eta (\infty )$ at $\gamma =1$ is within $0.061\,\%$ of the benchmark value but was obtained using a mesh having only $3709$ unknowns. The cost advantages of using a boundary element method are clear.

There is evidence in our results presented in both table 1 and figure 4 of extrudate swell ratios that continue to show small increases as the mesh is refined, so that we would expect the true converged values to be slightly higher at all values of the surface tension $\gamma$ than those shown here and bringing them into even closer agreement with the results of Georgiou & Boudouvis (1999) and Salamon et al. (1995). However, calculations for $M>1854$ are beyond what we can perform conveniently with the computing resources currently at our disposal.

Figure 5 shows, on a log-normal scale, the extrudate swell ratios calculated on the finest mesh and, although with the present code we are unable to compute the free surface at zero surface tension, the graph shows a clear trend to a limiting value of $\eta (\infty )$ as $\gamma \rightarrow 0$.Tanner (1988) presented a number of extrudate swell ratios at zero surface tension available in the literature at that time and, based on these, estimated that $\eta (\infty )$ was in the range $[1.188, 1.192]$. On their finest ordinary finite element mesh Georgiou & Boudouvis (1999) calculated the extrudate swell ratio at zero surface tension to be equal to $1.1869$. The corresponding values obtained by Taliadorou, Georgiou & Mitsoulis (2008) and Claus, Cantwell & Phillips (2015), for example, were $1.1878$ and $1.1891$, respectively. The trend of the results displayed in figure 5 is certainly consistent with these limiting values. As will be explained in § 5.2, when $\gamma =0$ we should have a separation angle equal to ${\rm \pi}$ and a surface that is locally described by $y = 1+O(r^2)$. This is not the case for any of the results presented in figure 5 and we take up this point in more detail in § 5.3. However, Georgiou et al. (1989) showed that the numerically determined extrudate swell ratios for the planar extrudate swell problem are not sensitive to small variations of the exponents of the singular solution, so that we would not expect the extrudate swell ratio at $\gamma =0$ calculated by matching the far-field solution to $y = 1+O(r^2)$ to be significantly different from the values calculated by the above-cited authors.

Figure 5. Finest mesh ($N_1=N_2=N_3=618$). Log-normal plot of extrudate swell ratios $\eta (\infty )$ computed at values $1$, $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

4.2.2. Shape of the free surface

In figures 6–8 we plot the converged free surface function $y=\eta (x)$, given by the linear combination (3.25), as well as its first and second derivatives, over the range of dimensionless surface tensions considered in this article. The computations were done on our finest mesh ($M=1854$). Although the choice of coefficients $\{\beta _i\}_{i=1}^n$ in (3.25) is constrained by the requirement that $\eta (0)=1$, the other coefficients $\{\alpha _i, \gamma _i, \epsilon _{0,i}, \epsilon _{\infty,i}\}$ appearing in the basis functions $\eta _i(x)$ are determined, as was explained in § 3.2.1, so as to minimize the residual of the normal stress balance (3.18) in the interval $[x_0, x_\infty ]$. It follows that

(4.2)\begin{equation} \eta''(0) ={-}2\sum_{i=1}^n \beta_i\alpha_i (\epsilon_{0,i} - \epsilon_{\infty,i}) \gamma_i, \end{equation}

will be finite and that there is no expectation that

(4.3)\begin{equation} \eta'(0) = \sum_{i=1}^n \beta_i\alpha_i\epsilon_{0,i}, \end{equation}

will equal $\arctan (\alpha )$.

Figure 6. Finest mesh ($N_1=N_2=N_3=618$). Free surface $y=\eta (x)$ computed at values $1$, $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

Figure 7. Finest mesh ($N_1=N_2=N_3=618$). First derivative $y=\eta '(x)$ computed at values $1$, $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

Figure 8. Finest mesh ($N_1=N_2=N_3=618$). Second derivative $y=\eta ''(x)$ computed at values $1$, $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

Given the inadequacy of the representation (3.25) of the free surface in the neighbourhood of the separation point, we match the correct asymptotic form (2.19) as explained in § 3.2.2 to this far-field solution. In figure 9 we show the results of the matching at $x=x_0$ of (2.19) with (3.25), as computed on our finest mesh. Over this very small interval ($x\in [0,2r_c]$) both the graph of (2.19) and that of (3.25) resemble straight line segments. However, $\kappa (0)$, the curvature at the origin, is infinite.

Figure 9. Finest mesh ($N_1=N_2=N_3=618$). Matching of the inner solution (2.19) to the outer solution $y=\eta (x)$ at $x=x_0$ (see (3.2)) at values $1$, $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

In table 2 we show the values of the parameters $h_1$, $h_2$, $\alpha$ (in radians), $\lambda _1$ and $c_p$ computed on the finest mesh. We compare the values of $\alpha$ shown in the fourth column of the table with the few that can be found elsewhere in the literature in § 4.2.4. As mentioned in the introduction, following an earlier conjecture by Schultz & Gervasio (1990) that the free surface $y=\eta (x)$ near the separation point is of the form (1.1), Salamon et al. (1995) estimated the coefficients to be $a=0$, $b=0.176$, $c=0.0263$ and $n=1.43$ for $\gamma =1$ and $Re=0$. This gives $\lambda _1 = n -1 = 0.43$, not too far from our own computed value of $0.45321$.

Table 2. Values of the parameters $h_1$, $h_2$, $\alpha$, $\lambda _1$ and $c_p$ calculated from (3.28) to (3.32). The far-field solutions were computed on the finest mesh ($M=1854$).

4.2.3. Normal stress on the free surface

Evidence, in the case $\gamma =1$, of the successful solution of the least-squares problem (3.18) is shown in figure 10. As noted in the introductory paragraph of § 4, values of the objective function (3.18) for all surface tensions were always (sometimes significantly) less than $O(10^{-5})$. In figure 11 we show the free surface normal stress values $\gamma \kappa (x)$ computed for $\gamma = 0.001$, $0.01$, $0.1$ and $1$ and $x\in [x_0, x_\infty ]$. Since these computations are not valid in the immediate neighbourhood of the separation point we have matched the normal stress with the asymptotic solution at $x=x_0$ (see (3.29)) and show the results of this matching in figure 12, where the open circles indicate the matching point for each value of $\gamma$. The $O(r^{\lambda _1-1})$ behaviour of the normal stresses very near the singularity is clearly seen in the log–log plots of figure 13. The slope of the triangle hypotenuse is $\lambda _1 -1 = 0.45321- 1 = - 0.54679$ and is the gradient of the line segment corresponding to the $\gamma =1$ data for very small $r$, although as we may note from column 5 of table 2 there is very little difference in the values of the exponent of $r$ as the surface tension decreases from $1$ to $0.001$. Vanishing or infinite surface tensions would require that $\lambda _1=0.5$ (see Michael (1958) and the discussion in §§ 5.1 and 5.2). Since $-1<\lambda _1<1$ the components of viscous stress, although singular at $(0,1)$, therefore remain integrable and the components of the velocity vanish as $r\rightarrow 0$. The data shown in figures 10–13 was computed using the finest mesh $M=1854$.

Figure 10. Finest mesh ($N_1=N_2=N_3=618$). Signed curvature (2.16) and normal stress $\sigma _{nn}$ (2.17) on the free surface for $x\in [x_0, x_\infty ]$ at value $\gamma =1$ of the dimensionless surface tension.

Figure 11. Finest mesh ($N_1=N_2=N_3=618$). Normal stress $\sigma _{nn}=\gamma \kappa (x)$ on the free surface for $x\in [x_0, x_\infty ]$ at values $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

Figure 12. Finest mesh ($N_1=N_2=N_3=618$). Matching of the inner normal stress $\sigma _{nn} = \gamma \kappa (x)$ on the free surface to the outer normal stress at $x=x_0$ (see (3.2)) at values $1$, $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

Figure 13. Finest mesh ($N_1=N_2=N_3=618$). Log–log plot of the inner normal stress $\sigma _{nn} = \gamma \kappa (r)$ on the free surface at values $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

4.2.4. Separation angles $\alpha$

In figure 14 we present the values of the separation angles $\alpha -180^{\circ }$ for different values of $\gamma$ calculated on meshes with $M$ varying from $624$ to $1854$. The results of our finest mesh calculations have already been presented in radians in the fourth column of table 2 and upon conversion to degrees give $\alpha -180^{\circ }$ values of $9.36^{\circ }$, $11.01^{\circ }$, $11.16^{\circ }$ and $11.15^{\circ }$ when $\gamma =1$, $0.1$, $0.01$ and $0.001$, respectively. As mentioned in the introduction, the data from the experimental results reported in Batchelor et al. (1973), Nickell et al. (1974), Tanner (1986) and Tanner et al. (1985) showed separation angles to be between $189^{\circ }$ and $194^{\circ }$. The high-resolution finite element computations by Salamon et al. (1995) allowed for an estimate of $\alpha =(180+ \textrm {arctan}(0.176))\approx 189.98^{\circ }$ when $\gamma =1$. The values of $\alpha$ that we have computed therefore fall within the interval arising from previous experimental and numerical papers, although there is, admittedly, a paucity of such results available in the literature.

Figure 14. Meshes with $N_1=N_2=N_3$ varying from $208$ to $618$. Separation angle $\alpha - 180^{\circ }$ computed at values $1$, $0.1$, $0.01$ and $0.001$ of the dimensionless surface tension $\gamma$.

5. Discussion and conclusions

We now draw some conclusions from the analysis of § 2.2 of this paper and the numerical results that we have presented in § 4.2. In particular, in the light of our analysis and numerical results, we try to state concisely what can be said about the planar free surface assumed by Michael (1958), and the forms of the free surface both when the surface tension $\gamma =0$ and when $\gamma >0$.

5.1. Planar free surface

A planar free surface, as assumed by Michael (1958), cannot be predicted by the analysis of § 2.2 since, assuming $h_1=h_2=c_p/\gamma =0$, (2.35), (2.38)–(2.42), (2.44) and (2.50) lead to the conclusion that $\alpha ={\rm \pi}$, $\lambda _1=1/2$ and $\displaystyle \lambda _2=3/2$, which contradicts the assumptions that $p_1=\lambda _1$ and $p_1=\lambda _2-\lambda _1$ (and hence that $\lambda _2=2\lambda _1$). When $\gamma \not =0$, a planar free surface can be predicted by the theory of §§ 2.2.4–2.2.6 only when $c_p/\gamma =0$ (that is, when either $c_p=0$ (as assumed by Michael) or $\gamma \rightarrow \infty$, corresponding to ‘stick-slip’ flow (infinite surface tension)) and either $p_1\not =\lambda _2-\lambda _1$ or $p_1\not =\lambda _1$.

5.2. Vanishing surface tension

Suppose that $\gamma =0$. Then from (2.50) we have

(5.1)\begin{equation} (3\lambda_1^2+1)f_1'(\alpha) + f_1'''(\alpha) =0, \end{equation}

and this is equivalent to equation (5) of Michael (1958) and, as in that article, leads to the conclusion that $\alpha ={\rm \pi}$ and $\displaystyle \lambda _1=1/2$. As explained in § 2.2.1 we require $c_p=c_p(0)=0$. Furthermore, from (2.53) (with either $p_1\not =\lambda _1$ or $p_1\not =\lambda _2-\lambda _1$, so that $\lambda _2\not =2\lambda _1=1$) it follows that

(5.2)\begin{equation} \frac{(3\lambda_2^2+1)f_2'(\alpha) + f_2'''(\alpha)}{\lambda_2(\lambda_2^2-1)} =0, \end{equation}

and this requires, given the form (2.22) for $f_2(\theta )$ and the satisfaction of the wall conditions (2.31), that either $B_2=D_2=0$ or that $\cos (\lambda _2{\rm \pi} )=0$.

Suppose $B_2=D_2=0$. Then $f_2(\theta )=2C_2\sin (\lambda _2\theta )\sin \theta$ and therefore, since $f_2({\rm \pi} )=0$, we have from (2.35) that $h_1=0$.

Now suppose that $\cos (\lambda _2{\rm \pi} )=0$. Then this leads to the conclusion that

(5.3)\begin{equation} f_2''({\rm \pi})=(\lambda_2^2-1)f_2({\rm \pi}). \end{equation}

If $p_1=\lambda _2-\lambda _1$ then we would have from (2.43), (2.50) and (5.3) that

(5.4)\begin{equation} -2h_1\lambda_1f_1'({\rm \pi})(p_1 + \lambda_1)=0, \end{equation}

so that (since $f_1'({\rm \pi} )\not =0$) either $p_1=-\lambda _1 \Rightarrow \lambda _2=0$ (which is not possible from (2.23)) or $h_1=0$. However, if $p_1\not =\lambda _2-\lambda _1$ it follows from (2.35) that $h_1=0$.

We conclude that $\gamma =0$ leads to a separation angle equal to ${\rm \pi}$ and a surface that is locally described by $y = 1+O(r^2)$.

5.3. Non-zero surface tensions

For $0<\gamma <\infty$ and non-zero $h_1$ and $h_2$, (2.50) and (2.53) lead to the conclusion that $\alpha \not ={\rm \pi}$ and indeed will not even tend to ${\rm \pi}$, however small the value of the surface tension. This is because with both $p_1=\lambda _1$ and $p_1=\lambda _2-\lambda _1$, if $\alpha \rightarrow {\rm \pi}$, $\displaystyle \lambda _1\rightarrow 1/2$ and $\lambda _2\rightarrow 1$ and the term on the right-hand side of (2.53)

(5.5)\begin{equation} \frac{(3\lambda_2^2+1)f_2'(\alpha) + f_2'''(\alpha)}{\lambda_2(\lambda_2^2-1)} \sim 2({B}_2+{D}_2), \end{equation}

becomes unbounded (unless $h_1$ is identically zero, but this means that $\alpha ={\rm \pi}$ from (2.52)). Thus we identify the case of $\alpha ={\rm \pi}$ when $\gamma =0$ to be a singular limit.

The values of $\alpha$ that we have computed for values of $\gamma =0.001$, $0.01$, $0.1$ and $1$, fall well within the interval of values arising from previous experimental and numerical papers. Our analysis in § 2.2 shows that for non-zero surface tensions the normal stress and curvature are unbounded at the point of separation, consistent with the analysis and results of Anderson & Davis (1993), Schultz & Gervasio (1990) and Salamon et al. (1995). The numerical results presented in § 4.2 (where we determined $h_1$, $h_2$, $\alpha$, $\lambda _1$ and $c_p$ by matching the local solution (2.19) and its curvature to the corresponding values of the global solution $y=\eta (x)$ and satisfying (2.42), (2.50) and (2.53)) have pointed precisely to this scenario.